A model of replicating coupled oscillators generates naturally occurring cell networks Free

A network of cells joined through intercellular bridges (ICBs) is formed when a founder cell and its progeny undergo incomplete divisions (Chaigne and Brunet, 2022). Examples of networks built by this process include colonies of choanoflagellates (Fairclough et al., 2010; Larson et al., 2020) – the closest living relatives of animals – and both male and female germline cysts (Fawcett et al., 1959; Fawcett, 1961; Haglund et al., 2011; Lu et al., 2017). The first example demonstrates the potential role of this process in the evolution of multicellularity, whereas the second example illustrates its importance in animal reproduction. In both cases, the networks are typically ‘small’, having fewer than 100 cells. Despite the ubiquity of small cell networks in nature, a theory for their generation via coupled, transient oscillations is lacking.

We turn to perhaps the best-studied class of cellular networks: the ovarian germline cysts of invertebrates (Fig. 1), which are well suited for mathematical modeling because their final topology is often stereotyped within a given species, but variable across species. These networks form inside compartments of the ovary that function like assembly lines (Fig. S1). At one end of each compartment, a germline stem cell generates a daughter cell called a cystoblast. As the cystoblast moves away from the germline stem cell, the cystoblast and its progeny divide incompletely, resulting in a network of cells (germline cyst) connected by ICBs. A subset of these cells, typically only one, becomes oocytes. The resulting ovarian germline cysts can be divided into four categories based on topology (Fig. 1C-F). The first class of cysts is linear chains of various lengths. The cysts in the second class contain 2ⁿ cells arranged in a topology called ‘maximally branched’ and result from n rounds of synchronous divisions (if no divisions split existing edges; see Fig. S2). Cysts in the third class contain 2ⁿ cells, but the cells form linear chains connected to each other via a branched core. The last class contains branched cysts with a non-2ⁿ number of cells.

With this diversity of natural forms in mind, and building on experimental studies in the fruit fly Drosophila melanogaster (Fig. S1), we devised a model by which such networks are assembled step by step through the mitotic addition of nodes and edges. The model deterministically maps a single-cell initial condition to a final network, the topology of which depends on the value of several key parameters such as the strength of cell-cell coupling and the asymmetry of cell division. We demonstrate that the model is able to generate nearly all reported invertebrate cysts and identify their relative locations in the space of model parameters. This in turn allows us to characterize neighboring regions where a small change in parameters causes one final network to be converted to another. We discuss how our framework could guide genetic manipulations of network formation, while also providing insight into the origin of intra-species network variability. Finally, we discuss generalizations based on alternative rules for incorporating new nodes (Alsous et al., 2021; Diegmiller et al., 2023) (Fig. S2), which will expand the space of possible graphs that can be generated through a finite sequence of incomplete divisions.

We present a model of multicellular network assembly in which each node is a cell cycle oscillator. In brief, each full oscillation of a node (the mother cell) generates another oscillator (the daughter cell) to which it remains diffusively coupled (Kouvaris et al., 2012) via a network edge (ICB). The network reaches some final topology once all cells have stopped dividing. The model is outlined in Fig. 2, and Table S1 provides an overview of terms and expressions. Additional details and modeling assumptions are in Materials and Methods and the supplementary Materials and Methods.

The parameter z controls whether the above system will tend to a unique stable fixed point or a limit cycle by shifting where the nullclines for x and y intersect. Intersection on the left or right segments of g(x) (which defines the x-nullcline) implies a stable fixed point, whereas intersection on the middle segment leads to limit-cycle oscillations (see Fig. 2A and Materials and Methods).

Accordingly, the initiation and termination of oscillations can be controlled by moving the system into and then out of the oscillatory regime. Motivated by studies of D. melanogaster germline cysts (Mathieu and Huynh, 2017; Insco et al., 2009; Ji et al., 2017; McKearin and Spradling, 1990; McKearin and Ohlstein, 1995), we considered a time-dependent regulation h(t) in Eqn 1 with a ‘pulsed’ form. To implement this, we assumed for simplicity that z(t) is a triangular pulse (see Fig. 2A,B and Materials and Methods). Because z(t) positively regulates y, it promotes mitotic entry and then exit by translating the y-nullcline into and out of the oscillatory regime. The profile of z(t) will therefore strongly affect the number of transient oscillations of both x and y and, thus, the finite growth of the network (outlined later).

Oscillations of the cell cycle components x_a and y_a drive the incomplete division of a mother cell a into two cells (Tyson and Novak, 2008; Hinnant et al., 2020, 2017) that remain physically connected by an ICB. In our model, a cellular division event occurs whenever an indicator variable for the oscillations (in our case, x; see Materials and Methods) cyclically crosses a threshold (vertical dashed lines in Fig. 2A). After the criterion for a division event is satisfied, two network attributes must be specified: the expanded adjacency matrix and the post-division cell states.

First, to update the adjacency matrix at each division, we assumed that the newly generated cell connects to its mother without altering existing edges (Fig. 2D); thus, in-line divisions, which insert new progeny between existing cells (as illustrated in Fig. S2), do not occur. This assumption is based on studies in female germline cysts of D. melanogaster and other insects (Diegmiller et al., 2022; McGrail and Hays, 1997; Storto and King, 1989; Telfer, 1975) (see Discussion for alternative choices; Alsous et al., 2021; Diegmiller et al., 2023; Swiatek and Urbisz, 2019). The ‘post-division mother’ occupies the original position of the mother in the network, whereas the ‘post-division daughter’ connects via the new edge.

Second, to relate the state of the pre-division mother a to the states of the post-division mother a′ and daughter b′, we specified the maps and ⁠. When both maps are symmetric, the model only generates ‘maximally branched’ networks of 2ⁿ cells (Fig. 1D). As a minimal implementation of asymmetric division (de Cuevas and Spradling, 1998; McGrail and Hays, 1997; Diegmiller et al., 2023), we assumed copying of the oscillator variables x_a → (x_a, x_a), but allowed the regulatory variable z to be asymmetrically inherited as z_a → [(1−δ)z_a, (1 + δ)z_a] with |δ | < 1, where δ is the division asymmetry (Fig. 2D). A non-zero δ generates z differences between progeny and, thus, introduces heterogeneity in the network of oscillators by modifying their respective periods T(z) (see Materials and Methods and supplementary Materials and Methods for further details).

To illustrate how the model can generate naturally occurring networks from a single cell, Fig. 3 displays the model-predicted sequence of divisions for a parameter set where the final network matches the main variant of female Chrysopa perla (green lacewing, see Fig. 1F) germline cysts (Rousset, 1978). At time t = 0, there is one cell (cell 1) for which the states x and y reside at a unique stable fixed point. As the pulsed variable z increases, this fixed point shifts until becoming destabilized, giving way to limit-cycle oscillations (Fig. 2A). After one cycle (Fig. 2B), cell 1 produces cell 2 (Fig. 3A), while retaining more z than cell 2. Both complete their next oscillation nearly synchronously, generating a linear four-cell network. After this, the z variable declines. Cells 3 and 4 divide nearly synchronously, whereas cells 1 and 2 do not divide (Fig. 3A,B). The result is a linear six-cell network in which each cell undergoes one final division nearly synchronously.

Not only does the generated final network (Fig. 3A) agree with the primary C. perla variant, but the division history also matches the order of events inferred from fixed samples (Rousset, 1978), with cells 1 and 2 skipping the third round of mitosis before dividing in the final round. In our model, this ‘skipping’ phenomenon occurs due to a transiently desynchronizing effect of moderate diffusive coupling (explored below). As the z pulse limits the time window over which oscillations can occur, transient delays have a permanent effect on the structure of the network.

Finally, we note that other regions of parameter space (with much stronger cell-cell coupling) also produce this primary lacewing variant. However, they have qualitatively different network formation trajectories (see Fig. S3) and, importantly, do not neighbor secondary and tertiary lacewing variants in parameter space (explored below; see Materials and Methods and Fig. S4).

We next studied which of the natural networks (Fig. 1D) the model is able to generate by performing multi-dimensional parameter scans along pulse strength v, division asymmetry δ and diffusion c (see Fig. 4 and Materials and Methods). Each point in parameter space generates a particular final network, in the same manner as in Fig. 3. The resolution of the parameter sweeps allowed us to identify contiguous region(s) of parameter space that generate the same final network (Fig. S4). This allowed us to define distinct parameter regions as being ‘neighbors’ when their surfaces are in direct contact (at the resolution of our search). Given two neighboring regions, small variations in v, c or δ could allow one to generate either of the associated final networks.

A key example relating one network (the D. melanogaster cyst) to its neighbors is shown in Fig. 4A. Among the neighboring natural networks, we noted that two of them, the Gyrinus natator (whirligig beetle) and Habrobracon juglandis (parasitic wasp) networks, correspond to established genetic perturbations of the wild-type D. melanogaster network: they can be reached by mutations in either half pint (16 to 8 cells) (Van Buskirk and Schüpbach, 2002) or encore (16 to 32 cells) (Hawkins et al., 1996; Ohlmeyer and Schupbach, 2003). We also found that small variations in parameters could transform the D. melanogaster cyst into linear networks of various lengths, as well as non-2ⁿ branched networks, such as the two most common C. perla variants (see Fig. S5B).

Besides the reported natural networks from Fig. 1, many as-yet unreported networks neighbor the region of parameter space that generates D. melanogaster cysts. As examples, Fig. 4 shows the three with the largest contact boundaries (see Materials and Methods). Inspired by the half pint name, we named two of these -pint and -pint owing to the loss and gain of four cells, respectively, relative to the D. melanogaster network. We termed the remaining network one-short because it differs from the D. melanogaster cyst by one cell. Model-predicted division sequences for generating these three networks are depicted in Fig. S6A-C.

Aside from networks directly neighboring the D. melanogaster network, we identified regions for each of the linear cyst structures from two to 16 cells, one of the reported cysts in Osmylus fulvicephalus (giant lacewing) and the Nasonia vitripennis (jewel wasp) germline cyst (Fig. 4B). The model also produces a collection of linear cysts of larger lengths (not shown in Fig. 4B). Although all these networks appear in the parameter space, some, like that of D. melanogaster or Panorpa communis (scorpionfly), occupied large parameter regions relative to others, such as that of O. fulvicephalus or N. vitripennis. Some natural cyst structures – like that of Dactylobiotus parthenogeneticus (water bear), Linepithema humile (Argentine ant) or Bombus terrestris (buff-tailed bumblebee) – were not yet identified by this initial search. In the supplementary Materials and Methods, we describe an extension of the model that enables the generation of both the D. parthenogeneticus and L. humile networks by allowing the division asymmetry δ to depend on the number of previous divisions of the mother.

To develop an understanding for the relationships between the networks identified above (Fig. 4), we next considered the dependence of the final network on one parameter at a time (Fig. 5). We focused on the pulse strength v (Fig. 5A) and the coupling strength c (Fig. 5C) because of their generically non-monotonic effect on network size and their experimental relevance.

The dependence on v is most easily understood for the case of symmetric divisions (δ=0). Because the cells are identical in this case, the model can only generate networks of size 2ⁿ, where n is the number of rounds of division. Fig. 5B provides a heuristic explanation of the dependence of n on v. These rely on knowing the period of oscillations T(z), which we estimate in the supplementary Materials and Methods. When z is fixed, the number of cycles is the ratio of the time spent in the oscillatory region to the period T(z). For the case of changing z(t) – and thus a changing cycle period T(z) – the analogous computation depends on the integral of with respect to time t (see Materials and Methods). Fig. 5A shows the heuristic result and the simulation data: over a narrow range in v, n rapidly increases from 0 to 5, thereafter decreasing in a series of plateaus. Fig. S5A shows how these plateaus are altered in the case of relatively strong diffusion c and division asymmetry δ.

Cell coupling (c) has particularly interesting effects on final network structure: whereas the network topology itself is modified by the oscillations, the extent of oscillator synchrony depends on network topology. To understand the role of c, one can consider the case of two coupled oscillators – each with static z values – that are initially desynchronized (Fig. 5D). One can show (see supplementary Materials and Methods) that both cells will oscillate for sufficiently weak or strong coupling but, for intermediate values of c, they become trapped in a non-oscillatory collective fixed point. This effect is commonly referred to as ‘oscillator death’ (Bar-Eli, 1985, 1984; Ermentrout, 1990) in networks of fixed topology. Although the pulsed bifurcation parameter prevents our networks from experiencing indefinite oscillator death, the effect of transient desynchronization leaves its mark on the final network structure.

Fig. 5C shows a characteristic signature of the transient desynchronization associated with intermediate coupling. At weak coupling, the final network is that of D. melanogaster because the division asymmetry is weak enough to allow four rounds of mostly synchronous divisions. For intermediate coupling, there is a sharp transition from the D. melanogaster network to a set of small linear cysts with at most one branch. For strong coupling, division synchrony is recovered, and the final network returns to that of D. melanogaster. See Fig. S5B for an analogous plot showing a more gradual effect of diffusion on network structure.

Taken together, these data give a physical understanding of how, by variation of parameters v and c, the model generates a diverse set of naturally occurring topologies in invertebrate ovaries (Fig. 1). In the Discussion, we consider how these principles for network formation – combined with existing and future experimental data – can further constrain where a given species sits in parameter space and thus the model-predicted effects of perturbations.

The generation of cellular networks by incomplete divisions is a gamete-production strategy observed commonly in metazoans, and the resulting network topologies differ considerably across species. Despite the evolutionary significance of these networks, we lack a quantitative framework for network formation and inter-species variability. Motivated by studies of invertebrate ovaries in which such networks vary minimally within a species, we developed a model of network assembly based on replicating oscillators. A regulatory pulse of strength v promotes transient oscillations on network nodes (cells). Each oscillation causes an incomplete division that leaves the resulting cells coupled by diffusion with rate c. Asymmetric division (δ) of the regulator allows heterogeneity between the resulting cells. This model produces most known germline cyst topologies in invertebrates by variation of v, c and δ (Fig. 4). For networks formed by n rounds of synchronous divisions, our model explains how changing v determines n, whereas for networks formed only by asynchronous division patterns, our model reveals how the interplay between δ and c mediates such choreography. This analysis identifies where in parameter space natural networks reside and predicts parameter changes that transform one final network into another.

The network transitions predicted by our model can be tested in D. melanogaster using genetic techniques. Deciding which genes to manipulate requires knowledge of where the corresponding gene products fit into our model. Fortunately, existing experimental evidence relates some specific proteins to our model parameters, particularly for the regulatory pulse z. For example, a protein called Bag-of-marbles (Bam) has a pulsed temporal profile and directly affects Cyclin A levels (Ji et al., 2017), thus tuning the time window of proliferation. Changing Bam expression has corresponding effects on network size, but this has been mostly characterized in males (Insco et al., 2009). Thus, future experimental measurements of the time course of Bam expression in females (wild-type and mutant) are required to establish its effect on model parameters such as v. Although there are no obvious individual genes controlling either cell-cell coupling or division asymmetry, mutations affecting a shared network-spanning organelle (called the fusome, Fig. S1C) cause D. melanogaster to generate non-2ⁿ networks (Snapp et al., 2004; Lin et al., 1994; McGrail and Hays, 1997; de Cuevas et al., 1996; Yue and Spradling, 1992) and plausibly correspond to changes in c and δ. Future work linking specific mutations to shifts in model parameters will enable model-informed experimental realization of a diverse set of topologies – including as-yet unreported networks predicted by the model such as -pint, -pint and one-short – in D. melanogaster ovaries.

A key additional test of the model relates not to the final network topology, but rather to the model-predicted trajectory by which the network forms. These tests require only the identification of network intermediates (e.g. by analyzing fixed germariums) and are particularly useful for species in which networks form by asynchronous divisions (e.g. N. vitripennis, see Fig. S6D) (Eastin et al., 2020). To illustrate how the model-predicted history can be tested, we considered C. perla. Its primary variant has 12 cells and appears in three large regions of comparable size in our model parameter space. For two of these regions, the network forms by three rounds of synchronous divisions, then, only four of eight cells divide in the last round (Fig. S3); in the other region, the two central cells skip the third round of divisions, then join all others in the last round (Fig. 3). The latter region matches the experimentally inferred history (Rousset, 1978), and small parameter variations in that region – but not in the other – generate other reported C. perla variants (Fig. S4). The future experimental reconstructions of such division sequences in other species will help to constrain where each wild-type species resides in parameter space, thus refining the model-predicted response to experimental perturbations.

Our work also establishes a foundation for studying networks that form less stereotypically, such as the germline cysts in invertebrate testes (Diegmiller et al., 2023; Rasmussen, 1973). In our model, sperm cyst variability can be captured by modifying the effect of division on the network: compared to our previous assumption that divisions never modify existing edges, in sperm cysts, a cell can instead divide into an existing network edge (Diegmiller et al., 2023). Whenever a new cell is generated, the potentially stochastic choice between these two modes results in graph variability (Fig. S2). This variability causes, for each point in parameter space, a deterministic final network to become a distribution of networks. Rather than transforming one graph into another, parameter changes shift the probability density in this space of graphs. Such a framework would enable the study of germline cyst formation even in mammals, in which network fragmentation causes additional variability (Lei and Spradling, 2013). From choanoflagellate colonies to the germline cysts in metazoans, the evolutionary importance of cell networks demands a unified mathematical framework for their formation by oscillations, for which this work provides a key first step.

The assumptions underlying the presented model are based on experiments primarily in D. melanogaster. Some assumptions may be violated in some other invertebrates; nonetheless, it is useful to construct a mathematical model based on a well-characterized species and then to test whether such a model is sufficiently rich to capture variation across species.

The two most basic assumptions of the presented model are: (1) no cell death occurs during network formation and (2) the ICBs formed by incomplete division remain intact. Regarding assumption 1, although the generation of an oocyte – from this network of cells – eventually requires the death of many cells in the network, this death does not occur during cyst formation, but rather much later (Foley and Cooley, 1998). In fact, even in response to starvation of the fruit fly, germ cell death occurs only after the full network has formed (Drummond-Barbosa and Spradling, 2001). Regarding assumption 2, in some invertebrates, some cells initially connected by ICBs disconnect (a process called fission or fragmentation) late in the network formation or soon after the network formation (Büning, 1994). These cases are beyond the scope of this study.

In addition to assumptions 1 and 2, the model contains five additional assumptions: (3) which intracellular factors cause cell division, (4) how the state of one cell affects the state of another cell, (5) the incorporation of new cells into the network, (6) how mother-daughter asymmetries at divisions affect cell cycles and (7) how the arrest of divisions is regulated. As outlined below, some of these assumptions – such as 3, 5 and 7 – are more grounded in experimental data than others, such as 4 and 6.

For assumption 3, we assumed that the element that drives network growth is an oscillator since oscillating levels of activated and inactivated cyclin-CDK complexes control when a cell divides (Diegmiller et al., 2022; Alsous et al., 2021; Hinnant et al., 2017; Tyson, 1991; Tyson and Novak, 2008; Novak and Tyson, 1993; Csikász-Nagy et al., 2006). The existing quantitative characterization of cell cycles in the germline indicates that they are complex, including not only S and M phases, but also gap phases (Hinnant et al., 2017). In the absence of measurements of cell cycle regulators from live-imaged cysts to constrain the parameters of such a model, we treated each cell as a simple two-dimensional relaxation oscillator – motivated by cell cycles in early embryonic development of D. melanogaster and Xenopus laevis (Strogatz, 2018; Diegmiller et al., 2022; Yang and Ferrell, 2013; Ferrell et al., 2011; Mattingly et al., 2017; Deneke et al., 2016; Novak and Tyson, 1993).

For assumption 4, we assumed that the oscillations of one cell can affect those of another cell by diffusion over the network edges (ICBs) (Alsous et al., 2021; Doherty et al., 2021; Diegmiller et al., 2022). Cell cycle coupling has long been hypothesized to be mediated by the fusome (Fig. S1C). Localization of a cyclin (in this case, active Cyclin A-CDK1) to the fusome could promote a trigger wave (Gelens et al., 2014) along the fusome, if proteins promoting Cyclin A-CDK1 activation are present there (Lilly et al., 2000). Importantly, this coupling seems to require something like this reaction-diffusion scheme along the fusome (Lilly et al., 2000), as opposed to diffusion through shared cytoplasm alone (Snapp et al., 2004).

For assumption 5, we assumed that, as observed in the female germline clusters of many invertebrates, each division introduces a node-edge pair, without altering existing edges. Motivating this assumption is experimental work that showed that the fusome (Fig. S1C), by anchoring one pole of the mitotic spindle, prohibits divisions from modifying existing network edges (Diegmiller et al., 2022; McGrail and Hays, 1997; Storto and King, 1989; Telfer, 1975). See the Discussion for a proposed model extension to handle the case of divisions into existing network edges.

For assumption 6, asymmetric divisions can occur, leading to persistent differences in the oscillations of directly related cells. Although further experimental work is required, one hypothesized mediator of the asymmetry is unequal division of fusome volume, which plausibly leads to different effective concentrations of fusome-associated cell cycle regulators (such as deubiquitylases) (Diegmiller et al., 2023; de Cuevas and Spradling, 1998; Huynh and St Johnston, 2004). Our assumption that the oscillating variables are copied at division is consistent with experimental observations. For example, Cyclin A, the levels of which regulate the number of divisions, is diffuse in the cytoplasm at mitosis, so the difference in its concentration between the two resulting cells immediately after division is likely small (Lilly et al., 2000).

For assumption 7, the number of divisions is likely regulated by a pulse of a cell cycle regulator (Ji et al., 2017; Insco et al., 2009). Existing experimental evidence indicates that a complex involving two proteins, Bam and Ovarian tumor (Otu), acts to deubiquitylate Cyclin A. Furthermore, Bam is reported to have a ‘pulsed’ concentration (McKearin and Spradling, 1990; McKearin and Ohlstein, 1995; Insco et al., 2009) throughout cells of the network, increasing as the network starts to grow and decreasing as it stops, thus also changing Cyclin A degradation rates.

With respect to assumption 7, the oscillator variables x and y are affected by pulse z, but do not themselves appear in the differential equation for z. z follows its own prescribed dynamics. In agreement with this treatment, the initiation of the pulse of Bam depends on factors external to the cell (Chen and McKearin, 2003). Initial evidence also suggests a strong role for external factors, from some somatic cells surrounding the dividing cyst, in controlling the last divisions (Shi et al., 2021). Thus, our treatment of the pulsed variable is in accordance with existing experimental evidence.

There are five parameters of the single-cell relaxation oscillator. Three (ε, γ, z) appear directly in Eqn 2, and two control the shape of g(x) as explained in the next section. In total, the growing network model has eight free parameters: z becomes a function of two parameters (v and t_p, see below) and we also introduce the cell coupling parameter c and division asymmetry δ. In this work, we focused on the multicellular model and thus studied variation of v, c and δ.

From the intracellular dynamics in Eqn 2, note that the x- and y-nullclines are defined by y = g(x) and ⁠, respectively. When the nullclines intersect on the left (L) or right (R) branch of g(x), the system tends to the associated stable fixed point. Conversely, when the nullclines intersect on the middle branch of g(x), the fixed point is unstable and oscillatory behavior results. Thus, the corners of g(x), and can be used to define entrance and exit from the oscillatory region (annotated in Fig. 2A).

For the minimal single-cell oscillator considered here, we defined division events based on completed periods of its limit cycle. To detect these periods, we used repeated threshold crossing of the fast variable x. A complete cycle occurs when x crosses a threshold x_T from the left, then the right, then from the left once more in sequence. We used ⁠, which is the center of the oscillation region defined by g(x) (depicted by the vertical dashed lines in Fig. 2A).

We developed code to simulate the network formation process outlined in Fig. 2. The code was written in Python 3.8 and relies on various libraries including NumPy and SciPy.

Our approach included vectorized routines to numerically solve the (growing) system of ordinary differential equations (Eqn 3). In brief, given an initial network configuration [e.g. a single cell in state x⁽⁰⁾ at time t₀], we integrated the system for a fixed time window Δt. We then inspected the full trajectory for completed intracellular oscillation events. There were two cases: (1) no completed oscillations were detected; or (2) at least one oscillation was detected. If t* denotes the completion of the first oscillation that has not yet been detected (note that t₀ ≤ t* ≤ t₀+Δt), we discarded the portion of the trajectory after t* and expanded the graph by initializing a daughter cell (introducing N additional dynamic variables). The expanded network at time t* was treated as an initial condition for the next iteration.

We repeated this process of incrementally extending the network trajectory by Δt and inspecting for oscillations until the network stopped growing, which was detectable by x_a → 0 for all cells a ∈ {1, …, M}. The code included numerous other functionalities, such as graph isomorphism checks to compare the growing network against a library of named graphs (e.g. those of Fig. 1C-F).

To identify natural networks produced by the model, we first simulated network formation for each point in a three-dimensional grid of 100×160×181 parameter sets (v, c, δ) ∈ [0.005, 0.017]×[0, 8]×[−0.09, 0.09]. Consolidated data from this grid appear in Fig. 4B of the main text (explained further in the next subsection). The large dataset also formed the basis for individual network trajectories and one-dimensional parameter scans shown in other figures.

Specific parameter values used to generate different items in the text are listed in Table S2. Parameter values used throughout the text are: ⁠, m = 2, γ = 10⁻² and ε = 10⁻². Furthermore, all networks were initialized as a single cell in the state (x⁽⁰⁾, y⁽⁰⁾, z⁽⁰⁾) = (0, 0, 0).

Our goal was to identify region(s) in the grid introduced above that produced networks that were isomorphic to one of the reported naturally occurring ones. To that end, {A₁, …, A_S} denotes a set of S target networks (e.g. those in Fig. 1C-F), and is the final network produced at point (v_i, c_j, δ_k) in parameter space. From these, we constructed a tensor Λ_i,j,k, where Λ_i,j,k=p if is graph isomorphic to A_p, and Λ_i,j,k=0 if it does not match any target network.

For each target network A_p, we then used regionprops3 (MATLAB, voxel connectivity of 26) to identify the connected component(s) in our tensor (this is how we defined ‘regions’ of parameter space). We disregarded any component that contained only one voxel. To check if two regions generating distinct target networks (A_p, A_q) were in contact, we again used regionprops3 (as above).

Although the analysis outlined above often identified more than one region that generated a given target network A_p, in Fig. 4B, we plotted a single point representing a particular region. See Fig. S4 for an example of how the underlying regions appear for D. melanogaster and C. perla. In most cases, including that of D. melanogaster, the point corresponds to the centroid of the largest corresponding region. For each of the networks contacting D. melanogaster in parameter space (see Fig. 4A), we specifically chose a representative region contacting D. melanogaster. Finally, the size of each point in Fig. 4B corresponds logarithmically to the voxel count of the region it represents.

We thank Jasmin Imran Alsous, Rocky Diegmiller, Trudi Schüpbach, John Tyson, Elizabeth Ables, Allison Beachum and Steven Strogatz for helpful discussion related to the model and/or manuscript, and Lucy Reading-Ikkanda for illustrations. The computations in this work were, in part, run at facilities supported by the Scientific Computing Core at the Flatiron Institute, a division of the Simons Foundation.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.202187.reviewer-comments.pdf