The Structure of Ontogenies in a Model Protocell

Biological individuals are a special class of physical systems that counter the universal trend towards disintegration. In any given moment, a physical system has a structural configuration—the spatial arrangement of components from which it is constituted. The states and locations of these components unfold dynamically, and whereas most configurations tend towards a uniform equilibrium, biological individuals are a unique subclass of physical systems that persist as individuals. The theory of autopoiesis argues this viability comes from their closure of production; as a result of their intrinsic dynamics and material exchange with the environment, biological individuals produce and distribute the materials needed to stabilize themselves [18].

An individual and its environment form a coupled system, with the state of the individual in part determining the state of the environment, and the state of the environment in part determining the state of the individual. Thus, from the perspective of the individual, the exact state of the environment is unknown; instead, the environment appears as a probability distribution over possible perturbations to the individual's state. These perturbations, together with the individual's intrinsic dynamics, determine the individual's subsequent states. Thus, a perturbation can have one of three consequences: The individual can be unaffected (a robust transition), it can be changed to a different viable configuration (a plastic transition), or it can cross into the set of nonviable configurations and disintegrate (a destructive transition) [2].

An individual that remains viable for any length of time experiences a sequence of perturbations that induce a corresponding sequence of configurational changes. In this article, we refer to an unbroken trajectory through the set of viable configurations as an ontogeny; when the trajectory crosses the boundary of viability into a nonviable region, closure of production is broken and the ontogeny ends. Different sequences of perturbations have the potential to induce different ontogenies. If the set of possible perturbations is known, one can in principle map the entire network of possible transitions that an individual can undergo, by exposing copies of the same individual to all sequences of perturbations and recording the resulting sequences of configurational change.

Viability is often defined as a biological individual's capacity to sustain its function. This is traditionally treated as a binary notion: Either the individual is currently viable or it is not. For example, W. R. Ashby postulated a set of essential variables that need to be maintained within a viable region for a system to persist [3]. Autopoiesis also proposes a binary notion: Either a system persists through closure of production, or it disintegrates. More recently, the binary notion of viability has been extended to a more graded notion [5, 10]. This extension focuses not on whether a system is inside or outside of a viable region, but on its temporal distance to the boundary. In this article, the viability of any configuration is measured as the average number of environmental perturbations it is from disintegration, weighted by the probability of those perturbations.

Adaptivity is associated with increases in an individual's viability resulting from interactions with its environment [4, 10]. Some have suggested that adaptivity requires an explicit regulatory mechanism [8, 10]. However, as an operationally defined notion, adaptivity should be neutral with respect to the mechanism by which it is achieved. Instead, by adopting an operational definition of adaptivity, the mechanisms underlying adaptive change become an interesting target of inquiry, rather than presumed a priori. This generalization allows for examples of adaptivity that arise from more implicit, systemic processes rather than only those with explicit, dedicated mechanisms.

Developing an account of ontogeny, viability, and adaptivity requires theoretical models to demonstrate how these notions can be quantified. Simple models, in particular, allow for unambiguous definitions and computations. There is a long history of building such models. These include particle-based models of basic autopoiesis [24], particle-based models capable of self-replication [15], models of the origins of protocells [19], models exploring the interrelations of autopoiesis and behavior [11], reaction-diffusion systems that support emergent individuals [25], lattice-based models of metabolism-boundary co-construction [21], emergence of organization in lambda calculus [12], emergence of organization in the Game of Life [7], and models with metabolism, self-replication, and membrane synthesis [13]. However, despite the abundance of models, most served merely as demonstrations for a particular concept and were not analyzed further.

In this article, the relationship between viability, ontogeny, and adaptivity is investigated in the context of a model protocell that we recently proposed [1]. After reviewing the model (Section 2) and describing an environment in which we place it (Section 3), the article is organized as follows. Section 4 builds upon a framework for exploring ontogenies as a network structure [6]. Applying this methodology reveals a rich complexity of ontogenic structure, which we characterize through a combination of graph-theoretic, morphology-based, and statistical measures. In Section 5, viability is measured for all configurations across the ontogenic network, and increases in the measure allow us to identify trends of adaptive change. In Section 6, a second ontogenic network is derived from a new environment; comparing it with the original ontogenic network demonstrates a more general understanding of how ontogenies unfold through the space of viable configurations. Finally, in Section 7, the mechanism of a particular adaptive transition is analyzed in detail.

For a model to fully address the emergence of individuality, the modeled individual must emerge from lower-level processes rather than being hardcoded into the model. The model itself can only define primitives that do not have any of the individual-level properties of interest. To accomplish this within a spatial chemical system, individuals need to exhibit metabolism-boundary co-construction [2]. In other words, they need to synthesize their own components via a metabolism, and separate themselves from an environment by constructing and maintaining an emergent boundary. While such emergence has been explored in systems with abstract physics, such as the Game of Life [6], the model described here moves towards a more realistic chemistry.

The continuous-time spatial model of molecular concentration dynamics considered here was first introduced in [1]. Since it serves as an exemplar system for the analytical methods proposed in this work, we only offer a brief overview of this model. The physics of this model include the diffusion, repulsion, chemical reactions, and decay of four molecular species: membrane (M), autocatalyst (A), food (F), and water (W). The chemical reactions are such that A is produced from an autocatalytic reaction between A and F, while M molecules are produced by a reaction in which A catalyzes the conversion of F to M. Both A and M also decay at a constant rate. Considered individually, molecular concentrations diffuse across a two-dimensional lattice at a constant rate. However, the presence of repulsion between molecular types breaks the traditional symmetries. Based on the behavior of phospholipid bilayers, the model implements anisotropic repulsion between membrane molecules M and both A and W. This requires the introduction of an additional state variable, θ, which defines the orientation of M at every lattice site and behaves with its own dynamic of alignment with neighboring membrane orientations. The combination of alignment dynamics and anisotropic repulsion supports the self-assembly of membrane sheets. Here, we utilize a 40 × 40 toroidal lattice. Each lattice point has the five state variables (M, A, F, W, and θ), resulting in an 8000-dimensional coupled dynamical system. For specific details of the model's implementation, including the forms of the reaction, diffusion, and repulsion equations, and the relevant parameter values, refer to [2].

The model displays two distinct types of equilibrium points: the uniform equilibrium point with zero concentrations of both A and M (Figure 1a), and stable inhomogeneities in which positive concentrations of A and M are maintained (Figure 1b). A configuration is considered stable when the average temporal derivatives of A and M over the entire lattice satisfy the stabilization criteria (Equation 2 in the Appendix). Even though these configurations are equilibrium points, they are still chemically active, with positive diffusion and reaction rates. It was demonstrated that, in the latter class of equilibrium points, M was necessary to contain A at concentrations high enough that A could continuously construct both A and M faster than their molecular decay. Thus, these configurations exhibit metabolism-boundary co-construction [2].

Throughout this work, we study the ontogeny of a particular configuration labeled SC (stable configuration), shown as the rightmost snapshot of Figure 1b, and on the right side of Figure 2. Due to θ's initialization and subsequent dynamics, SC has broken symmetries across both its horizontal and vertical axes (with the horizontal broken symmetry more visibly apparent than the vertical broken symmetry).

From the individual's point of view, an environment is a probability distribution over a set of perturbations. We first consider an environment ε₁, which consists of perturbations that increase the concentration of either autocatalyst (A) or membrane (M) in the local neighborhoods of nine distinct focal points. These locations are specified relative to the given configuration according to an algorithm that guarantees placement on the configuration (Figure 2). The increase in concentration is determined by a Gaussian function of the form , where |x_ij − x_f| is the Euclidean distance (in numbers of lattice sites) from the focal point x_f to the given cell x_ij, α is the magnitude of the function, and a constant σ² = 2.0. Four different magnitudes of α = [0.5, 1, 1.5, 2] are used. In total, the environment ε₁ consists of 72 possible perturbations (9 locations × 2 molecular species × 4 magnitudes), each of which occurs with a uniform probability.

Once a perturbation is applied, it instantaneously displaces the system in state space. The system dynamics then unfold towards a limit set, resulting in one of three classes of possible outcomes: the uniform state, the same viable configuration, or a different viable configuration. Examples of these three possibilities are illustrated in Figure 3.

To determine the identity of the resulting configuration, C, a comparison metric is used within the 3200-dimensional state space of M and A molecular concentrations (40 × 40 lattice sites for two molecule types). The structural distance from each previously observed configuration, C′, is found by taking the summed absolute difference between the configurations' states. If this value is sufficiently low, then we consider the configurations equivalent. The comparison metric is defined in Equation 3 in the Appendix.

Typical environments are the source of repeated perturbations that induce a sequence of changes to an individual's structure. A single trajectory through the set of viable configurations is an ontogeny. Different perturbations have the potential to induce different plastic transitions, resulting in different ontogenic trajectories. The structure of all possible ontogenies defines an ontogenic network, in which the viable configurations and death state constitute the network's nodes, and each environmental perturbation is a directed edge (see the Appendix for extended definition). Given a specific individual's configuration in a specific environment, the full ontogenic network is obtained by exhaustively characterizing the configuration's response to every environmental perturbation. The process is then repeated for all subsequent configurations until closure is achieved (i.e., every transition results in either a previously characterized configuration or death) [6].

The full ontogenic network formed by SC's repeated exposure to environment ε₁ is a multigraph with a set of reachable configurations as its nodes, each of which is the source for 72 directed edges. Following each perturbation, the system was given sufficient time to stabilize before the next perturbation was applied. The network was generated by a breadth-first search that exposed all configurations to ε₁ as they were discovered. As the frontier of the search expanded from SC, it uncovered an exponentially increasing number of configurations. Due to computational limitations,¹ the search was terminated at a uniform depth of 17 edges from SC, resulting in a total of 268 discovered configurations. The combination of a breadth-first search with truncation at a uniform depth fully characterizes the local search space around SC. Of these configurations, 114 are at a depth of 17, and their outgoing transitions were left unsearched. These configurations will be termed the unsearched frontier, and are excluded from this section's analysis. This leaves 154 configurations whose transitions were fully characterized.

As a first step towards characterizing the structure of this network, we focus on the relationships between the 154 searched configurations. We reduce the full ontogenic network by first combining edges according to equivalence classes of transition, and second by considering the death state with all destructive transitions separately from the rest of the network. The resulting reduced ontogenic network (ON₁) is shown in Figure 4. SC's initial asymmetry is propagated through all configurations in ON₁. In theory, mirrored initial conditions can yield ontogenies identical to ON₁, but with mirrored configurations.

Remarkably, this pairing of a simple configuration and environment generates an extensive ontogenic network, rich with features. The system's many stable configurations (represented in Figure 4 by the morphological structure of A and M concentrations) are shown as nodes, while the death state is shown by the surrounding black ellipse. The directed edges in ON₁ reflect the probability for a transition to occur given a random perturbation from the environment, with robust transitions shown as green self-loops, plastic transitions shown as directed purple edges between configurations, and destructive transitions shown in gray.

A statistical analysis of ON₁'s searched configurations reveals that they vary widely in their response to perturbations from ε₁. It is interesting to note that all configurations have nonzero probabilities of both destruction and survival in this environment. For each configuration, one can calculate its robustness (the probability of undergoing a robust transition, Equation 6), its plasticity (the probability of undergoing a plastic transition, Equation 7), and its fragility (the probability of undergoing a destructive transition, Equation 8). (Equations are found in the Appendix.) Specifically, for the configurations of ON₁, robustness occurs in the range [0.0, 0.79], plasticity occurs in the range [0.0, 0.82], and fragility occurs in the range [0.10, 0.96]. The outdegree distribution of the network reflects the number of different plastic options available for viable configurations. This distribution has a range of [0, 9] with a mean of 3.49, which reflects the fact that the average viable configuration can undergo plastic transitions to three or four other configurations following a perturbation. Further, 30% of plastic transitions are bidirectional—a number that indicates the network has many more bidirectional edges than found in a random graph. Indeed, this hypothesis is supported by a p-value ≪0.001 when comparing ON₁ with an ensemble of 1000 random graphs with the same number of nodes and degree distribution.

The graph-theoretic structure of the viable configurations in ON₁ reflects several interesting characteristics of ontogenies starting at configuration SC in environment ε₁. First, the number of viable configurations increases exponentially with minimum path length from SC. Second, ON₁ has several distinct clusters as highlighted by InfoMap network community detection [22] (Figure 5a). These graph-theoretic clusters reflect sets of configurations for which transitions are more likely to remain within the set than leave it. Third, there are several configurations that function as bottlenecks for the network in the sense that ontogenies must pass through those configurations to reach different areas of the network. These configurations are determined by high values of betweenness centrality [20].

Strongly connected components (SCCs) found in the current ON₁ demonstrate irreversibility, branching, and attractors (Figure 5b). SCCs are sets of configurations in which all configurations are mutually reachable [20]. There are 32 SCCs in ON₁. The presence of multiple SCCs indicates irreversibility—if the system moves from one SCC to another, it cannot return. Branching is illustrated by diverging paths from the SCCs. If an individual exits an SCC along one branch, the alternatives can no longer be explored. Finally, those SCCs with no outgoing connections are attractors in the sense that if an ontogeny enters one of these sets, it is guaranteed to remain there until death.

The analysis of ontogenies is further enriched by considering their morphologies—their unique spatial arrangements of molecular concentrations. Recall that each configuration defines a point in the 3200-dimensional state space of M and A molecular concentrations. Due to their high dimensionality, these configurations are projected onto a 2-dimensional manifold using the IsoMap manifold identification technique [23] including the whole set of other configurations as neighbors. The resulting projection is shown in Figure 6. This projection is indicative of morphological clustering with many points tightly grouped and relatively large distances between the groups. Specific partitions of configurations can be identified using K-means clustering [17]. The resulting eight clusters are shown using different colors in Figure 6. Comparison of the exemplar configurations from each cluster, identified by their central position relative to the cluster's mean, demonstrates how the clusters vary along several qualitative dimensions: the size of the configuration, the shape of the outer membrane, the thickness of the outer membrane, and the number and arrangement of internal membrane structures.

Further graph-theoretic analysis reveals how morphological features change as ontogenies unfold in ON₁ as a function of individual plastic transitions, graph-theoretic path length from SC, and clustering. First, plastic transitions are much more likely to occur between morphologically similar configurations. A Bayesian BEST test [16] between the distribution of all inter-configuration distances using the Euclidean metric and the distribution of configuration distances for those linked by a plastic transition shows that the average difference was 3.25, which fell in a 95% highest-density interval (HDI) of [3.16, 3.34]. Because 0 falls outside the HDI, there is credible evidence that the distributions have different means. Furthermore, the distribution of credible effect sizes is slightly left-skewed, with a mode of 2.18 and a 95% HDI of [2.11, 2.25], which again suggests a credible effect. Second, a reliable pattern of growth occurs in ON₁; the size of morphologies increases with increasing graph-theoretic path length from SC. Finally, there is a strong correspondence between the morphological clusters and the InfoMap clusters previously identified in ON₁ as reflected by a normalized mutual information [9] value of 0.67 between the two clusterings. Therefore, as an ontogeny unfolds and an individual falls into a graph-theoretic cluster, it also tends to maintain its morphological features by remaining within a corresponding morphological cluster.

To proceed in the case of SC in ε₁, we return to the full ontogenetic network and classify each of the 72 perturbations by their probabilities of inducing a robust, destructive, or plastic transition when applied to a random configuration from the set of viable configurations in ON₁. The resulting classification is visualized in Figure 7a according to barycentric coordinates for the three probabilities. Here, each point (arrow or circle) represents one of the 72 perturbations, and the inverse distance between the point and each vertex reflects the associated probability of a transition in that equivalence class; a point directly on a vertex denotes 100% of transitions falling in that category.

Inspection of Figure 7a reveals that perturbations can vary widely in their consequences. Notably, all perturbations but one are destructive to at least one configuration, but none of the 72 perturbations are destructive to all configurations. There is a gradated tradeoff between perturbations' probabilities of inducing a robust transition versus inducing a destructive transition. A small subset of perturbations are associated with a large tendency to induce plastic transitions. These results can be further subdivided according to the location and magnitude of the applied perturbation as shown in Figure 7b. In the first row of the subfigure, three normalized histograms are shown, which categorize perturbations based on outcome, molecular type (membrane in blue, autocatalyst in red), and magnitude (light to dark color). As one would expect, robust transitions are primarily associated with small perturbations, while destructive transitions are associated with large perturbations. Plastic transitions occur more frequently for perturbations to membrane concentrations than to autocatalyst concentrations, regardless of magnitude. The second row of the subfigure illustrates three additional histograms categorizing perturbations based on outcome and location. These diagrams indicate that plastic transitions result more often from perturbations to the center and north east locations, while destructive transitions are slightly biased to perturbations on the north west.

Viability is an individual's ability to maintain itself, and is a consequence of well-matched configurations and environments. While some sequences of perturbations applied to an individual yield long ontogenies, the same individual exposed to different perturbations can result in shorter-lived ontogenies. The viability of a configuration is here defined as the expected life span over all of its possible ontogenies [2]. Two notions of adaptivity follow from this definition. Local adaptivity is the change in viability resulting from a plastic transition; an adaptive transition increases a system's viability, and a maladaptive transition reduces it. Global adaptivity is the tendency for an individual's viability to increase with longer ontogenies; positive global adaptivity is the tendency for viability to increase over time, and negative global adaptivity is the tendency for viability to decrease over time.

The above notion of viability can be quantified by considering the ontogenic network as a transition matrix for a Markov chain. As a consequence of this, death becomes an absorbing state and most configurations are transients. Finding the average number of transitions from each transient state to an absorbing state is a well-studied problem [14], and gives our measure of viability (Equation 9 in the Appendix). The lower bound on this measure is 1, which occurs when all perturbations bring the configuration to death within one step. The existence of an immortal configuration would give an upper bound of infinite viability.

In order to calculate viability for ON₁, all configurations on the unsearched frontier are assigned to transition directly to death. This gives a conservative lower bound for the actual viability of all searched configurations. The analysis of ON₁ reveals a wide range of viabilities, from 1.0588 to 6.3046 with an average of 3.3516, reflecting the fact that the average configuration is 3 random perturbations away from death. Given the viabilities for each configuration, the adaptivity of a transition is calculated from the difference in viability across the transition (Equation 10 in the Appendix). Our analysis identified a surprising abundance of adaptive transitions (Figure 8). Indeed, 40% of the edges are locally adaptive, which means just under half of the plastic transitions produce an increase in viability.

It is important to note that robustness, plasticity, fragility, and viability are all distinct notions, and highlight different features of configurations. This is demonstrated by the fact that configurations identified as extreme examples of these measures are all distinct from one another (Figure 9). The most viable configuration is not the most robust, and the least viable configuration is not the most fragile. Unlike robustness or fragility, which are instantaneous measures, viability is a consequence of how configurations are embedded within the full ontogenic network. The viability of a configuration depends on the structure with which it is connected to other configurations with plastic transitions, and the viabilities of those configurations.

Remarkably, ontogenies from SC also display global adaptivity. A correlation analysis found that a configuration's graph-theoretic path length from SC is positively correlated with its viability (r-value of 0.4622, p < 10⁻⁸). This means that longer-lived ontogenies beginning at SC will generally experience an increase in their viability. Note that viability is a forward-looking measure, not based on how long it took to survive up to a given configuration, but rather on how long it is expected to survive from that configuration. This global trend is configuration specific; each configuration in ON₁ can have a different global adaptivity, and can even be globally maladaptive.

The structure of ontogenies depends on the particular pairing of an individual and its environment. If the individual's initial configuration is changed, we would expect to obtain a different ontogenic network. Likewise, if the environment is changed, we would also expect a different ontogenic network with its own unique properties. Such networks are all instances of a general set of possible ontogenic networks. To gain a broader understanding of this set, it is useful to consider multiple ontogenic networks and compare their structures. In this section, we compare ON₁ with one additional ontogenic network, ON₂ (Figure 10), according to their graph-theoretic properties, morphologies, and the features of interactions between the protocell and its environment.

The second ontogenic network, ON₂, is obtained by exposing SC to a new environment, ε₂. The environments ε₁ and ε₂ are identical in all but one feature: Whereas the locations of ε₁'s perturbations are determined relative to the protocell (as depicted in Figure 2), the locations of ε₂'s perturbations remain fixed in space (at positions denoted in the left of Figure 2). The distinction between environments with spatially relative and those with spatially absolute perturbations is important. Environments with spatially relative perturbations, such as ε₁, track the individual so that it cannot avoid the perturbations. In contrast, environments with spatially fixed perturbations, such as ε₂, allow individuals to adjust their position relative to those perturbations.

The procedures by which the two ontogenic networks were obtained were identical except for two minor differences. First, due to computational limitations, the search of ON₂ was terminated at a depth of 7 transitions from SC rather than 17. Second, the stabilization criteria for two of the configurations in ON₂ was slightly relaxed by reducing the convergence window to 700 time steps.

There are three distinctions that can be drawn between the graph-theoretic structures of ON₁ and ON₂ that reflect key variations in the resulting ontogenies. First, ON₂ contains many more configurations (387) than ON₁ (268), even though the search depth of ON₂ is substantially less than that of ON₁. Given its higher rate of growth with search depth, if ON₂ had been searched to the same depth as ON₁, its size would far exceed the size of ON₁. Second, in ON₂, all branches rooted at SC extend to the unsearched frontier. In contrast, in ON₁ only one branch from SC extends to the unsearched frontier. Third, ON₂ exhibits more irreversibility than ON₁; ON₂ contains 102 strongly connected components, as compared with 32 in ON₁. Since there are many small SCCs in ON₂, when a plastic transition occurs, the protocell is less likely to have the option to return to a prior configuration.

Configurations in ON₂ are on average less robust, less plastic, and more fragile than in ON₁ (Table 1). Configurations in ON₂ are also less viable, and the transitions are less adaptive (only 34% of ON₂'s edges are adaptive, compared to 40% of edges in ON₁). However, these conclusions may be confounded both by the reduced search depth and by the increased size of ON₂'s unsearched frontier.

A comparison of ON₁ and ON₂'s morphologies shows that they explore different regions of morphology space. The 641 total configurations uncovered by either ON₁ or ON₂ can be differentiated by their structures. Indeed, the IsoMap decomposition for all 312 searched configurations (Figure 11) clearly illustrates this separation. The distinction between the configurations in ON₁ and ON₂ is also visible in the comparison of morphological features: We see more configurations with only a single internal membrane structure, there are some morphologies with three internal membrane structures, and many have diagonal membrane structures (one of these is shown in Figure 12). As with ON₁, these configurations cluster according to morphological similarity. Interestingly, there were only 14 configurations shared by the two networks. Of these 14, four include SC and its three immediate progeny, while the others occur within a search depth of 5 in ON₁ and 4 in ON₂. Only six plastic transitions exist within both networks amongst these 14 configurations, and only two transitions are shared that do not involve SC. This comparison suggests that as the search procedure expands beyond the unsearched frontiers, the reachable morphologies will continue to diverge in their features.

The two ontogenic networks also differ in their patterns of morphological change and in the consequences of these changes. For example, ON₁'s configurations show a reliable increase in size, while the size of configurations in ON₂ varies widely. The configurations in ON₂ range in size from slightly smaller than SC to approximately the same size as the largest configuration uncovered in ON₁. The possibility of engulfment is a consequence of increased size in ε₂. Because the perturbations are fixed in space, their effects become increasingly internalized as a configuration grows in size (Figure 12). In contrast, engulfment is impossible in ε₁. As a configuration grows in size in ε₁, the perturbations simply move further away. This ability of an individual to modify the impact of an environmental perturbation on it is a first step towards two-way environment-individual interaction.

This model provides an excellent opportunity to investigate the mechanisms underlying adaptivity. Previous formulations of adaptivity have assumed an explicit regulatory subsystem [3, 10]. This mechanism is designed to monitor the system's internal state relative to its boundary of viability and use this information to bring the system to more viable states. Models that implement this a priori assumption cover only a subset of possible mechanisms of adaptivity. In contrast, this article's model demonstrates an emergent type of adaptivity. Analyzing the processes involved here can provide insights not previously possible.

As a first step towards understanding the mechanisms of emergent adaptivity, we examine a minimally adaptive scenario embedded within ON₁. An adaptive transition requires an environment of at least two perturbations, one that increases viability, and another that reduces it. Taken to its extreme, one perturbation would bring an initial configuration to death, and the other perturbation would bring it to a second, immortal configuration. Multiple instances of this exact scenario are found embedded within ON₁ when environment ε₁ is restricted to only two perturbations.

The chosen example of a minimally adaptive scenario is shown in Figure 13a. The scenario begins with configuration α, which is subjected to the two perturbations. The perturbations are placed at different locations, but otherwise have the same magnitude and are both to the membrane field. On the top branch, α disintegrates, whereas on the bottom branch it undergoes a plastic transition to configuration β. Then β is subjected to the same two perturbations, both of which result in robust transitions. In this minimal context, β can survive all perturbations, whereas α can survive only half of them. The transition from α to β is therefore adaptive.

One way to analyze this scenario is to utilize dynamical systems theory to characterize the system's phase space (Figure 13b). For visualization, the high-dimensional state space is projected onto the top two principal components. This makes trajectories appear to overlap, even though no overlap is possible in the full dynamics. Configurations α, β, and death are equilibrium points, each surrounded by a basin of attraction. The real basins of attraction are not visualizable, so faux basins of attraction are added to represent a hypothetical division of phase space. The two classes of perturbations (dashed red lines) instantaneously displace the system within the state space. Whereas perturbations to α move the system either to death's basin of attraction or to β's basin of attraction, when applied to β they displace the system within the same basin of attraction. Given this analysis, adaptivity is explained by the congruence between the two configurations' basins of attractions and the available perturbations. The position of the equilibrium points and the shapes of their basins are such that the perturbed states fall outside α's basin, yet remain within β's.

A more detailed investigation of mechanism needs to address the specific physico-chemical interactions that take place during this adaptive transition. A spatial analysis is here approached by taking sequential snapshots of transient configurations throughout the minimally adaptive scenario and identifying the critical differences in their morphologies (Figure 13c). First, we look at the adaptive transition from α to β (Figure 13c(1)) and the divergence from α that occurs throughout this transition (Figure 13c(2)). The sequence begins with α, at which there is no difference. Next, the perturbation displaces the membrane field, seen as a circular difference in the top right. Further down the sequence, the membrane concentration spreads around the boundary. While the northeast side of the configuration remains at increased concentrations, the rest of the boundary is slightly lower from its initial concentrations; this trait stabilizes at β. The most obvious difference between stable configuration β and α is a local increase in membrane concentration, which makes β rounder and thickens its boundary (arrow I).

Next, we examine how these morphological differences allow β to survive a perturbation that α does not survive (Figure 13c(3)). The difference between these two sequences shows where the relevant divergences begin (Figure 13c(4)). The sequence begins with the initial difference between α and β. When the perturbations are applied to both configurations, they do not appear in the difference plot, because the perturbations are spatially aligned. Later, the behavior begins to diverge. The obvious difference between their initial structures (arrow I) is not where the fatal divergence begins. Instead, a region on the bottom of α (arrow II), which bulges out and to the right, is the important feature. This growth draws α's membrane downward, reducing the membrane concentration in the southeast boundary (arrow III). This opens up a tear in α's membrane, from which autocatalyst pours out and ultimately brings about disintegration, whereas β's morphology allows it to stay intact and restabilize.

In contrast with previous formulations of adaptivity, we see that no explicit regulatory mechanism is found during this emergent adaptive transition. Adaptivity is the result of distributed processes, and is better explained by their emergent spatiotemporal dynamics. The analysis demonstrates that local changes in chemical distributions, such as thickening membranes, have consequences for subsequent behavior. In dynamical systems terms, adaptivity is determined by the shape of configurations' basins of attraction, and how interactions with the environment move a system through the phase space. A transition is adaptive if it brings the individual to a configuration with a more accommodating basin.

Ontogenies result from the interplay between an individual's configurations and its environment. Here, we mapped the space of all possible ontogenies for a model protocell in a specific environment by exposing it to sequential perturbations and measuring its responses. The resulting ontogenic network contained rich structure as uncovered by a combination of statistical and graph-theoretic methods, including bottlenecks that constrain ontogenic change and branches that reflect irreversible transitions. The set of configurational morphologies was found to cluster according to morphological similarity, and these clusters roughly corresponded to those identified in the network topology. Exploration of the same protocell exposed to a second environment generated a different ontogenic network with distinct properties. This second network displayed a different branching structure, demonstrated increased irreversibility, and expanded towards a different region of morphology space with unique morphological features.

The viability of individual configurations was quantified as the average expected life span of the protocell. Two notions of adaptivity then followed from this definition of viability: Local adaptivity measured the change in viability resulting from plastic transitions, and global adaptivity measured long-term trends in viability. Surprisingly, the first ontogenic network (ON₁) had an abundance of locally adaptive transitions, and even displayed global adaptivity. Finally, the mechanisms of a minimally adaptive scenario were analyzed, demonstrating how adaptivity can be achieved by distributed processes without requiring an explicit regulatory mechanism.

The combination of this model and analytical techniques provides a foundation for studying the emergence of viability, ontogeny, and adaptivity in more biologically realistic systems. One natural extension along these lines re-conceptualizes viable configurations as members of dynamically richer limit sets. Another recognizes that environments also include sequential perturbations that occur on similar or faster time scales compared to an individual's intrinsic dynamics; these environments would prevent the system from stabilizing at a limit set before the next perturbation arrives and require a new, more continuous conceptualization of ontogenic change. Further, biological individuals and environments are structurally coupled [18], allowing an individual's behavior to induce correlations with its environment which affect future interactions. All of these factors need to be considered in a generalization of adaptivity to real-world biological individuals.

A full ontogenic network is a multigraph (ON = (C, E, R)) with a set of nodes C, a set of edges E, and a set of edge labels identified by the set of perturbations R. This network is generated by applying perturbations (R) with a probability distribution (P) to configurations in C. The search begins with an initial configuration C₁ and the death state C₀. By applying perturbations r ∈ R to C₁, new configurations in C are uncovered. The same set of perturbations are applied recursively to all uncovered configurations until closure is achieved (i.e., every transition results in either a previously characterized configuration or death).

We would like to thank Simon McGregor and Pedro A. M. Mediano for pointing out an error in an earlier draft of this paper. This work was supported by NSF grant IIC-1216739 to R.D.B., by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute, and by the Indiana METACyt Initiative. The software is freely available at https://github.com/eagmon/co-construction_base.