text:
\section{Introduction}
Operating close to the critical point at a phase transition between order (subcritical) and disorder (supercitical) has long been associated with optimal performance of complex systems. Several biological systems, such as gene regulatory networks \citep{balleza2008critical, ramo2006perturbation}, neural networks \citep{tkacik2015thermodynamics, schneidman2006neural_cultures}, collectively behaving cells \citep{halley2009collective_cells, depalo2017collective_cells} or swarms \citep{cavagna2010scale, chate2014insect} have been shown to operate close to a critical point.
 Criticality has been associated with an ability to solve complex tasks \citep{villegas2016noise_subcritical},  flexibility towards changes in the environment and good evolvability \citep{aldana2007robustness_evolvability} in complex living systems \citep{kauffman1993origins}.
 All these properties provide an adaptive advantage in natural environments, leading to the assumption that evolutionary dynamics push living systems close to the critical regime.
 \par
 \begin{figure}
\centering
\includegraphics[width=.95\columnwidth]{Figures/Methods_figure_cut_smaller.png}
\caption{
Snapshot in time of population dynamics and schematic representation of the control network. \textbf{A:}
Environment with 50 organisms (red circles with trails) foraging for food (green dots).
\textbf{B:} A network with $12$ neurons. Four sensory, four hidden and four motor neurons. All potentially allowed edges are displayed. The exact topology and edge weights are subject to evolution by the EA.
}
\label{fig:methods}
\end{figure}

 On the other hand, it has been suggested that the ubiquitous presence of noise in nature pushes living systems into a more robust subcritical regime.
 For example, in an evolutionary model of Random Boolean Networks (RBNs) decreasing the system size, making the task less complex, or introducing noise to the system pushes the optimal regime further into the subcritical range \citep{villegas2016noise_subcritical}. Similarly \citet{pauli_ramo2007_critical_information_propagation_subcritical_noise} observed, that information propagation is maximized in critical RBNs, however the optimal regime shifts slightly into the subcritical regime under the presence of noise.
 A related phenomenon has been observed on neuromorphic chips, where simpler tasks lead to optimal behavior in the subcritical regime, whereas harder tasks require progressively more critical dynamics \citep{cramer_control_2020}.
 Finally, for some applications, the combination of systems at different distances to criticality is shown to lead to optimal results \citep{zierenberg2020tailored}. The supercitical state has been universally observed to perform poorly \citep{villegas2016noise_subcritical, kauffman1993origins}.

The benefits of criticality for the evolvability of living systems have been associated with the genotype-phenotype coupling.
Specifically, it has been shown  \citep{dejong2006evolutionary_unified_book} that a tight genotype-phenotype coupling leads to optimal evolvability.
Due to this coupling, the dynamical regime has an impact on the properties of the fitness landscape. In an RBN model, the super- and subcritical regimes were shown to disturb the genotype-phenotype coupling \citep{kauffman1993origins} and lead to either very rugged or overly flat fitness landscapes.
A rugged fitness landscape means that the evolutionary dynamics are just a random search -- inefficient in high dimensions~\citep{kauffman1987waiting_times_double}, whereas
a very flat landscape strongly dampens the optimization process. Both phenomena lead to a complexity catastrophe, where an increase of system size leads to a failure to discover satisfying solutions with evolutionary search.
Critical RBNs result in an intermediately rugged fitness landscape which allows for efficient hill climbing search and is less prone to the complexity catastrophe.

We study how the dynamical regime of populations of evolving organisms influences their ability to solve a task. Our investigation is conducted in a simple foraging game of scalable difficulty.
This system allows us to analyze the changes in the dynamical regimes during the evolution.
In addition, we propose a potential answer to the question of which dynamical regime demonstrates the best performance and stability with respect to changes in the environment.

\section{Methods}

\begin{figure*}
\centering
\includegraphics[width=\textwidth]{Figures/heat_capacity_compressed.png}\vspace{-.8em}
\caption{The dynamical regime of a network can be calculated by scaling the inverse temperature with a factor $c_\beta$ to measure its heat capacity and measuring its distance to the corresponding peak.
Heat capacity (Eq.~\ref{eq:Heat_Cap}) of the Ising  networks (Figure~\ref{fig:methods}) for 50 initially subcritical ($\beta_{\mathrm{init}} = 10$, blue), critical ($\beta_{\mathrm{init}} = 1$, green), and supercritical ($\beta_{\mathrm{init}} = 0.1$, red) organisms as a function of $c_\beta$.
For each organism it reaches the maximum (marked by a dot) at individual values $c_\beta = c_\beta^\mathrm{crit}$.
Dynamical regime $\delta = \log (c_\beta^\mathrm{crit})\approx - \log (\beta_{\mathrm{init}}).$
The displayed populations are unevolved and the resulting dynamical regimes closely correspond to their respective $\beta_{\mathrm{init}}$.}
\label{fig:heat_cap}
\end{figure*}

\fs{We investigate a 2D environment where organisms controlled by individual neural networks forage for food}.
Each organism gains energy by eating food particles and consumes energy by moving. The organisms eat the food particles as soon as they run over them.
We can increase the difficulty of the task by requiring the organism's velocity to be below a certain threshold when running over food in order to be able to consume it.
The fitness of an organism is determined by its average energy throughout its lifetime, and an evolutionary algorithm (EA) optimizes the network that controls the organism to maximize its fitness.

\subsection{Organism}

The organisms in our model are controlled by an Ising neural network (INN) that has been previously used by \citet{aguilera2017ising_neural_network} as well as \citet{khajehabdollahi2020sinas_paper}.
The Ising network consists of $N$ neurons that can be in one of two states $ s_i \in \{-1, 1 \}, i = 1, \dots N$. There are sensory neurons that only receive input from the sensors and motor neurons that control the agent. All other neurons are hidden units.
The connectivity of the network is described by the weight matrix $J \in [-2, 2]^{N \times N}$
and the adjacency matrix $A \in \{0,1\}^{N \times N}$, as shown in Figure~\ref{fig:methods}B.
The adjacency matrix is such that no direct connections between sensor and motor neurons exist.
Following the Ising model, each network state has an associated energy
\begin{eqnarray}
E(s_1,\ldots,s_N)= -\sum\limits_{\{i,j \mid A_{ij}=1\}} J_{i,j} s_i s_j
\label{eq:ising_net_energy}
\end{eqnarray}
where $\{i, j \mid A_{ij}=1\}$ is the set of connections.
The network stochastically minimizes the energy by following Glauber dynamics:
At each network iteration, all non-sensor neurons are updated in a random order and the state of neuron $i$ changes from $s_i$ to $-s_i$ with probability:
\begin{eqnarray}
p_i &=& \frac{1}{1+e^{\beta \cdot \Delta E_i}},
    \label{eq:spin_flip}\\
\Delta E_i &=& E (s_1,\mydots,s_i,\mydots,s_N) - E (s_1,\mydots,-s_i,\mydots,s_N), \nonumber
\end{eqnarray}
where $\beta$ is the inverse temperature of the network ($\beta = 1 / (T \cdot k_B)$,  $k_B$ is the Boltzmann constant which we set to one and omit for simplicity) and $\Delta E_i$ is the change in the energy of the network that is caused by the spin-flip of the $i_{th}$ neuron (changing its state $s_i$ to $-s_i$).
The energy change $\Delta E_i$ is determined by the connectivity matrix $J$ and the states of neighboring neurons.
A negative energy change leads to a greater likelihood of a flip. The parameter $\beta$ controls the likelihood of energetically unfavorable flips.
A larger $\beta$ leads to deterministic network behavior dominated by the connectivity, whereas a smaller $\beta$ leads to more random behavior.
In principle, many iterations would be required to converge to  equilibrium. For practical reasons, we perform a fixed number of iterations (10) to process a new sensor input and create motor commands.

An organism has four input neurons that receive information about the angle $\theta_\mathrm{food}$ and distance $d_\mathrm{food}$ from the closest food particle as well as its own velocity $v$ and energy $E$.
Moreover, each organism has four output neurons that control linear and rotational acceleration (2 neurons each) and $N_h$ hidden neurons (Figure~\ref{fig:methods}B).
For most simulations we take $N_h = 4$, to test scalability of learning we also consider a network with $N_h = 20$.
At the beginning of each simulation, an organism is provided with an amount of initial energy $E^\mathrm{init}=2$.
Movement reduces energy and consuming food particles increases it.
We consider two versions of this environment:  In the \emph{simple task}
organisms consume food when passing over it.
In the \emph{hard task} organisms have to slow down and almost stop to be able to consume food.
Unless stated differently, a simulation lasts for a \emph{lifetime} of $t=2000$ time steps after which the evolutionary algorithm (EA) is applied, and the task is simple.
50 INN-controlled organisms are placed in a $2$D environment with periodic boundaries and ever-respawning food particles, conserved to a value of 100 (Figure \ref{fig:methods}A).

\subsection{Evolutionary algorithm}
The evolutionary algorithm applied to the INNs consists of a combination of elitism, mutation, and mating.
At the end of the 2D simulation described above, the fitness of each organism is defined as their mean energy throughout their lifetime.
Subsequently, the 20 fittest organisms are selected for reproduction.

The next generation of organisms is produced by applying a combination of copying, mutation, and mating procedures to the selected individuals.
The copying algorithm transfers some of the fittest individuals into the next generation unchanged.
The mutation algorithm adds or deletes edges in $A$ (connections not present in Figure~\ref{fig:methods}B cannot be added), re-samples a random edge weight in $J$ from a uniform distribution $\mathcal{U}(-2, 2)$, and perturbs the inverse temperature with multiplicative noise $\beta' = \beta \cdot \Delta \beta$ where $\Delta \beta \sim \mathcal{N}(1,\,0.02)$.
Finally, the  mating algorithm takes a weighted average of the connectivity $J$ and inverse temperature $\beta$ from two parents to produce an offspring.
In most of our simulations, the EA iterates for 4000 generations.

\begin{figure*}\centering
\begin{tabular}{ c c c }
\textbf{A.} Simple Task & \textbf{B.} Simple Task & \textbf{C.} Hard Task\\
$12$ Neurons & $28$ Neurons & $12$ Neurons \\
 \includegraphics[height=3.5cm]{Figures/Gen_vs_fitness_panels_no_scatter/Gens_vs_fitness_panel_1.pdf} & \includegraphics[height=3.5cm]{Figures/Gen_vs_fitness_panels_no_scatter/Gens_vs_fitness_panel_3.pdf} & \includegraphics[height=3.5cm]{Figures/Gen_vs_fitness_panels_no_scatter/Gens_vs_fitness_panel_5.pdf} \\
 \includegraphics[height=4.0338cm]{Figures/Gen_vs_fitness_panels_no_scatter/Gens_vs_fitness_panel_2.pdf} & \includegraphics[height=4.0338cm]{Figures/Gen_vs_fitness_panels_no_scatter/Gens_vs_fitness_panel_4.pdf} & \includegraphics[height=4.0338cm]{Figures/Gen_vs_fitness_panels_no_scatter/Gens_vs_fitness_panel_6.pdf}

\end{tabular}\vspace{-.5em}
\caption{Critically initialized populations can be successfully evolved in different circumstances, whereas for subcritically initiated populations, a harder task or an increased system size can lead to the breakdown of evolutionary dynamics.
For each panel 10 initially critical ($\beta = 1$, green, top row) or initially subcritical ($\beta = 0.1$, blue, bottom row) populations evolve for 4000 generations.
The dashed horizontal line at fitness = 2 in the subcritical panels correspond to the organisms' initial energy. It can be seen in the 12 neuron, hard task panel (bottom right) that the network is unable to achieve a fitness above its starting value.}
\label{fig:fitness_plots_vs_generation}
\end{figure*}

\subsection{Defining the dynamical regime of an organism}

\renewcommand{\thealgorithm}{}
\begin{algorithm}[tbp]
\caption{}
\label{alg:model}
\begin{algorithmic}[1]

\State {\color{ourred} \textbf{Evolutionary Algorithm}}
\Comment{\textit{Evolves $J$, $A$, and $\beta$}}
\For{$\text{generation} = 1 \textbf{ to } \text{Generations}$}
\State {\color{ourred} \textbf{Foraging Game in $2$D environment}}
\hspace{\algorithmicindent}\For{$t = 1 \textbf{ to } \text{Time Steps}$}
\State $\text{Update Sensor Neurons}(\text{organism})$
\State {\color{ourred} \textbf{Update INN}}
\hspace{\algorithmicindent}\hspace{\algorithmicindent}\For{$\text{iter} = 1 \textbf{ to } \text{Network Iterations}$}
\hspace{\algorithmicindent}\hspace{\algorithmicindent}\hspace{\algorithmicindent}\For{$\text{non-sensor neuron} \textbf{ in } \text{INN}$}
\State $\text{Potential spin-flip(neuron)}$
\EndFor
\EndFor
\State $\text{Read Motor Neurons}(\text{organism})$
\State $\text{Move in 2D environment}(\text{organism})$
\EndFor
\State $\text{Evolve}(\text{population of organisms})$
\State $\text{Reset 2D environment}$
\EndFor
\end{algorithmic}
\end{algorithm}

We use the heat capacity from statistical physics to derive a measure of an organism's dynamical regime (sub-, super-, critical).
In our finite system, we estimate the putative divergence point by changing the inverse temperature $\beta$ multiplying it with a scaling constant $c_\beta$.
This change of temperature influences how likely the state of the neurons will flip (Eq.~\ref{eq:spin_flip}), and thus change the equilibrium distribution of energies $E$ (Eq.~\ref{eq:ising_net_energy}). Now we search for a $c_\beta$ that delivers the maximal value of the heat capacity $C_H(c_{\beta})$, defined as:
\begin{eqnarray}
    \label{eq:Heat_Cap}

    C_H(c_{\beta}) = \frac{1}{T^2} \textrm{Var}(E) = c_{\beta}^2  \beta^2  \textrm{Var}(E).
\end{eqnarray}
We define the $c_\beta^{\mathrm{crit}} = \underset{c_{\beta}}{\mathrm{argmax}}\, C_H(c_{\beta}) $. An analogous procedure was used in~\citep{tkacik2015thermodynamics}.

We define the distance of the network from the critical point by the logarithm of the scaling factor required to bring the network to criticality, $\delta = \log(c_\beta^{\mathrm{crit}})$.
For unevolved organisms (first generation, Figure \ref{fig:heat_cap}), the relationship between $\beta_{\mathrm{init}}$ and $\delta$ can be approximated by:
\begin{eqnarray}
\delta \approx - \log \beta_{\mathrm{init}}, \: \beta_{\mathrm{init}} \in [-0.1, 10].
\label{eq:beta_delta_approx}
\end{eqnarray}
On a technical side, for the stability of the numerical procedure for finding the equilibrium distribution to estimate $C_H(c_{\beta})$ we initiate the motor and hidden neurons in the state $(s_1,\ldots,s_N)$ delivering minimum of $E(s_1,\ldots,s_N)$, while keeping sensory inputs fixed to the observed values.
During the evolution, the distance from the critical point can (and will) change.

\section{Results}

\subsection{Convergence of evolution}
Populations of different initial states follow distinct evolutionary strategies but are all able to solve the standard foraging task.
We observe evolution for 4000 generations, concentrating on the populations initiated at the supercritical ($\beta = 10$, $\delta \approx 1$), critical  ($\beta = 1$, $\delta \approx 0$) and subcritical ($\beta = 0.1$, $\delta \approx -1$) regimes.
Critical populations begin to rapidly gain fitness from the first generation in all EA realizations (Figure~\ref{fig:fitness_plots_vs_generation}A).
The gradual and stable increase of fitness of the initially critical population suggests that successful hill climbing of the fitness landscape is taking place.
In contrast, for subcritical populations, fitness mainly evolves via random jumps, and only half of the simulations reach the same fitness as the critical populations after 4000 generations (Figure~\ref{fig:fitness_plots_vs_generation}A).
Such fitness dynamics indicate a random search strategy, which often leads to a population getting trapped in local maxima for extended periods of time.
Confirming the previous observations by~\citep{khajehabdollahi2020sinas_paper}, we see that supercritical populations, after an initial random period, follow the same path as the critical ones.

\begin{figure*}\centering
\begin{tabular}{ c c }
\textbf{A.} Simple Task & \textbf{B.} Hard Task\\
 \includegraphics[height=1.8in]{Figures/delta_vs_generations_all_in_one_georgs_final_simple_task.pdf} & \includegraphics[height=1.8in]{Figures/delta_vs_generations_all_in_one_georg_final_hard_task.pdf}
\end{tabular}\vspace{-.5em}
\caption{Changes in the distance to criticality over the course of evolution.
Populations initiated at various distances to criticality ($\delta$ between -1.5 and 1) and evolved on a \textbf{A:} simple task and \textbf{B:} hard task.
The color indicates their fitness at generation 4000. Populations with $\delta < -1$ remain at fitness 2 for the hard task, signifying no evolutionary progress.
}
\label{fig:delta_vs_generation}
\end{figure*}

\begin{figure}
\centering
\includegraphics[width= 0.95\columnwidth]{Figures/simple_vs_hard_task.pdf}\vspace{-.5em}
\caption{The dynamical regime of the initially critical population in the harder task remains closer to the critical regime, than in the simple task throughout evolution.
Main plot: For each condition 10 populations are initiated in the critical regime end evolve for 4000 generations (dark - hard task, light - simple task).
Inset: histograms of dynamical regimes at generation 4000 of 75 populations per task, vertical lines denote means, $p < 10^{-13}$. Difference between means in inset and main panel are due to the small sample size for the main panel.
}
\label{fig:simple_vs_hard_task}
\end{figure}

For successfully evolvable populations, moderate changes in the complexity of the control network should not destroy the ability of the EA to reach a good fitness.

We test the differences in evolvability for initially critical and subcritical populations by changing the size of the network from 12 to 28 neurons.
Same as for smaller networks, the initially critical populations rapidly evolve in all realizations.
By generation 4000 they even reach a slightly larger fitness than populations of smaller networks that had evolved for the same amount of time (Figure~\ref{fig:fitness_plots_vs_generation}B).
At the same time, subcritially-initialized populations do not reach even half of their normal fitness.
We observe the same difference between the dynamical regimes when we increase the task's complexity requiring organisms to slow down almost to zero velocity in order to consume food particles (Figure~\ref{fig:fitness_plots_vs_generation}C).
In this harder task, the evolved populations' maximal fitness is expected to be lower than for the simple task.
For the initially critical populations, we still observe the same hill climbing dynamics. However, the initially subcritical populations stay at an energy level of exactly two.
This signifies that they do not use the originally supplied energy for moving and remain static throughout all 4000 generations, trapped in a local optima.

Overall, we see that although in simple tasks all populations can converge to approximately the same fitness, there exists a significant difference between the initially subcritical and initially critical populations.
Namely, a convergence of the EA for critical populations resembles hill climbing. It is stable (all populations follow very similar fitness growth) and behaves similarly for larger networks or more complex tasks.
At the same time, for subcritical populations, the EA resembles random search, which is stochastic and fails to find solutions in high-dimensional cases or for more complex tasks.

\subsection{Evolution of the dynamical regime}
\begin{figure}
\centering
\includegraphics[width= 0.95\columnwidth]{Figures/Generalizability.png}
\vspace{-.5em}
\caption{Populations initialized at criticality are always able to find solutions whose behaviors generalize in time beyond their training condition, whereas subcritical populations often overfit to the same conditions.
{\bf A:} 54 populations of each type (triangles - initially critical, circles - initially subcritical).
After 4000 generations of evolution all critical populations reach a high fitness (indicated by color) and nearly perfect (with one exception) generalizability $\gamma_t$ (Eq.~\ref{eq:generalizability_param}).
The initially subcritical organisms split into badly generalizing cluster 2 (19 populations), generalizable cluster 3 (28 populations), and 7 populations that we did not assign to any cluster. This split is independent of their attained fitness.
{\bf B:} Energy and velocity as a function of time for representative examples of the organisms from cluster 2 and clusters 1,3 (marked in panel A). Organisms in cluster 2 reach their maximal fitness (sometimes higher than in clusters 1 and 3) at the end of their lives.
The dashed orange line denotes the training lifetime of 2000 time steps.}
\label{fig:generalizability}
\end{figure}

\fs{Next, we investigate how evolution changes the state of the populations.}
To do so, we select a wide range of initial dynamical regimes ($\delta \in [-1.5, 1]$) and examine how the dynamics of populations initialized in each of these regimes change throughout evolution.

Regardless of their initial dynamical regime, all populations converge to the subcritical regime, however, with different distances from the critical point (Figure~\ref{fig:delta_vs_generation}).
We also observe that strongly subcritical populations ($\delta < -1 $) generally achieve lower fitness levels in the simple task and are unable to solve the hard task.

Both near-critical and supercritical populations rapidly change their dynamical regimes and by generation 4000 reach an intermediately subcritical regime $\delta^*$.
Strongly subcritical populations with $\delta\ll\delta^* $ remain at their initial regimes, demonstrating a lack of evolutionary mobility, whereas subcritical populations initialized at different $0 \geq \delta \geq \delta^*$  can still approach $\delta^*$ (Figure~\ref{fig:delta_vs_generation}).

Task complexity changes the dynamical regime found in the evolutionary limit. The hard task requires a smaller distance from criticality.
We check the evolution of the dynamical regime in both simple and hard tasks (Figure~\ref{fig:simple_vs_hard_task}).
We utilize the observation that all populations with an initial regime $\delta > \delta^*$ converge to similar values. Thus, we consider only initially critical populations.
We obtain the distribution of dynamical states by considering 75 independent runs of evolution in both tasks after 4000 generations.
The harder task results in a smaller distance from the critical regime (large $\delta$), and the difference is significant ($p<10^{-13}$).
The difference is observed throughout the evolution (Figure~\ref{fig:simple_vs_hard_task}, the main panel based on ten populations).

\fs{Initiating the evolution close to the critical regime is important for unknown task complexity.}
We observe that the dynamical regime never changes towards larger values, but the subcritical convergence point can be at different distances from criticality.
Thus, only starting at the critical point guarantees that the optimal dynamical state can be reached by evolution.

\subsection{Generalizability}
\fs{ For successful biological systems robustness towards environmental change is the paramount feature, therefore, it should be used to distinguish the successfully evolved artificial organisms.}
We propose a simple measure to investigate how the model behaves outside of its explicit training conditions.
Specifically, for a population trained for the organism's lifetime $t_{\mathrm{train}}$ we define generalizability as the speed of growth of the average fitness  if the organism's lifetime is extended to $t_{\mathrm{extend}}$ . Formally:
\begin{eqnarray}
\gamma_t = \frac{\langle E_{t_\mathrm{train}} \rangle / t_\mathrm{train}}{\langle E_{t_\mathrm{extend}} \rangle / t_\mathrm{extend}}.
\label{eq:generalizability_param}
\end{eqnarray}
The \emph{stable generalizability}, $\gamma_t = 1$ corresponds to linear growth whereas sublinear behavior $\gamma_t \ll 1$  indicates possible overfitting to the particular organism's lifetime $t_{\mathrm{train}}$.

\fs{We consider initially critical ($\delta\approx0$) and initially subcritical ($\delta\approx-1$) populations evolved for 4000 generations and then test their performance for an extended lifetime of 50\,000 time steps (instead of the 2000 in training).}
 As reported in previous sections, the critical populations converge to $\delta \approx -0.55$, and they all have a similar fitness after training.
Interestingly, when increasing the organisms' lifetime, the fitness of the critical population continues to grow linearly, signifying almost perfect generalizability.
About half of the subcritical populations reach the same fitness level.
However, the subcritical populations split up into two clusters: one with generalizability close to 0 and another with generalizability close to 1 (Figure~\ref{fig:generalizability}A).
Surprisingly, there is no difference in fitness between these two clusters.
When we look more precisely at the individual organisms from the two clusters, we discover that the strategy of the organisms in the non-generalizable cluster is to increase the velocity permanently until the end of their training lifetime (Figure~\ref{fig:generalizability}B).
However, moving with such a high velocity is not compatible with the energy influx from feeding, and they break down shortly after the end of their training lifetime, this demonstrates that these organisms overfit the training conditions.
For all generalizable populations (clusters 1 and 3), the velocity fluctuates around a stable level, and the energy of the organisms grows linearly.

\fs{Overall, the initialization in the critical regime results in almost perfect generalizability of evolved populations, whereas initially strongly subcritical popuations risk overfit their training conditions.}

\subsection{Effect of genetic perturbations on the fitness}
\fs{Next, we examine the stability of the evolved organisms to genetic perturbations.}
We apply genetic perturbations of different magnitudes to the evolved organisms of initially critical $\beta_{\mathrm{init}}=1$ and subcritical $\beta_{\mathrm{init}}$ populations.
We perturb all weights of the connectivity matrix by randomly adding or subtracting a number $f_\mathrm{pert}$ and then evaluate the fitness of the resulting organism.
We find that fitness rapidly declines with perturbation magnitude for both populations, however the subcritical ones decline faster (Figure~\ref{fig:phenoty_genotype}A).
We evaluate the fitness decline by the slope of an exponential function fitted to the fitness.
For the hard task it is $\alpha = -2.26$ and for the simple task it is $\alpha = -5.03$, which is more than double the decline rate, indicating a much higher sensitivity of subcritical systems to perturbations.

\fs{The EA is a source of constant genetic perturbations that are necessary in the beginning of evolution but can become detrimental later}.
We consider the individual effect of the evolutionary operators (copy, mutate, and mate) on the resulting fitness of the organisms.
The variability of fitness for copying simply reflects the natural variability in community fitness rankings and organism behaviour.
However, both mating and mutation in fully evolved subcritical populations typically results in a fitness close to 2 -- signifying totally unfit organisms (Figure~\ref{fig:phenoty_genotype}B).
At the same time, initially critical organisms retain diverse fitness values after mutation and mating, some being close to the optimum. This indicates that the originally critical populations retain their evolvability as opposed to the rigid search performed by strongly subcritical populations.

\begin{figure}\centering
\includegraphics[width= 0.95\columnwidth]{Figures/genotype_phenotype_mapping_A.pdf}
\includegraphics[width= 0.95\columnwidth]{Figures/density_mutations_B.pdf}
\caption{Initially critical populations show larger genotypic stability than initially subcritical ones.
\textbf{A:} The phenotype (fitness) as a function of the genotypic perturbation (changes in connectivity) for two initial conditions: $\beta_{\mathrm{init}} = 1$ and $\beta_{\mathrm{init}} = 10$.
Dashed black lines indicate exponential fit with exponent -2.26 (critical) and -5.03 (subcritical).
\textbf{B:} The histograms of the fitness values for nearly fully evolve agents (born between generation 3500 and 4000), categorized according to the last evolutionary operator (copy, mutate, or mate) that was applied to them.
The $\beta_{\mathrm{init}} = 10$ agents are much less likely to remain fit when their genotype is changed by either mutation or mating.
}
\label{fig:phenoty_genotype}
\end{figure}

\begin{comment}
\begin{figure}\centering
\begin{tabular}{@{}c@{}}
A\hfill\vadjust{} \\[-.1em]
\includegraphics[width= 0.95\columnwidth]{Figures/genotype_phenotype_mapping.png}\\
B\hfill\vadjust{} \\[-.1em]
\includegraphics[width=0.95\columnwidth]{Figures/density_mutations.png}
\end{tabular}\vspace{-.5em}
\caption{Initially critical populations show larger genotypic stability than initially subcritical ones.
\textbf{A:} The phenotype (fitness) as a function of the genotypic perturbation (changes in connectivity) for two initial conditions: $\beta_{\mathrm{init}} = 1$ and $\beta_{\mathrm{init}} = 10$.
Dashed black lines indicate exponential fit with exponent -2.26 (critical) and -5.03 (subcritical).
\textbf{B:} The histograms of the fitness values for nearly fully evolve agents (born between generation 3500 and 4000), categorized according to the last evolutionary operator (copy, mutate, or mate) that was applied to them.
The $\beta_{\mathrm{init}} = 10$ agents are much less likely to remain fit when their genotype is changed by either mutation or mating.
}
\label{fig:phenoty_genotype}
\end{figure}
\end{comment}

\section{Discussion}
\fs{We demonstrate that in various scenarios evolving populations of agents converge to a moderately subcritical state with the resulting deviation from criticality depending on the tasks difficulty. }
This might appear to be a contradiction to the previous studies, suggesting that operating close to criticality is optimal for natural systems \citep{mora2011biological_poised_at_critic, munoz2018colloquium, roli2018dynamical}.
However, a recent body of research showed that for simple tasks, operating at some distance to criticality might be an optimal solution for the sensitivity/stability tradeoff~\citep{hidalgo2014information,tomen2014marginally,villegas2016noise_subcritical,cramer_control_2020}.

We observe that the distance from criticality affects an agent's ability to solve complex tasks and to robustly evolve generalizable behavior.
Specifically, we observe that slightly subcritical populations are evolvable for different complexities of the control network and task, whereas strongly subcritical populations fail in both cases.
Given that these properties are crucial for adaption in natural environments, we propose that living systems operate in the subcritical regime in close proximity to the critical point.
Moreover, we show that the optimal regime moves closer to criticality as we increase the task difficulty, which suggests that the optimal distance from criticality varies.
Those findings are confirmed by \citet{cramer_control_2020} as well as \citet{villegas2016noise_subcritical}, who showed that the optimal distance from criticality in the sub-critical regime decreases for higher task complexity or larger system size.
We further observe that populations can only become more subcritical during evolution and fail to decrease their distance to criticality even when this would have eventually led to superior behavior.
As it is \emph{a priori} unknown which distance from criticality will be optimal when evolving for a new task, starting at the critical point could be the only way for the evolutionary process to descend to the optimal regime.
However, in the long run this would require some sub-populations to always maintain closeness to criticality.
How this can be achieved for neuronal networks is a subject of vivid research (for a review, see~\citep{zeraati2020self,buendia2020feedback,kinouchi2020mechanisms}) and for the embodied Ising agents it remains open for further investigation.
Maintaining evolvability in simpler systems was, for instance, observed by switching between different rough energy landscapes \citep{wang2019evolving}.
The inhomogeneity of the environment and coevolution can also contribute to the preservation of the critical regime~\citep{hidalgo2014information}.
 Overall, the maintenance of evolvability throughout evolution is an important question beyond the embodied Ising agents studied here.

Our results extend and partly revise the earlier findings of \citet{khajehabdollahi2020sinas_paper}, that reported a superior evolvability of critical populations and an approximate convergence to criticality during evolution.
We confirm that the critical regime allows reliable evolvability, additionally, we extend our understanding by considering a set of tasks and architectures in the model.
However, our more precise procedure to infer the dynamical regime and the fine sampling of initial conditions uncover additional intricate dynamics.
Namely, that critical populations of Ising-agents converge to the subcritical regime, and the distance to criticality depends on the task complexity.
We also propose a new way to investigate capabilities of the resulting organisms by defining generalizability and genetic stability measures.
Both measures reveal the benefits of staying close to the critical state beyond a simple fitness comparison.

\subsection*{Acknowledgments}

This work was supported by a Sofja Kovalevskaja Award
from the Alexander von Humboldt Foundation, endowed by the Federal Ministry of Education and Research. We acknowledge the  support from the BMBF through the T\"ubingen AI Center (FKZ: 01IS18039B).

\footnotesize
\bibliographystyle{apalike}