Evolutionary Implications of Self-Assembling Cybernetic Materials with Collective Problem-Solving Intelligence at Multiple Scales

Evolution works with an active rather than a passive substrate, i.e., with biological cells with agendas of their own [1]. Thus, at every stage of development during morphogenesis, collective decisions are made at vastly different length- and time scales within an organism, guiding the formation of the mature phenotype. We aim to model exactly this process via the neural cellular automata (NCAs) described in Section 3.1 and employ evolutionary algorithms (EAs) to evolve the parameters of such NCAs so the latter perform well on a target morphogenesis task, see Section 3.2. Without loss of generality, we consider an Nx×Ny=8×8 Czech flag pattern (as a more complex version of the classic French flag problem of morphogenesis [115,116]) as the target pattern for our in silico morphogenesis experiments, with a fixed number of Nj=64 cells in total, NG=3 distinct cell types (colored blue, white, and red, respectively) and NH=1 hidden state, which renders the dimension of the NCA cell state NC=4. We use a square grid of cells with N=8 neighbors per cell and with fixed boundary conditions (see Appendix C for details).

Starting from a genotype xj defined in Equation (3), we perform a number of tD=25 developmental steps per morphogenesis experiment to “grow” a phenotype pj, described by Equation (4), based on which the fitness score r(pj) is evaluated following Equation (5) (see Figure 2B for an illustration of this process). During this entire process, we limit the magnitude of the numerical cell state values ci(tk) at all time steps tk to the interval lc=[−3,3], and, analogously, limit the magnitude of the proposed actions ai(t) of each uni-cellular agent to the interval la=[−1,1]. This is achieved by clipping the numerical values of ci(tk+1) after a cell state update described by Equation (1), and the ANN outputs ai(t) to the respective limits lc and la. The noise level ξc defined in Equation (1) is counted in units of the action limits max(la) and is thus sampled from a Gaussian distribution with zero mean and standard deviation ξc independently for each of the NC=4 cell state elements, thus affecting the cell state updates during development; the actual numerical values for the hyperparameters above turned out to be well suited for the problem at hand, especially to reasonably compare and discuss simulation results for the means of this contribution but are not crucial for the more general aspects of the evolutionary implications of multi-scale intelligence discussed here.

To study the effects of different types of decision-making machinery within a cell, we utilize two different architectures for the NCA artificial neural networks (ANNs), a Feedforward (FF) and a recurrent ANN inspired by gene regulatory networks [117,118,119,120] (RGRNs). (The terminology FF and RGRN stems from the respective agents’ Feedforward and Recurrent Gene Regulatory Network ANN controller layers (see Appendix A for details)). Notably, the RGRN-agent architecture augments cells with an internal memory that is independent of their states in the NCA and thus can not be accessed by the cells’ neighbors. To balance the length of the structural genome x(S) and functional genome x(F) defined in Equation (3), the two ANN architectures, FF and RGRN, are chosen such that the number of parameters NFF=192 and NRGRN=164 is roughly the same as the number of initial cell states Nj×NC=64×4=256. Thus, the ANNs utilized here—and detailed in Table A1 of Appendix A—are tiny compared to Ref. [101].

For each experiment of evolving the parameters of an NCA, i.e., for each independent run of the EA, we typically utilize a population X of NP=96 individuals and a maximum number of NM=2000 generations. As the EAs ultimate fitness criterion, we consider the average Fj=〈r(pj)〉NE of NE=8 statistically independent fitness scores r(pj) of corresponding morphogenesis simulations starting from an individual j’s genotype xj and resulting in a corresponding phenotype pj after tD developmental steps; the developmental program described via Equation (1) is imperfect due to the developmental noise applied to the cell state updates and can thus lead to different, noise-induced phenotypic realizations. Typical values used here for the corresponding reward factors defined in Equation (5) are rT=0.25 and rS=0.5. We consider the problem solved if a fitness of Fj=max(nj(C))=Nj=64 is reached, but since we reward individuals to maintain the target pattern over time (via rT), the maximum possible fitness score after tD developmental time steps is max(rj(pj))=70.25 in this example. Further details about the hyper-parameters of the EA and afforded computational resources can be found in Appendix D.

We aim in this contribution to investigate the evolutionary implications of biologically inspired multi-scale competency architectures [1,94]. Thus, we compare two qualitatively different evolutionary processes, both with the objective of morphogenetic pattern formation but whose genomes either (i) directly encode the phenotypic features of a two-dimensional target pattern (cf. Figure 1A), or (ii) encode the cellular competencies of a multi-scale architecture that gives rise to the same phenotypic features (cf. Figure 1C). Notably, different definitions of direct and indirect encodings in multi-agent systems have been used in the literature [54]. Here, we specifically distinguish between structural parameters xj(S)={ci(0)}j in the search space that directly encode features of the phenotype, i.e., specific initial cell types gi(0)≈g^i and functional parameters xj(F)=θj that indirectly, or rather functionally, encode the target pattern by parametrizing the intercellular communication and intracellular information processing competencies of the NCA that facilitate the self-orchestrated pattern formation process.

If no ANN at all were present in our model, i.e., θj={}, and in the absence of noise ξc=0, we would re-establish a direct mapping between the genotype and phenotype as ci(0)=ci(tD), and thus a direct encoding of the target cell type pattern could be achieved gi(0)=gi(tD). However, by default, we allow each cell to successively regulate its own cell state towards a target homeostatic value via an iterative perception–action cycle defined by Equation (1) and, moreover, to communicate in that way with neighboring cells. More specifically, each cell updates its cell state solely based on its own and the states of its adjacent neighbors, which, in turn, update their states based on their respective local environment. We explicitly avoid direct environmental feedback to the cells’ perception (such as an individual or collective reward signal) but fully restrict the NCA to intercellular communication (via cell state updates) and intracellular information processing. These uni-cellular agents thus exhibit a certain level of problem-solving competencies that can be utilized for the challenge at hand, in our case, for a collective system-level objective of forming a specific two-dimensional target pattern [95,101,102].

With the explicit partitioning of the genome into a structural part, i.e., xj(S), and a functional part, i.e., xj(F), we can study the effect of direct vs. multi-scale, or the competency-driven encoding of phenotypic traits in the process of evolution, and, moreover, quantitatively tackle the question of whether competent parts affect the process of evolution and evolvability. In any case, the initial cell state pattern is given by the structural part of the genome. Thus, in the absence of noise and without any active functional part in the genome, the set of initial cell states directly represents the final pattern, while otherwise, cell states can either be modified passively by noise in the system or actively through actions by the cells during the developmental stage. Thus, we employ CMA-ES [103] to either evolve the (i) structural, or both (ii) the structural and functional parts of the genome of an NCA simultaneously with the shared objective of self-assembling an 8×8 Czech-flag pattern in tD=25 developmental time steps in the presence of noise ξc=0.25 (cf. Section 3.2 and Section 4.1 for details). More explicitly, in case (i), we restrict the cell state update of the NCA by disabling all cell actions ai(tk)→ai*(tk)=0 but formally keep the functional part of the genome in the evolutionary search. In turn, we allow the NCA in case (ii) to afford both the structural and functional parts of the genome, thus giving the evolutionary process the opportunity to prioritize one over the other. We thus bias the evolutionary process in case (i) to effectively search the space of direct phenotypic encodings, while keeping the search space dimensions balanced in both cases. The results of this experiment are presented in Figure 3.

We can see in Figure 3 that both the evolution of the (i) direct and (ii) multi-scale encoding schemes of the target pattern can be achieved with the presented framework, and a fitness threshold of Fj=64 is reached after ≈300–600 generations, thus solving the problem. However, depending on the encoding scheme (i) or (ii), we can identify clear qualitative differences in the strategy and the “efficiency” of the evolutionary process, i.e., how many generations it takes to reach a certain fitness threshold and eventually converge (cf. top and bottom panel of Figure 3, respectively). The respective fitness score of the direct case (i) grows steadily and almost monotonically over successive generations until the threshold of Fj=64 is reached after 668 generations for that particular run, and the EA converges at a maximum fitness of maxFj(i)=70.25 after 942 generations (see Section 4.1 for details on the threshold fitness values). In contrast, the evolutionary process of the multi-scale case (ii) undergoes significant leaps as reflected by the corresponding fitness score, which can increase rapidly if a suitable innovation, i.e., a favorable crossover or mutation event in the functional parameters θj, occurs; the initial standard deviation of the fitness of the entire population is significantly larger compared to the direct case (i), yet the threshold fitness of Fj=64 is reached in 428 generations, and the EA converges after 679 generations (although at a lower maximum fitness of maxFj(ii)=69).

The results presented in Section 4.2 are based on selected evolutionary optimization runs that are representative of related experiments with similar parameterizations. However, one should keep in mind that such results are always susceptible to chance in initial conditions or mutations in the EA but also to developmental noise; moreover, the hyperparameters of the evolutionary search or even the specific ANN architectures can influence the evolvability of such NCA systems. Thus, we present in Section 4.3 below a more statistically significant analysis of the evolutionary implications of direct and multi-scale encodings under various conditions of the cellular agents’ competency levels and the developmental noise.

Our separation of the genotype xj into a structural xj(S)={ci(0)}j and into a functional part xj(F)=θj moreover allows us to extract the structural (or genotypic) fitness along an entire evolutionary history: we define the structural fitness as the fitness score r(pj*) of a phenotype pj* with evolved structural genes {ci(0)}j but with disabled agency ai(tk)→ai*(tk)=0. Notably, in the direct case (i), we have pj=pj*, which is illustrated in Figure 1A and reflected in the top panel of Figure 3; the structural fitness of the multi-scale case (ii) is explicitly visualized in the bottom panel of Figure 3. In the latter case, the structural fitness remains essentially detached from the phenotypic fitness pj*≪pj during the entire evolutionary history (which also explains the convergence to a suboptimal maximal fitness level of max(Fj)=69 in this particular NCA solution, as the final Czech flag pattern first needs to be assembled from the corresponding imperfect initial cell configurations xj(S)). This all suggests that in contrast to (i), the EA in (ii) can make the most use of exploring the functional part of the genome, i.e., the space of behavior-shaping signaling and information processing [1], and, in turn, that the mere presence of competent parts drastically changes the search space accessible to evolution [3]; to show this explicitly, we present in Appendix E an illustration of the evolution of the morphogenesis process.

Interestingly, we still observe a slow but steady increase in the structural fitness in the long term in case (ii), owed to the small additional reward signal rT reinforcing the cellular agents to maintain the target pattern over time. This can most efficiently be achieved if the agent starts from a perfect set of initial cell types, representing a particular sub-space in the parameter space that might not necessarily be easily accessible to the EA at all stages during the evolutionary search. However, we would like to stress that such a slow transfer of problem-specific competencies from an agential, highly adaptive functional part xj(F) to a rather rigid structural part xj(S) of the genome could be a manifestation of the Baldwin effect [14].

While remaining neutral with respect to the system level fitness score, this competency transfer seems to affect the entire population presented in Figure 3 as reflected by the successively increasing population-averaged structural fitness score. Notably, and as detailed in Appendix F, we identified an associated decrease in the robustness against increasingly noisy cell state updates of the corresponding solutions with larger structural fitness. This suggests a reduction in uni-cellular competencies and might relate to the “paradox of robustness” discussed in Refs. [51,52,53,121,122,123,124,125,126,127]. Through a computational lens, such a competency transfer would also allow, as soon as the structural part of the genome is reliable enough, to repurpose the system’s competency to adapt to other independent tasks, and thus may facilitate the, in biology, ubiquitous effects of adaptability and polycomputing in related systems [128].

This all illustrates that an agential material [1,94], or more precisely, a substrate composed of competent parts, can have significant effects on the process of evolution and evolvability, especially for morphogenesis tasks. We thus conclude that, if competent parts are available, evolution prefers exploiting competency over direct encoding—if the environment requires competency at all (see discussion in Section 4.3). This leads to the conclusion that “competency at the lowest level greatly affects evolution and evolvability at the system level”.

Here, we investigate the effects of varying different levels of competency at the cellular level of a multi-scale competency architecture on the evolutionary process of morphogenesis. More specifically, we introduce the decision-making probability (I) PD, which constrains the ability of each cell individually to perform cell state updates in the environment: PD defines the probability at which a proposed cell state update of each individual cell in the NCA is executed (or otherwise omitted). Thus, varying the decision-making probability from PD=0 to PD=1 smoothly transitions the system’s behavior from a direct encoding scheme without competency to an increasingly reliable multi-scale competency architecture (cf. Figure 3).

Another, somewhat hidden, level of competency we already discussed in Section 3.1 is each cell’s ANN architecture: An RGRN agent with internal memory can acquire and execute tasks differently than a simpler FF agent without any feedback connections except for its cell state ci(tk). Comparing the evolutionary implications of such functionally different ANN architectures is, however, not trivial, and is thus kept to a minimum here.

However, we parameterize both FF and RGRN agents such that their controller part of the ANNs (cf. Figure 1C, Section 3.1, and Appendix A) are (II) stacks of R redundant copies of the same controller ANN, each copy with its own set of parameters, which take the same pre-processed aggregated sensor embedding as input, and whose individual outputs are averaged into a single action-output of a cell. Inspired by redundancy in error-correcting codes [105,106], we thus allow cells with higher values of this redundancy number R, i.e., with many alternative routes through the controller part of the ANN, to—in principle—integrate environmental signals more generally compared to R=1, thus affecting the cells competency.

While scaling from PD=0 to PD=1 smoothly increases a cell’s competency to reliably regulate its cell state, increasing R enhances the computational capacities of the uni-cellular agents. Henceforward, we interpret (I) PD and (II) R as two competency levels in our system, which we can vary (I) continuously and (II) discretely.

Analogous to Section 4.1 and Section 4.2, we thus utilize CMA-ES to evolve the genotypic parameters of an NCA to self-assemble the 8×8 Czech flag pattern under different conditions (I–II), and expose the cells to different noise-levels (III) ξc during cell state updates defined in Equation (1).

In Figure 4A,B, we present the corresponding fitness scores of a maximum of 2000 generations of CMA-ES for different noise levels ξc∈[0,0.5], averaged over different values of the decision-making probability PD∈{0,12.5%,25%,50%,100%} for both FFagents and RGRN agents. Moreover, for each realization of ξc and PD, we utilize experiments with different redundancy numbers R∈{1,2,4,8,16} and employ 15 statistically independent EA runs for each parameter combination ξc, PD, and R, and thus arrive at 75 statistically (and functionally, with respect to an agent’s ANN architecture) independent fitness trajectories per (PD, ξc) combination; see Section 4.1 and Appendix D for more details on the EA parameters. In Figure 4C, we present the average number of generations it takes to solve the problem (to reach a fitness threshold of Fj=64) for each combination of PD and ξc, aggregated over the agents’ ANN architectures, FF or RGRN, and the respective redundancy numbers R for 15 statistically independent EA runs each; in Figure 4D, we present the data from Figure 4C but separately for both ANN architectures.

We observe in Figure 4 that, depending on these two parameters PD and ξc, for no- or very low noise levels ξc≈0, the evolutionary search is most efficient, i.e., finds the solution in the fewest number of generations, on average, for low values of the competency level PD≈0. Thus, in these situations, direct encoding (achieved via PD=0) seems to be preferable to competency-driven encodings with PD>0 (as indicated by the bottom red arrow in Figure 4C); this is partly owed to the specific definition of the cell types gi(tk) given by Equation (2), making a noiseless search very simple for the EA. However, for more realistic, noise conditions ξc>0, the situation changes drastically. With increasing the noise level, the evolutionary efficiency of NCAs with higher competency levels is significantly greater compared to those with low competency levels, especially for the direct encoding scheme (as indicated by the green arrows in Figure 4C); for noise levels of ξc=0.375 and 0.5, the EA does not even find solutions for the direct encoding case with PD=0 in 2000 generations, as cell state updates become increasingly necessary to counteract the noise in the system. There is a clear trend of increasing the evolutionary efficiency in our in silico morphogenesis experiments by increasing the competency level for increasingly difficult environments with high noise levels.

Thus, we conclude that scaling competency has a strong effect on the process of evolution, and in realistic situations (with moderate to high noise), competency may greatly improve the evolutionary efficiency and evolvability of collective self-regulative systems.

It might be noteworthy that for evolving the 8×8 Czech flag pattern, essentially no qualitative difference in the evolutionary efficiency between FF agents and RGRN agents with the given number of parameters was observed. Also, the evolutionary implications of utilizing a number of R>1 redundant copies within the controller ANNs of the cells of an NCA is much less pronounced, compared to the results depicted in Figure 4 as can be seen in Figure A8 of Appendix H. However, for more advanced problems such as assembling a 9×9 smiley-face pattern (see Appendix G), RGRN agents seem to outperform simpler FF agents significantly in terms of evolutionary efficiency. Moreover, a larger redundancy number of R≥4 is required by the evolutionary process to more efficiently evolve the functional parameters of an NCA compared to a direct encoding scheme, hinting at a capacity bottleneck of the deployed ANNs.

A careful analysis of the results shown in Figure 4 reveals that the largest competency level of PD=1 does not result in the highest evolutionary efficiency for any presented noise level. On the contrary, populations with slightly lower competency levels of PD=0.5 or even PD=0.25 perform best at noise levels ξc∈{0.25,0.375,0.5} and 0.125, respectively (as indicated by the green and red arrow ends in Figure 4C). In fact, cells with an initially random genome (comprising the ANN and initial cell state parameters) that are forced to make “uninformed”, i.e., initially random, decisions at every time step can interfere with the performance of the EA, as even initially perfect cell state configurations will be destroyed during such a randomized developmental stage. We suspect that this leads to corresponding delays in the evolutionary search compared to situations where populations can better rely on the structural part of the genome. Indeed, populations with “overconfident” actions can be trapped in local optima for many generations at all stages of the EA, which, in our system, may only be resolved by very specific but random mutations of the functional part of the genome (as we show later through Figure 5 in Section 4.4). This is reflected in Figure 4A,B by the large deviations in the average fitness trajectories for large PD values.

The insights from the above lead to the questions of whether there is a “natural” or optimal competency level, with respect to the decision-making probability PD, or whether a mutable competency level can be utilized by the evolutionary process to improve the efficiency of guiding a population towards high-fitness regions in the parameter space. Thus, we include the decision-making probability as an additional competency gene xj(C) into the NCA genome xj→xj=xj(S)∪xj(F)∪xj(C), cf. Equation (3), and we perform in silico morphogenesis evolution experiments of the 8×8 Czech flag pattern for different noise levels ξc, analogous to Section 4.3. We analogously limit the numerical range of the competency gene xj(C) to the interval [−3,3], and extract the corresponding decision-making probability via PD,j=12(tanh(xj(C))+1). Notably, for the experiments shown in this sub-section, we use L2 regularization on the genotypic parameters xj=(xj,1,⋯,xj,Nx) through subtracting rL2×∑i=1Nxxj,i2 from the fitness score defined in Equation (5), with rL2=0.01 (the L2 regularization applied to xj does not introduce a bias between the minimal PD,j=0 and maximal PD,j=1 competency levels, as the L2 regularization is applied to xj(C), not PD,j; both PD,j and the L2 regularization are symmetric with respect to the sign of xj(C)).

In Figure 5A,B, we present the evolved competency level for different noise levels after fitness thresholds of 64 and 70 are crossed, respectively, for 10 independent lineages per noise level for an RGRN architecture. The problem is considered solved at a fitness of 64, but since we reward the NCAs to maintain the target pattern over time via rT in Equation (5), a higher maximal fitness score of 70.25 can be reached after tD developmental steps for a sufficiently long evolution. Thus, we here relate Figure 5A to the evolutionary stage of having achieved the process of morphogenesis, and Figure 5B of having achieved morphostasis. For both cases, we essentially see two strategies emerging (see also Figure 5C–E): (i) one, where competency is maximized very early during the evolutionary process that then remains near the maximally possible value of PD=1, and (ii) a hybrid strategy where a significantly lower competency level is assumed that still allows to solve the problem.

Notably, strategy (i) is predominantly pursued at high noise levels, where large cell state fluctuations in the environment favor informed actions by the cellular agents. In contrast, the second strategy (ii) emerges more frequently in lineages evolved at low noise levels where, especially at very low noise levels ξc≈0, most of the evolutionary processes result in solutions that avoid competency altogether, and a direct encoding scheme (PD = 0) is evolved. Intermediate competency levels evolve in the corresponding intermediate noise regime. Following the trend of evolving morphogenesis (by crossing a fitness score of 64) to morphostatsis (by converging to the maximal fitness value of ≈70) in Figure 5A through B, we see that the two strategies, (i) and (ii), “sharpen” during the course of the evolutionary process such that PD predominantly converges to the minimally or maximally possible values of 0 and 1, depending on the environmental conditions.

We also illustrate the evolved competency level of the particular lineage at all noise levels in Figure 5A,B, at which the respective fitness threshold is crossed in the minimal and maximal number of generations (and on average) amongst all 10 independent lineages per noise level. This clearly reveals that evolutionary processes that follow a more direct encoding strategy (ii) can evolve the problem at hand efficiently—if this is permitted by the developmental noise. However, when increasing the noise level, the evolutionary process can afford to evolve—or put differently, increasingly relies on evolving—the multi-cellular intelligence of the NCA to perform morphogenesis and morphostasis, thus following a third strategy (iii) that integrates both strategies (i) and (ii) in a non-trivial way. We observe in Figure 5A that the most efficient strategy for evolving morphogenesis seems indeed to be such a hybrid approach (iii), where a minimally necessary competency level is utilized at a specific noise level such that the corresponding evolutionary process can, again, be very efficient in solving the task.

Moreover, this also holds for the stage where morphostasis is reached, cf. Figure 5B: lineages that efficiently evolved to solve morphogenesis in our experiments also (typically) evolve to solve morphostasis efficiently. To emphasize this, we present in Figure 5C–E the “temporal dynamics” of the population-wise highest fitness and the corresponding competency level per generation for all lineages at selected noise levels ξc={0,0.125,0.25}; we also present for all corresponding lineages that have been evolved at these selected noise levels the genotypic competency level PD,j against the corresponding phenotypic fitness scores rj, and we find an apparent yet non-trivial relation between these two quantities: typically, an initial rise in fitness rj in early generations is associated with a decline in PD,j which is more pronounced at lower noise levels. For intermediate noise levels 0<ξc≪1, we find that PD,j often assumes a minimum (i.e., a minimally required yet finite competency level) when the evolutionary process reaches a fitness level of ≈64. We suspect that this allows the evolving morphogenetic process to establish good starting configurations based on changes in the structural genome, which can most efficiently be performed at a minimal(ly necessary) competency level given a certain developmental noise level in the environment. However, the competency is then quickly pulled towards a maximum level of PD,j=1 when the EA converges at a maximum fitness score of ≈70, at the morphostasis stage. For large noise levels, e.g., ξc=0.25 as depicted in Figure 5E, the competency level rises with the corresponding fitness score in a much more monotonic way, emphasizing the necessity of the corresponding NCAs to utilize the cellular competency to solve the problem already at an early stage of the evolutionary process.

Curiously, we also see lineages that settle at the highest possible competency levels throughout their evolutionary history, even in conditions without noise as can be seen in Figure 5C: here, an initial “frozen accident” may cause an entire lineage to maintain high competency levels due to a lack of diversity in the corresponding gene, although this is not even necessary to solve the task. However, these high competency levels early on during the evolutionary process can cause the population to stagnate at sub-optimal regions in the parameter space for many generations if the corresponding policy of the cells is sub-optimal but rigid to strategy changes via small mutations in the genome. The population seems “trapped” until a favorable mutation or crossover event occurs in the functional part of the genome of an individual that guides the entire population towards higher fitness scores, eventually solving the problem. We suspect that this is also the reason for the lower evolutionary efficiency of the “most competent” configurations (with PD=1) compared to the slightly less competent cases (with PD=0.5) of the experiments depicted in Figure 4 [129].

Thus, we conclude, that if the evolutionary process can afford to evolve its own competency level, there seems to be a trade-off—during the entire course of the evolutionary process—between “going direct” or “going competent”, depending on the developmental noise. Moreover, the randomly initialized starting conditions may favor either direct or multi-scale encoding strategies, which may not only affect the “final” competency level that the evolutionary process converges to but can also greatly influence the efficiency of the evolutionary process itself. In general, the most efficient strategy for evolving morphogenesis tasks seems to be a non-trivial trade-off between finding a suitable initial cell state configuration that then allows the competency-based self-assembly of the target pattern to “kick in” and solve the task efficiently.