## Abstract

We introduce a spatial model of concentration dynamics that supports the emergence of spatiotemporal inhomogeneities that engage in metabolism–boundary co-construction. These configurations exhibit disintegration following some perturbations, and self-repair in response to others. We define *robustness* as a viable configuration's tendency to return to its prior configuration in response to perturbations, and *plasticity* as a viable configuration's tendency to change to other viable configurations. These properties are demonstrated and quantified in the model, allowing us to map a space of viable configurations and their possible transitions. Combining robustness and plasticity provides a measure of viability as the average expected survival time under ongoing perturbation, and allows us to measure how viability is affected as the configuration undergoes transitions. The framework introduced here is independent of the specific model we used, and is applicable for quantifying robustness, plasticity, and viability in any computational model of artificial life that demonstrates the conditions for viability that we promote.

## 1 Introduction

Autopoiesis is a proposed formal property of living systems, defined as “a network of processes of production (synthesis and destruction) of components, such that these components: 1) continuously regenerate and realize the network of processes that produces them, and 2) constitute the system as a distinguishable unity in the domain in which they exist” [35]. This formulation provides a unique characterization of life that shifts the focus from the material properties of autopoietic systems to their systemic organization. In particular, it emphasizes *viability*—a type of localized spatiotemporal stability. Living systems retain their autopoietic organization despite material and energetic exchange with the environment, but these systems are precarious, and with disintegration their organization is lost (i.e., they die).

The relationship between an autopoietic system and its environment is known as structural coupling; the system undergoes perturbations that have the potential to disrupt its viability, yet a successful autopoietic system temporarily avoids disintegration by engaging with the environment in reciprocal interactions that conserve its identity [22]. Theorists of enactive cognition have argued that this coupling provides grounds for defining a living system as a locus of both behavior and value. Behavior arises through the system's coupled interactions with its environment, and value arises through the necessary maintenance of its own viability [35]. Within this view, autopoietic systems that are capable of reconfiguring to increase their probability of ongoing viability are adaptive, while those that decrease their viability are maladaptive [10].

In order to characterize the viability of autopoietic systems, we propose to decompose the notion into two complementary properties: robustness and plasticity. *Robustness* is the tendency for a system to retain its configuration despite internal or external perturbations. It is used in systems biology to refer to a system's capacity to maintain its function despite perturbation, and there are several research efforts to uncover the mechanisms that bring about biological robustness [20, 2, 33]. Developmental biology uses robustness to refer to a morphology's resistance to variation and its ability to arrive at the same outcome despite varying environmental conditions [4]. *Plasticity* refers to a system's malleability in response to perturbation. For example, neural plasticity is the ability of synapses to reconfigure in response to either their own intrinsic activity or extrinsic signals. The use of the term in developmental biology refers to morphological or behavioral change in response to genetic or environmental influences [36, 4]. Within both systems biology and developmental biology, these terms have not previously been interpreted with regard to the autopoietic notion of viability. In this light, robustness can be defined as a viable configuration's tendency to return to its prior configuration following perturbation, and plasticity can be defined as a viable configuration's tendency to change to other viable configurations in response to perturbation. Together, robustness and plasticity characterize an autopoietic system's capacity to withstand perturbations and remain viable.

In this article, we make use of the decomposition of viability into robustness and plasticity to explore the space of viability for a minimal computational model of an autopoietic system. Several such models have been developed over the last 40 years [23], and have high-lighted different aspects of autopoiesis, such as genesis [30, 34], self-replication [37, 25, 18], evolvability [18], and behavior [29, 13, 16]. Here, we focus on the conditions of viability that arise in stationary autopoietic systems, setting aside the influence of motility for future study. Following many computational implementations of autopoiesis, our model focuses on a spatial system consisting of a metabolism for chemical construction and a boundary that counteracts chemical diffusion and enables the metabolism. While the requirement of a physical boundary has been a source of controversy in the autopoietic literature with respect to its role in identity [31], there is no controversy that for molecular autopoiesis, the tendency for diffusion must by countered by a self-constructed boundary that keeps the system together as an integrated whole. For it to be autopoietic, this boundary must be constructed and maintained by the system.

In the following background section, we divide current models into four broad classes according to their treatments of space and boundary. In Section 3, we propose a lattice-based model of concentration dynamics that supports the emergence of localized spatial inhomogeneities that engage in metabolism–boundary co-construction. These instances of minimal autopoiesis provide a test bed for developing methods to explore the space of viable configurations. Section 4 demonstrates several properties of autopoiesis that arise in this model, including the co-construction of metabolism and boundary, the precariousness of these configurations, and their capacity for self-repair. Finally, in Section 5, we transition to a systematic study of a viable configuration's robustness and plasticity by characterizing its response to perturbations. This analysis allows us to map the space of viable configurations in response to ongoing perturbations, characterize the configuration's basin boundary, and quantify robustness, plasticity, and viability in this domain.

## 2 Computational Autopoiesis

Computational models of autopoiesis differ in their treatment of spatial containment and its relation to metabolism. In this section, we classify previous models according to four broad categories that highlight different aspects of this relationship. The first type accounts for reaction kinetics but not spatiality, the second accounts for basic spatiality but not boundary construction, the third has an explicitly defined boundary, and the fourth consists of spatial models with distinct constructed boundaries.

We begin with arguably the most minimal type of model, which emphasizes reaction kinetics and ignores the spatial dimensions found in molecular systems. This is seen in the hypercycle [14], autocatalysis [19], chemical organization theory [12], and lambda calculus models [15]. This class of models typically assumes well-stirred conditions, in which concentrations are homogeneous across space. Reaction-based models are often more focused on the catalytic closure of a given chemical organization, which emphasizes how molecular sets can synthesize themselves through the reactions of their molecular components. A type of precariousness and robustness has been suggested by Barandiaran and Egbert [3] in a non-spatial model composed of ordinary differential equations. For this class of models, death is simply a stable attractor at which all relevant molecular concentrations are zero, and viable states are those that correspond to a stable equilibrium at positive values of the relevant molecular concentrations.

The second type of autopoietic model incorporates space but does not demonstrate the construction and maintenance of a boundary. For example, Virgo [32] argues that spots seen in the Gray–Scott reaction–diffusion system have the property of individuation. These spots are stable, spatially localized patterns. However, they have no distinct spatial boundary, making it difficult to argue that the spots define themselves as distinct unities in space that are self-regulating with regards to a medium. Similarly, gliders in the game of Life are members of this second class because they construct themselves as spatial entities without distinctly defined boundary components [5].

The third type of model includes spatial boundaries that have been built in a priori [7]. As an example, the model presented by Egbert and Di Paolo [13] has a membrane that is predefined as a circle of connected springs. The sections of this boundary can grow or shrink based on the local number of membrane molecules. This demonstrates boundary–metabolism co-regulation: The boundary contains the required metabolic molecules and keeps them from spilling into the environment, and the boundary requires metabolism because metabolism constructs the membrane molecules. However, the boundary itself does not emerge from the lower-level physics.

The fourth type consists of spatial models that demonstrate constructed membranes. This includes Varela, Maturana, and Uribe's original particle-based model of autopoiesis [30], which is defined on a two-dimensional lattice, each cell of which contains at most one particle at a time. Three types of molecules are included: substrates, catalysts, and links (membrane molecules). Through these molecules' reactions, a spatial boundary arises and contains the internal metabolism that continues to construct it. Similar models have demonstrated growth, change of shape, oscillation, self-reproduction, the formation of flexible membrane, and movement [37, 29, 8]. This class has been extended by modeling particle interactions in a continuous space, and by taking a macroscopic perspective on concentration dynamics rather than modeling the movement of individual molecules [24]. The model presented in this article falls into this fourth class.

There are many lessons we can learn from these past efforts to model autopoiesis. For a model to fully address the precariousness of an autopoietic organization, this organization must emerge from lower-level processes rather than being hardcoded into the model. We suggest that most models of autopoiesis to date have various shortcomings with regard to this requirement. The first type of model, which focuses on stoichiometric relations but ignores spatiality, implicitly hardcodes the distinction between system and environment. The second type of model lacks distinct boundary components that would contain the constituents; therefore such models cannot fully account for metabolism–boundary co-construction. The third type of model hardcodes a boundary, and so it does not account for the boundary's emergence. The fourth type of model, if implemented carefully, is the only class that has the potential to demonstrate both catalytic closure and maintenance of an emergent boundary.

Most models of autopoiesis to date were developed purely as proofs of concept for the autopoietic logic and its consequences, and were not analyzed further. Some exceptions include [3, 7, 5]. In contrast, we strongly believe that once constructed, models of autopoiesis should be subjected to a thorough analysis. It is only through such systematic study that the rich consequences of autopoiesis can be fully articulated.

## 3 Anisotropic Spatial Model

We consider a minimal autocatalytic model consisting of three parts: (1) a set of molecules and the corresponding set of constructive and destructive reactions, (2) spatial diffusion dynamics biased by repulsion, and (3) an anisotropic repulsion potential. These parts are combined in a two-dimensional lattice model derived from a previous model investigated by Ono and Ikegami [24,–26]. We simplified the reaction kinetics, used a rectangular lattice, and implemented continuous orientation fields. In addition, our model is completely deterministic.

*x*

_{ij}, the equations of motion take the form of a system of coupled differential equations:where

*t*is time,

*R*_{m}(

*x*

_{ij},

*t*) indicates rate contributions due to the mass-action assumptions of our chemical reactions, indicates a diffusive term biased by repulsion and orientation field θ, and

**Γ**

_{m}(

*x*

_{ij},

*t*) denotes potential external contributions to the system. Each of these terms will be addressed below.

### 3.1 Autocatalytic Stoichiometry

*A*, the membrane

*M*, the food

*F*, and water

*W*. Each of these molecular types reacts as follows. Autocatalyst

*A*catalyzes the production of additional

*A*and consumes food

*F*in the process, with a rate

*k*

_{AA}. Autocatalyst

*A*also catalyzes membrane

*M*from food

*F*at a rate

*k*

_{AM}. Both

*A*and

*M*decay at rates

*k*

_{A}and

*k*

_{M}, respectively. Water

*W*is an inert substance and is conserved throughout the simulation. These reactions are summarized by the following stoichiometric equations:with all reverse reaction rates set to 0. Conservation of mass can be assumed via the introduction of an inert waste product that is instantaneously removed from the system. Finally, instead of modeling each individual molecule, we coarse-grain over space and consider concentrations of molecules. The mass-action assumptions allow us to convert these stoichiometric relations into dynamical equations for the concentrations:Unlike the aforementioned Ono–Ikegami model [24], the reaction rates are all constant. There are also no constraints on the sum of all concentrations at a specific point. However, because these concentrations represent physical quantities, the model is only defined for positive values of

*A*,

*M*,

*F*,

*W*.

### 3.2 Diffusion and Repulsion

*N*

_{x}×

*N*

_{y}lattice. Each molecular concentration follows diffusive dynamics resulting in conserved concentration flows between the lattice sites:where

*d*

_{m}is the diffusion constant for the

*m*th-type molecular concentration, is the Moore neighborhood of each lattice site

*x*

_{ij}, and is the diffusive bias of concentration

*m*from site

*x*

_{ij}to site

*x*

_{αβ}. This direction-dependent diffusive bias reflects the tendency for molecular concentrations to flow away from repulsive forces and seek minimum values of repulsive interactions. Therefore, the diffusive bias is given as a function of the gradient of a repulsion potential between the sites:where is the gradient of the repulsion potential in the direction of

*x*

_{αβ}. Finally, this repulsion potential for molecular concentration type

*m*at lattice site

*x*

_{ij}takes the formwhich has contributions from all μth-type concentrations in the inclusive Moore neighborhood of

*x*

_{ij}. Within this framework, the specific interactions between concentration types

*m*and μ are determined by a sum of their interaction strength, , weighted by the molecular concentration at each site. This interaction strength is a function of the relative location within the neighborhood (including self-interactions), and the orientation field θ, introduced in the next section for those interactions involving

*M*,

*A*, and

*W*.

### 3.3 Anisotropic Membrane Potentials

Spatial containment relies on the presence of a boundary between elements. By its very definition, such a boundary must behave with a reduced dimension from that of the embedding space. For example, the boundary of a three-dimensional sphere is a two-dimensional surface. Motivated by the dynamics of lipid sheets that form in water, we achieve an approximation to boundary behavior by endowing the membrane *M* with asymmetric repulsion interactions with both autocatalyst *A* and water *W*. As shown in Figure 1a, traditional isotropic interactions cause the membrane concentrations to behave like oil in water, forming large symmetric clumps without any containment properties. The introduction of anisotropic interactions, whose schema is depicted in Figure 1b, results in the membrane forming sheets or layers, as seen in Figure 1c.

These anisotropic interactions are formulated via the introduction of an orientation field that describes the local orientation θ(*x*_{ij}, *t*) ∈ (0, π) of the membrane concentration at each lattice site *x*_{ij}. This orientation variable specifies the direction of the major axis for an elliptical repulsion potential between the membrane concentration *M* and both the *A* and *W* concentrations as shown in Figure 1b. The interaction strength in each of the eight neighboring cells is calculated as the normalized area *C*^{θ}(*x*_{αβ} − *x*_{ij}) of the directionally appropriate π/4 octant of the ellipse, with the −π/8 to π/8 octant directed towards the central right neighbor *x*_{i+1,j}. The lengths of the major and minor axes of this ellipse are given by and , respectively. Because membrane is immiscible with autocatalyst and water, there is a higher repulsion strength between concentrations on the same lattice site than between concentrations on adjacent lattice sites; thus when *x*_{αβ} ≠ *x*_{ij}.

*x*

_{ij},

*t*) at lattice site

*x*

_{ij}aligns with their neighbors with an angular velocity proportional to the concentration of membrane

*M*(

*x*

_{αβ},

*t*) in the neighboring cell

*x*

_{αβ}. Thus, the local orientation has the following dynamics:where ψ(

*x*) = sin(2

*x*) was chosen so as to smoothly stabilize aligned orientations and destabilize anti-aligned orientations.

### 3.4 Full Model

The full model consists of four differential equations in the form of Equation 1 (one for each molecular concentration type ) and one differential equation for the orientation in the form of Equation 7 for each site *x*_{ij} ∈ *N*_{x} × *N*_{y} in the lattice. Thus, the model is technically a 5*N*_{x}*N*_{y}*-*dimensional coupled, nonlinear dynamical system. Since many of the phenomena in which we are interested occur in nonequilibrium, the Γ_{m}(*x*_{ij}, *t*) term allows for general flows of concentrations in or out of the system. This will be particularly useful for modifying the food concentration *F*, which plays the role of an energy supply for the chemical reactions but does not participate in repulsive interactions. Additionally, we note that while the water concentration *W* is reactively inert, it contributes to repulsive interactions and fills what would otherwise be empty space. The spatial distributions of the autocatalyst *A* and membrane *M* concentrations are of primary interest to the following analysis.

The system we consider unfolds on a *N*_{x} = *N*_{y} = 40 square lattice with toroidal boundary conditions, giving a total of 8,000 ordinary differential equations. To keep a steady supply of food entering the system, the concentration of *F* is driven towards a saturation value of *s*_{F} at a rate *k*_{F} via the additional linear term Γ_{F}(*x*_{ij}, *t*) = *k*_{F}(*s*_{F} − *F*(*x*_{ij}, *t*)). Food was only supplied outside a circle of radius 14 cell units, and thus only by diffusion can it reach the region in which configurations appear. No additional interactions with an external source are assumed; thus Γ_{A}(*x*_{ij}, *t*) = Γ_{M}(*x*_{ij}, *t*) = Γ_{W}(*x*_{ij}, *t*) = 0. The general model supports a plethora of possible behaviors for the spatial distribution of the molecular concentrations. Although a systematic study of these parameters was not performed, our parameter selection was informed by two criteria: Repulsion should allow membranes to form in stable sheets, and the tradeoff between diffusion and production rates should support the requirement of containment for ongoing catalysis. The specific values for all parameters used for our simulations are given in the appendix.

As is typically the case in reaction–diffusion models, a previous incarnation of this model [1] used the forward Euler method to advance the states of the system over time. However, due to the numerical sensitivities associated with this approach, in this article we have opted instead to utilize a fourth-order adaptive Runge–Kutta method, as implemented by MATLAB's ode45 ordinary differential equation solver [27].

## 4 Model Behavior

### 4.1 Initiation and Time Evolution

Given any initial spatial configuration for the molecular concentrations and orientation field, the system unfolds deterministically, typically towards a uniform state in which *A* and *M* concentrations have been depleted. To facilitate finding a viable stable configuration, with nonzero values of *A* and *M*, the system was initialized in configuration *IC* (Figure 2). This configuration consists of a small circle of autocatalyst with a diameter of 5 lattice cells at a uniform concentration of 0.9. The circle of autocatalyst is surrounded by a ring of membrane that is 4 lattice cells thick with a concentration of 0.9. The membrane orientation field is initiated as concentric circles centered on the configuration. The center point of the cell has an orientation of 0, which breaks the system's initial rotational symmetry. Food is initiated as a uniform field with a concentration of 0.18; water is initiated with a concentration of 0.1 where there is membrane or autocatalyst, and 1.0 everywhere else. This cell-like initial configuration, with autocatalyst contained by a membrane boundary, facilitates the subsequent evolution towards a viable configuration. No systematic study of initial configurations was performed.

We are primarily interested in the stable configuration *SC* (Figure 2), which arises from the initial configuration after a transient of approximately 16,000 time steps, and then persists for over 100,000 simulated time steps. Because of its persistence over this long time scale, which is much longer than the time scale of its underlying cell processes, we assume that it has reached a stable attractor. The exact configuration of the system at *t* = 25,000 will be resimulated and tested for the remainder of this article. In *SC*, the membrane is thinner than the initial configuration, and there is a distinct bridge formation of membrane cutting across the lower half of the cell and splitting the autocatalyst into two regions. Due to the directional repulsion of the anisotropic membrane molecules, cells with low molecular concentration appear directly to the left and the right of the bridge. This broken symmetry arises in part from the initial broken symmetry of the center site's orientation.

### 4.2 Containment and Metabolism

Of particular importance to models of autopoiesis is the establishment of viable metabolism–boundary co-construction. While the boundary's dependence on metabolism is explicitly made through the reactions in our model, the reverse dependence needs a bit more justification. Here we demonstrate the dependence of the autocatalyst on its containment within the membrane by separating them (Figure 3). The membrane field, water field, and orientation field are each displaced upwards by 16 lattice cells from their original position, while the autocatalyst field and food field are kept in place. This division was selected because the membrane requires orientation and water to retain its shape, and the autocatalyst requires food to continue catalyzing its products.

Whereas *SC* would have lasted for over 100,000 time steps, the displaced configuration is clearly in the process of disintegrating after 200 time steps have passed, and the final vestiges of the relevant molecules disintegrate around *t* = 1,000. These results illustrate the strong interdependence of boundary and metabolism in the process of co-construction. The processes of containment and molecular catalysis counter the effects of diffusion and decay. Without the membrane, the autocatalyst diffuses away. Even though membrane continues to be synthesized, this does not occur fast enough to restore containment. Conversely, the boundary needs the autocatalyst for its synthesis. Soon after separation, the internal bridge elongates and the membrane thins. Within several hundred time steps, only vestiges of the boundary survive, and soon after they are gone. Only through coupled co-constructive processes can the autocatalyst and membrane molecular species persist.

### 4.3 Self-Repair and Disintegration

Figure 4 illustrates *SC*'s precariousness and capacity for self-repair. Two sequences are shown in which a tear perturbation was applied to its membrane, with a larger tear in the bottom sequence. A tear is a perturbation in which a small section of membrane is removed by setting all the concentrations to a value of 0 and randomizing the orientation of those lattice cells. At *t* = 1, we can see the missing membrane section in both sequences.

In the top sequence, which has a smaller tear, the membrane flows back into the removed section and the system fully restabilizes by *t* = 500. The orientation, which is not illustrated in the figure, also realigns and allows the membrane concentration to flow around the boundary. Interestingly, this event results in a different stable configuration than before the perturbation. The cell is elongated on the right side, where it suffered the perturbation, resembling residual outgrowth.

In the bottom sequence, with the larger tear, the autocatalyst begins to leak out, as shown at *t* = 10. Membrane also flows into the tear, but the autocatalyst repels it, and therefore the hole cannot be repaired. By *t* = 500, most of the autocatalyst has leaked out, and the membrane is partially disintegrated. The leaking autocatalyst creates new membrane molecules, which self-assemble additional strings. However, with the containment condition broken, the metabolism cannot persist and the boundary decays. Most of the concentrations of *M* and *A* have been depleted by *t* = 1,000, and soon thereafter the system returns to its terminal, uniform state.

## 5 Robustness and Plasticity

In Section 4, we observed conditions that allow the membrane and the autocatalyst to sustain one another through co-construction and repair themselves following tear perturbations. We now quantify the robustness and plasticity of viable configurations. Here, robustness is simply defined as the capacity to retain configurational stability despite perturbations to the system. Plasticity is defined as the capacity to restabilize at different viable configurations following perturbations. We study these properties by systematically applying perturbations and observing which conditions drive the system to disintegration and which allow it to remain viable.

Quantifying robustness and plasticity immediately highlights three challenges: (1) how to quantify viable stability, (2) how to characterize a system perturbation, and (3) how to characterize a system's response to a perturbation.

Viable configurations are associated with stable, heterogeneous configurations of autocatalyst and membrane. On the other hand, the death state is a stable attractor for which the relevant concentrations are homogeneous. We identify stability with extended periods of time in which all of *A* and *M*'s time derivatives are near zero. Thus, a structural configuration is stable when the average temporal derivatives of *A* and *M* over the whole system satisfy for all *t*_{k} ∈ [*t*, *t* + 5000] and a sufficiently small *ϵ*_{1} = 0.05. A finite-difference approximation was used in calculating this quantity.

There are many types of system perturbations that could be explored. We classify the space of possible perturbations by four main dimensions: *type*, *location*, *time*, and *stochasticity*. *Type* refers to the five different state variables: the molecular types (*A*, *M*, *F*, and *W*) and the orientation field of the membrane (θ). Perturbations can be applied to any combination of these variables. *Location* refers to the positions within the lattice at which these molecules are perturbed. Perturbations can be applied to spatially localized regions, or globally to all points throughout the system. *Time* refers to the temporal structure of perturbations throughout a simulation. *Stochasticity* refers to whether a specific perturbation is applied, or whether it is randomly selected from a class of perturbations.

Ideally, each class of perturbations will be characterized by a small number of parameters, allowing one to explore the system's response as these parameters are varied. For example, the tear perturbation, introduced in Section 4.3, is an instantaneous spatially localized perturbation, applied to *M* and θ and parameterized by the spatial extent of the tear. Another example is global noise, which instantaneously perturbs every variable in the system by a direction and amplitude. The space of possible perturbations is vast.

Finally, we must characterize a system's response to a perturbation. Here, we measure how long it takes for a system to reach stability following a perturbation. Once it stabilizes, we determine its structural distance from the original configuration by a structural comparison function. Specifically, the distance between two stable configurations *C* and *C*′ is found by taking the summed absolute difference between the two configurations' states, . If this sum is near zero (viz., *d*(*C*, *C*′) < *ϵ*_{2} with *ϵ*_{2} = 1.0), then we consider the system to have restabilized at the original configuration and demonstrated robustness. In contrast, larger values indicate a new configuration has been reached. If this new configuration is viable, then the system has demonstrated plasticity. If it is in the death state, at which molecular concentrations are spatially uniform, then the system has disintegrated.

### 5.1 Gaussian Blur Perturbations

The Gaussian blur is an ideal perturbation for systematic study because of its low dimensionality and strong influence over the system's state. It is a global perturbation applied instantaneously to the four molecular species (but not the orientation field θ). Gaussian blur is a common technique from image processing used to smooth an image, reducing noise and detail. For our model, applying the blurring function is akin to instantaneously shaking the system. More shaking leads to more blurring, as molecules between farther regions in space become mixed. Mathematically, the Gaussian blur is a convolution of the concentration fields with a Gaussian function, resulting in a Gaussian-weighted spatial average. By increasing the variance (σ^{2}) of the Gaussian function, we can make the blurring increasingly severe.

Figure 5 shows the time to restabilization following a Gaussian blur perturbation applied to *SC* as σ^{2} was systematically varied. The figure shows that the restabilization time generally increases as the blur becomes more severe. Interestingly, there are several intervals of σ^{2} that destabilize the system so that it can no longer recover and instead evolves to the stable death attractor. This occurs for a short interval at 0.25 < σ^{2} < 0.28, for a longer interval at 2.20 < σ^{2} < 2.68, and for all σ^{2} > 3.36.

Figure 5 also shows the potential of the system to reconfigure following a perturbation, demonstrating its plasticity. In this case, seven new stable viable configurations were reached following the Gaussian blur perturbation. For small σ^{2}, the system returned to *SC*. At σ^{2} > 0.28, we see a discrete change occur; the system takes longer to stabilize, and when it stabilizes, it is at a new configuration *C*_{1}. This configuration looks similar to *SC*, but has a larger internal space for the autocatalyst, and the membrane bulges outwards in some places. As σ^{2} is increased beyond 1.82, configuration *C*_{2} is reached, and later we see configurations *C*_{3} through *C*_{7}. In *C*_{2}, *C*_{3}, *C*_{4}, and *C*_{7}, a broken symmetry exists in the internal membrane stripes. Presumably there are other stable configurations identical to these, but with reversed broken symmetry. For configurations *C*_{4} through *C*_{7}, different internal structures formed in addition to the original internal bridge.

### 5.2 Spatially Localized Gaussian Perturbations

*A*,

*M*,

*F*, and

*W*) in the region around a given focal point. The function is of the formwhere |

*x*

_{ij}−

*x*

_{f}| is the distance from the focal point

*x*

_{f}to the given cell

*x*

_{ij}, α is the amplitude of the function, and σ

^{2}is the variance. For a specific focal point, a subclass of perturbations is defined by varying α and σ

^{2}. An α of 0 produces no change to the concentrations. Negative α reduces the concentrations, and positive α increases the concentrations of the affected cells surrounding the focal point. If a perturbation lowers any concentration value below 0, the value is set to 0.

Figure 6 depicts four examples of spatially localized Gaussian perturbations applied to *SC*, which differ according to their focal points (the bright green dots in the inset configuration diagrams). In these plots, *SC*'s restabilization is shown in green, stabilization to a different configuration is shown in purple, and perturbations that lead to death are shown in red. The restabilization time is indicated by the shade, where darker shading indicates longer restabilization times. Since α = 0 perturbations have no effect on *SC*, the corresponding light green line indicates zero time steps to regain stability. The shading also demonstrates that it takes increasingly longer to restabilize as a boundary between different configurations is approached.

A qualitative comparison of the four perturbation figures reveals that the location of perturbation significantly alters the system's response. The size, shape, and location of the red, green, and purple areas change considerably between the profiles, indicating that the system's tendency towards disintegration, robust destabilization, or plastic change is highly dependent on the focal point for these localized perturbations. Moreover, we can see that a (σ^{2}, α) perturbation applied in one location might lead to death, while the system can recover from the same perturbation applied to a different location. For example, the lower regions of Figure 6a and b are mostly red, while the same regions in Figure 6c and d are green. Figure 6c also has an island of plastic change in the top right corner. This region denotes severe perturbations that, if applied to any of the other three locations, lead to death. These perturbation profiles reveal that a system can respond very differently depending on the class or subclass of perturbation.

### 5.3 Mapping the Network of Viable Configurations

Plasticity implies that the system can transition to different configurations following a perturbation. In the previous perturbation studies, we saw evidence that such transitions were prevalent in our model. Specifically, when they were systematically varied, parameterized perturbations, applied once to a single configuration *SC*, resulted in the system restabilizing in new viable configurations. These new configurations have different spatial dependences, and are likely to have distinct perturbation response profiles. Likewise, each of these can be explored and evaluated in terms of robustness and plasticity. Here, we begin an exploration of the structure of the space of viable configurations by mapping the possible transitions between configurations due to a specific class of sequential perturbations: the Gaussian blur perturbation (introduced in Section 5.1).

Starting from *SC*, successive Gaussian blur perturbations, with 0.0 ≤ σ^{2} ≤ 4.0 and sampled with a resolution of 0.01, were applied to the system. The perturbations were timed so that the system demonstrated it had restabilized (according to the criteria specified in the introduction to Section 5) before the next perturbation was applied. Configurations *C*_{1} to *C*_{7}, which were seen in Section 5.1, reappear here as the configurations directly reachable from *SC*. These configurations were then tested with the same set of Gaussian perturbations, and the investigation terminated when all combinations of perturbations led to death or returned the system to a previously characterized configuration.

The resulting space of configurations can be represented by the network shown in Figure 7a. Within this network, nodes correspond to viable configurations, and directed edges indicate possible transitions following a perturbation. The proportion of perturbations that cause each transition are indicated by that transition's edge weight (line width). This network representation is an extended application of the methodology used to map transitions in response to perturbations between a glider's configurations in the Game of Life cellular automata [6].

Expressed as a network, the intricate structure of the space of viable configurations with respect to the Gaussian blur perturbation becomes apparent. Self-loops (green edges connecting a configuration to itself) reflect the robustness of a specific configuration. If we focus on a single configuration, the sum of the weights for all purple edges defines the instantaneous plasticity of that configuration. Similarly, the sum of the weights for purple and green edges defines a quantification for the instantaneous viability of a configuration.

If we consider sequences of perturbations, it becomes clear that these instantaneous notions must be extended to reflect the potential number of perturbations before the system enters the death state—a property of all possible paths in the network. To this end, the network is interpreted as the backbone representation for a Markov process describing the transitions between configurations given Gaussian blur perturbations with random σ^{2} (here we assume a uniform distribution for σ^{2}). This Markov process has an absorbing state corresponding to the death state of the model [17]. The mean survival time, or expected number of perturbations before disintegration when starting at each of the possible configurations, is shown in Figure 7b. For this specific example, *SC* has the longest mean survival time and all plastic changes tend to decrease viability. This might not be a surprise, as most living systems that we are familiar with in the real world are also likely to become weaker if their contents were to be blended. There is one exception to this, the transition from configuration *C*_{9} to *C*_{13}. Inspection of the configurations' mean survival times (1.197 and 1.234 respectively) reveals that this transition slightly increases the system's chances of surviving a random Gaussian perturbation. Other classes of perturbations, or other initial configurations, might exhibit a larger tendency to contain transitions that increase viability.

### 5.4 Exploring a Configuration Basin

Technically, the system is an 8000-dimensional dynamical system, and configuration *SC* is an attractor of the dynamics. As with any attractor, *SC* has a basin of attraction—a set of system states that return to it. Perturbations to *SC* move the system to a new point in the state space. If the system returns to *SC*, then the perturbation has remained within its basin of attraction. Those perturbations that resulted in new configurations or death push the system into the basins of attraction for those other attractors. The phase portrait of the dynamical system captures the decomposition of the full state space into all of the possible basins of attraction, and provides a complete characterization of the model in the presence of any type of perturbation.

This phase portrait is illustrated by a conceptual sketch shown in Figure 8a. The idealized basin of attraction for *SC* (shown in green) shares many properties with the actual basin. A large portion of the basin is bounded by the physical constraint that concentrations must be nonzero (represented by the brown boundary in the state space). The basin is also contiguous with several other basins of attraction, including a large death state (shown in red) and other configurations (shown in purple). While general, unstructured perturbations would move the system to any point in the state space, structured perturbations could only access a subspace of the possible states (illustrated by the light yellow region). Given a phase portrait and unstructured perturbations, the proportion of perturbations that result in a configuration corresponds to the relative volume of a configuration's basin. For structured perturbations, the proportion of perturbations that result in a configuration corresponds to the relative volume of the intersection between the basin and the subspace accessible by the perturbation. Likewise, the probability of transitioning between configurations given a randomly chosen direction to perturb the system corresponds to the relative proportion of the shared surface area between any two basins.

Due to the massive size of this system, it would be extremely difficult to exactly define the full boundary for *SC*'s basin of attraction. Here, we resort to a restricted statistical characterization of the boundary. Specifically, we concentrate only on the 177 lattice sites that fall within a circle of radius 7 centered on the membrane's focal point. This limits our search to an 885-dimensional subspace of the full model's state space. We also preemptively allow for the fact that most randomly selected directions will cross the physical boundary associated with non-negative molecular concentrations, by restricting our search to only explore directions in which concentrations are added to the 467 dimensions with near-0 concentration values (less than 0.1).

Such a search starts by randomly choosing an 885-dimensional unit vector according to Marsaglia's algorithm [21]. This vector is scaled by a distance *d* and used to additively perturb the states of *SC*. If the system returns to *SC*, then the perturbation remained within its basin of attraction. We searched along each of these randomly given directions to find the maximum value of *d* such that we remained within the basin of attraction for *SC* (with a minimal resolution of 0.01). We also took one step just beyond the boundary and recorded the new configuration in which the system stabilized.

The statistical distribution for 2000 samples of the maximum values of *d* that remain within the basin of attraction for *SC* is shown in Figure 8b. The mean of the distribution is 0.8993, and its variance is 0.0781. It has several interesting properties. The fact that the support of the distribution starts at *d* ≈ 0.4, and not *d* ≈ 0.0, is consistent with describing *SC* as relatively robust. This reflects a nontrivial minimum distance between *SC* and its basin boundary (nearly 20% of the maximum found distance of 2.0). The multimodal nature of the distribution may reflect the fact that the boundary is not convex and instead undulates with the 884-dimensional equivalent of scattered crests and valleys. This exploration also reveals that there are only three other configurations that share significant proportions of *SC*'s basin boundary (46%, 26%, and 17% of the samples transitioned to each configuration, respectively). The immediately contiguous death basin boundary is relatively much smaller (only 2% of samples resulted in death). A few samples also found the physical concentration boundary (9%). Surprisingly, plastic transitions are far more likely than death.

Interestingly, none of the configurations found by the Gaussian blur or localized Gaussian perturbations were observed here. There are several possible reasons for this. One possibility is that the random directions were undersampled. Another possibility is that the basin boundary was not fully explored due to shadowing in nonconvex volumes. For example, in Figure 8a, the area indicated by the striped green pattern cannot be revealed by rays extending from the configuration's attractor. A third possibility is that *SC*'s basin of attraction does not directly share a boundary with those configurations. Instead, they were uncovered by large perturbations that jumped over the basins of *SC*'s neighboring configurations.

The robustness of *SC* to unstructured perturbations is associated with the volume of its basin of attraction. In order to approximate this volume, we used a Monte Carlo method to sample 10,000 points within a restricted 885-dimensional hypersphere of radius 2.0 surrounding *SC*. Within this restricted hypersphere, 47% of the samples returned to *SC*. From this, an approximation of the volume can be calculated. Some of *SC*'s basin likely falls outside of the sampled region, and is therefore not reflected in this approximation.

Somewhat counterintuitively, the probability of transitioning between one configuration and another can be quite different than that of the reverse transition. This phenomenon has been studied in the context of evolution by Conrad [9] and Stadler et al. [28]. Here, we see that the phenomenon occurs because the contiguous surface area of the two configurations' basins corresponds to different proportions of each basin's total surface area. Thus, it is more likely that transitions occur from the smaller basin to the larger one than the other way around.

## 6 Discussion

This article has introduced a spatial model of concentration dynamics, which, through the processes of diffusion, repulsion, and chemical reactions, supported the emergence of stable spatiotemporal inhomogeneities that we called viable configurations. These configurations were shown to demonstrate metabolism–boundary co-construction, precariousness, and self-repair. Tools from dynamical systems theory were used to map the configurations' response to structured classes of perturbations. A network summary of the space of viable configurations featured three important types of responses to perturbation: robust (recovery of the same configuration), plastic (transition to a different viable configuration), and disintegration. This analysis illustrated the intricate structure underlying the space of viable configurations and provided a concrete basis for investigating autopoiesis and its consequences.

The spatial model discussed here was derived as an abstraction of idealized molecular dynamics. In the course of its derivation, several approximations resulted in its formulation as a continuous-state lattice model. These approximations could be further relaxed to derive a model fully continuous in state and space; the resulting system would be an instance of a nonlinear, partial integrodifferential equation. An important feature of general spatiotemporal systems is that the relevant configurations emerge as sustained, localized spatial inhomogeneities instead of prespecified system–environment distinctions. These configurations are precarious in the sense that their spatial co-dependences can dissolve, leaving the system in a homogeneous terminal state.

Dynamical systems theory provided a convenient framework for the analysis of our spatiotemporal lattice model's response to perturbations. The model's system of dynamical equations defined a phase portrait. The specific configurations studied here were identified as equilibrium points of the concentration dynamics. It is important to note that these points are chemically still active: Concentrations are continuously being depleted and regenerated through the balance of reaction and diffusion dynamics. Perturbations moved the system around the phase space and allowed it to transition between different basins of attraction. In particular, the system's response to perturbation can be related to the relative volume and surface area of attractor basins. Section 5.4 provided a preliminary investigation of a single configuration's basin. The system's movement through the phase space was also shown to be highly dependent on the class of perturbations considered. Accordingly, conclusions about the perturbed phase space should be conditioned on the class of perturbations. Many of these conclusions can change for different perturbation classes.

Grounded in this dynamical systems view, robustness, plasticity, and disintegration were characterized as different properties of a viable configuration's response to perturbation. The structure of possible transitions induced by a class of perturbations can be summarized as a transition network, as shown in Section 5.3. Such networks can have many interesting properties that capture attributes of the system's viability. For example, configurations without transitions to the death state would indicate strongly viable configurations, while configurations whose only outgoing edge is a self-loop would be completely robust or invariant to the perturbation. Cycles could also appear in the network structure, capturing the possibility of a system returning to a previous configuration following a series of perturbations. Asymmetric bidirectional transitions would indicate that the underlying dynamics can perform a forward transition more easily than the reverse transition. Finally, the network might have a large branching structure, allowing it to make long chains of contingent distinctions about its environment.

The model and analytical framework presented here can be used to ground several of the ongoing discussions on value, adaptivity, and behavior in enactive cognitive science. This field of research emphasizes the dynamic interaction between an autopoietic system and its environment. Here, a situation's value to the individual is defined by the extent to which it affects the individual's viability [11]. By interpreting a situation as a perturbation, this now becomes measurable; a transition induced by a given perturbation can increase or decrease a configuration's viability, and therefore the perturbation can be interpreted as good or bad for the configuration. Adaptivity is also focused on an extension of viability that captures a system's tendency to be altered away from disintegration, or towards increasingly viable states. By associating adaptivity with increases in quantifications of viability, and identifying related attributes of a model and its dynamic structure, we can explore questions related to the minimal functional requirements of adaptivity and classify systems accordingly. Finally, the effect of behavior on an autopoietic system's viability can be explored in this model through the addition of motility. Such an extension would require methods for characterizing the intertwined relations between a viable configuration, its emergent behavioral tendencies, and the location-dependent perturbational classes to which it is exposed.

## Acknowledgments

We would like to thank Eduardo Izquierdo, Matthew Egbert, and the anonymous reviewers for feedback on an earlier version of this article, and Zach Haga for his efforts in collecting data. This work was supported in part by NSF grant IIC-1216739 to R.D.B., as well as NSF IGERT fellowships to E.A. and A.J.G. This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute, and in part by the Indiana METACyt Initiative. The software is freely available at https://github.com/eagmon/co-construction_base.

## References

### Appendix: Model Parameters

See Table 1.

Property . | Symbol . | Value . |
---|---|---|

Diffusion | d_{A} | 0.25 |

d_{M} | 0.25 | |

d_{F} | 0.25 | |

d_{W} | 0.25 | |

Reaction rates | k_{A} | 0.005 |

k_{M} | 0.005 | |

k_{AA} | 0.25 | |

k_{AM} | 0.3 | |

Orientation dynamics | k_{θ} | 0.01 |

Food dynamics | k_{F} | 0.8 |

s_{F} | 0.18 | |

Repulsive interaction—same site | 0.3 | |

0.3 | ||

9.0 | ||

9.0 | ||

Repulsive interaction—neighbor site | 0.0125 | |

0.0125 | ||

3.0 ∗ C^{θ}(x_{αβ} − x_{ij}) | ||

3.0 ∗ C^{θ}(x_{αβ} − x_{ij}) | ||

Repulsive interaction—ellipse | 20 | |

1 |

Property . | Symbol . | Value . |
---|---|---|

Diffusion | d_{A} | 0.25 |

d_{M} | 0.25 | |

d_{F} | 0.25 | |

d_{W} | 0.25 | |

Reaction rates | k_{A} | 0.005 |

k_{M} | 0.005 | |

k_{AA} | 0.25 | |

k_{AM} | 0.3 | |

Orientation dynamics | k_{θ} | 0.01 |

Food dynamics | k_{F} | 0.8 |

s_{F} | 0.18 | |

Repulsive interaction—same site | 0.3 | |

0.3 | ||

9.0 | ||

9.0 | ||

Repulsive interaction—neighbor site | 0.0125 | |

0.0125 | ||

3.0 ∗ C^{θ}(x_{αβ} − x_{ij}) | ||

3.0 ∗ C^{θ}(x_{αβ} − x_{ij}) | ||

Repulsive interaction—ellipse | 20 | |

1 |

## Author notes

Contact author.

Cognitive Science Program, Indiana University, Bloomington, IN 47406. E-mail: agmon.eran@gmail.com (E.A.); ajgates@indiana.edu (A.J.G.); rdbeer@indiana.edu (R.D.B.)

School of Informatics and Computing, Indiana University, Bloomington, IN 47406.

Ecology and Evolution Unit, Okinawa Institute of Science and Technology, Okinawa, Japan. E-mail: valentin.churvay@oist.jp (V.C.)