Introduction

Development is the process by which the intricate complexity of multicellular organisms is constructed from a single fertilized egg or some simple vegetative structure (Fusco and Minelli, 2023). Not many other natural processes lead to so much complexity in such a short time. Development entails a natural process of pattern formation: specific cells and cell types end up in specific positions in space (i.e., anatomy). This process of pattern formation can be seen as the generation of spatial information (i.e., information about where each cell and cell type are) from previous information (e.g., information within the fertilized egg). This latter information includes the DNA but also spatial information in the form of compartments with different proteins and RNAs in different spatial locations within the oocyte. In most species, development cannot proceed if this latter spatial information is removed experimentally (Gilbert and Barresi, 2023). This initial spatial information within the oocyte arises from spatial asymmetries in the mother’s gonads or from the environment (Gilbert and Barresi, 2023). In this sense, pattern formation in development does not usually start from a spatially homogeneous initial condition but from an initial condition that has some simple spatial heterogeneities. In that sense, we use the term pattern transformation instead of pattern formation (Salazar-Ciudad et al., 2003).

Development can be seen as a sequence of transformations between initial developmental patterns and latter developmental patterns (what we call resulting patterns) over developmental time. By developmental pattern, or simply pattern, we mean a specific distribution of cell types in space or a specific distribution of gene product concentrations over space. For example, the earliest such patterns would be the zygote and the latest the adult. In most animals, early patterns arise from the division of the fertilized egg into different cells that, thus, inherit different parts of the fertilized egg and different proteins and RNAs (Gilbert and Barresi, 2023). Some of these latter molecules act as transcription factors that lead to the expression of further genes, i.e., the transcription and translation of genes and the eventual synthesis of their gene products (Gilbert and Barresi, 2023). Other gene products are secreted and diffuse in the extracellular space and bind to specific receptors in distant cells (Gilbert and Barresi, 2023). Here, we call these molecules extracellular signals but in the literature they are also called morphogens, growth factors, paracrine factors, etc. (Gilbert and Barresi, 2023).

Cells can respond to extracellular signals by changing gene expression and, often, by secreting additional extracellular signals and regulating cell behaviors such as cell division, cell contraction, cell adhesion, cell death, etc. (Salazar-Ciudad et al. 2003; Gilbert and Barresi, 2023). The former type of response leads to further changes in gene expression over cells (i.e., further pattern transformations), while the latter leads to cell movement and, consequently, to changes in the distribution of cells and gene expression over space (i.e., further pattern transformations) (Salazar-Ciudad et al. 2003).

How cells respond to extracellular signals depends on the signal receptors and signal transduction pathways they express (Gilbert and Barresi, 2023). Signal transduction pathways are actually networks of molecular interactions that integrate incoming signals to determine cell responses. Here, we use the term gene network to refer to these networks of interactions, even if they also include molecules that are not gene products (Gilbert and Barresi, 2023).

A fundamental question in developmental biology is how pattern transformations occur. In animals, pattern transformation involves gene networks, signaling by extracellular signals and the regulation of cell behaviors (e.g. cell division, cell contraction, cell adhesion) and the mechanical properties of cells and tissues (Salazar-Ciudad, et al., 2003). As a preliminary step to understand pattern transformation by these processes, we restrict ourselves to the pattern transformations occurring through gene networks and extracellular signaling alone (not including membrane-tethered signals or mechano-transduction). We ask:

1. Which are the gene networks topologies that can lead to pattern transformation?

2. Can we identify all these topologies and classify them into classes leading to similar pattern transformations?

Among pattern transformations we are only interested in those that are non-trivial. By non-trivial pattern transformations we mean transformations in which:

P1. The resulting pattern is stationary and heterogeneous in space.

P2. There is at least one gene product that, in the resulting pattern, has a new spatial distribution of concentration maxima and minima over space (i.e. critical points of the concentration distribution over space, see S1 of the SI for details). By new we mean that no gene product had this spatial distribution of concentration maxima and minima in the initial pattern (see Fig.1). Thus, non-trivial pattern transformations include the emergence of new maxima or minima, the disappearance of existing ones, or their spatial displacement from one location to another. This excludes the simple widening of already existing concentration maxima, since, in this case, the position of the maxima does not change (see S1 in SI for further details). We only consider resulting developmental patterns that are stable in time. We say that a gene product has been patterned if the spatial distribution of its maxima and minima is new (i.e. not present in any gene product in the initial pattern).

Figure 1:

There are several previous theoretical studies exploring how gene products can be wired to lead to pattern transformations through cell signaling (Salazar-Ciudad et al., 2000, 2001; Cotterel and Sharpe, 2010; Marcon et al., 2016; Zheng et al., 2016; Jimenez et al., 2017; Diego et al., 2018 ; Leyshon et al., 2021). Some of these studies are restricted to networks of three gene products (Cotterel and Sharpe, 2010; Zheng et al., 2016) while others are restricted to studying a specific class of gene network topology (Marcon et al., 2016; Zheng et al., 2016; Jimenez et al., 2017; Rand et al., 2021; Leyshon et al., 2021). In a previous study, we addressed the same questions considered here and identified two classes of gene network topologies capable of pattern transformation (Salazar-Ciudad, 2000). However, this previous study did not show why the identified topology classes are the only possible ones. In addition, due to its purely numerical approach, this previous study failed to identify one class of gene network topologies capable of non-trivial pattern formation. Here we take a more general analytical mathematical approach to identify and characterize all possible classes, regardless of the number of gene products they entail.

It is worth noting that there is an abundant theoretical literature on the topic of gene networks in cell biology (e.g., metabolism, gene regulation for basic cellular functions, etc.). However, the topic of this article is fundamentally different from those. Here we are not dealing with gene networks within a single cell, but with systems of cells with identical gene networks that are coupled through extracellular signaling. Specifically, we are interested in those genes networks leading to non-trivial pattern transformation (which is an inherently spatial question).

In this study, we consider three broad types of initial patterns (i.e., initial conditions): homogeneous-with-noise initial patterns, spike initial patterns and combined spike-homogeneous initial patterns. In a homogeneous-with-noise initial pattern, some gene products have zero concentration everywhere while others have the same non-zero concentration everywhere except for some small random fluctuations due to noise at the molecular level. In a spike initial pattern, the concentration of all gene products is zero everywhere except in a region (i.e. the spike region made of a small number of cells) where some gene products have the same non-zero concentration. In the combined spike-homogeneous initial pattern, some gene products have zero concentration everywhere while others have the same non-zero concentration everywhere except in the spike, where the concentration is larger (see Fig. 2). Any other initial pattern can be constructed by combining spikes of different heights (i.e. different gene product concentrations) at different positions. We then study how gene product interactions can be wired into networks with different topologies to transform these initial patterns into other (i.e. resulting patterns) non-trivial ones.

Figure 2:

It is worth noting that these three basic initial patterns correspond to spatially discontinuous functions: in homogeneous-with-noise initial patterns, white noise is discontinuous by definition; in spike and combined spike-homogeneous initial patterns, there is a concentration discontinuity between cells on the edge of the spike and nearby cells outside the spike. However, once extracellular signal diffusion begins, these sharp boundaries are smoothed into differentiable gradients, where critical points can be properly defined (e.g., at the center of the initial spike).

Methods: The model

This article considers a set of N_c cells, each with an identical gene network of N_g gene products. For the purpose of clarity we explain our results as if cells have simple spatial arrangements (e.g., a line in 1D, a square lattice in 2D, etc.) and we later discuss them for more complex arrangements of cells in space. We refer to such spatial arrangement as the domain of our system. We consider two types of gene products: intracellular gene products and extracellular signals. The extracellular signals are gene products that are secreted by cells and diffuse in the extracellular space while intracellular gene products are gene products that are not. Even though we only use the general term gene product, our model applies, in principle, to any molecule whose concentration changes as a consequence of the concentration of other molecules in cells.

We consider that extracellular signals diffuse in the extracellular space immediately apical to cells and, therefore, that diffusion occurs in a space with the same spatial dimension than the arrangement of cells. In a 2D lattice of cells, for example, we consider that diffusion occurs in a 2D extracellular space immediately apical to cells. We do not consider any spatial information at scales smaller than a cell (i.e. cells are considered as points). Thus, we treat the extracellular space as an arrangement of contiguous regions, one per cell, and a single concentration value per gene product and region (i.e. concentration variations within each region are not considered, see Fig.3). We consider systems made of many small cells in which diffusion can be approximated by classical continuous equations.

Figure 3:

For intracellular gene products, g_i(t, x) denotes the concentration of gene product I inside cell at position x at time t≥0. For extracellular signals, g_i(t, x) denotes the concentration of gene product i (i.e., signal) in the region of the extracellular space immediately apical to the cell at position x.

The dynamics of gene product concentrations in each cell in our model obey the following system of N_g reaction-diffusion equations (Murray 2002),

With initial pattern (i.e. initial condition),

And boundary conditions,

where ∂_t is the time derivative ; x is space in (n=1,2,3 dimensions) ; is the vector of gene product concentrations (from gene 1 to N_g) at the cell at time t and position x. is a vector of powers of the gene product concentrations and m is a vector of positive integer exponents; governs the interactions between gene products and determines how much of each gene product i is being synthesized in a given instant of time. f is, thus, a function that takes a vector of concentrations as input and returns another vector of concentrations as output ; ∇² denotes the Laplace operator (i.e., in in 2D; in in 3D). is an N_g × N_g -diagonal matrix with real positive diffusion coefficients for each gene product along its diagonal. By definition extracellular signals have d_i>0 while for intracellular gene products d_i=0. Similarly, is a N_g × N_g-diagonal matrix with the real degradation coefficients for each gene product along its diagonal, μ_i≥0. The degradation coefficient of a gene product is the default rate at which it would be degraded by proteases if no other factors affect its degradation. Notice that equation (1) has N_g components and applies to each cell. For the cases of extracellular signals, i.e. their concentration in a cell can affect that in distant cells through the diffusion term (last term in (1)).

In (2), g₀(x) represents the initial pattern (i.e. initial spatial condition). In (3), are the boundary conditions of the system (zero flux boundary conditions). ∇ denotes the nabla operator (i.e., ∂_x g(t, x) in 1D; (∂_x +∂_y) g(t, x) in 2D; and (∂_x +∂_y +∂_z) g(t, x) in 3D.); and is a unit vector normal to the boundary of the cell arrangement.

The last term in (1) is Fick’s second law of diffusion (Fick, 1855) describing the contribution of extracellular diffusion to changes in the concentration of extracellular signals. Our model does not consider membrane-bound signals (see Salazar-Ciudad et al., 2000 for a similar study that does).

The first two terms in system (1) can be understood as a reaction term. f defines a directed graph via its Jacobian matrix J (i.e., the matrix of first derivatives of f with respect to each gene product concentration), where each node stands for a given gene product, and each edge describes an interaction between gene products (Fig.2B).

Following the standard terminology in molecular biology, we say that a gene is being expressed in a cell at a give time if its gene product is being produced by that cell (i.e. the first term of (1) is larger than zero for that cell at that time). Notice that, because of the diffusion term, extracellular signals can have positive non-zero concentration around cells that are not expressing them.

We say that gene product j directly regulates gene product k if the corresponding element of the Jacobian matrix of f is non-zero (i.e., J_kj ≠ 0). This regulation is positive if J_kj>0 (i.e., g_j increases with g_k); or negative if J_kj<0 (i.e., g_j decreases with g_k). Similarly, we say that k is downstream of j if there exists a chain of regulations between gene products going from j to k. Correspondingly, j is said to be upstream from k. Notice that a gene product can be both upstream and downstream of a set of gene products, thus forming a regulatory loop. We call extracellular loops those regulatory loops in which at least one gene product is an extracellular signal.

Our results are only valid for functions f and reactions terms (i.e., f (g)− Mg^m) that satisfy some specific requirements. These are very broad biological requirements that are likely to be fulfilled by many developmental gene networks.

Requirements on the reaction term and f

R1

f is continuous and continuously differentiable (at least locally), so that its Jacobian matrix can be properly defined. Within each cell, we assume that the number of copies of each gene product and its rate of production are large enough for its concentration to be treated as continuous, as in many other models of pattern formation (Salazar-Ciudad et al., 2000, 2001; Cotterel and Sharpe, 2010; Marcon et al., 2016; Zheng et al., 2016; Jimenez et al., 2017; Diego et al., 2018 ; Leyshon et al., 2021). Notice that the degradation term in (1) is also continuous and continuously differentiable, − Mg^m (t, x).

R2

f is such that f (g(t, x)) − Mg^m (t, x) is non-linear in at least one of its components. In other words, at least for one gene product, the dependence of concentration on that of other gene products (or its own) is not linear. This is because it is well known that no non-trivial pattern transformation is possible in purely linear systems (Kondepudi & Prigogine, 2014; Murray, 2002).

R3

f is an explicit function of g, but not an explicit function of time or any time derivative of g. This means that the regulation of a gene product at a given moment depends only on the concentration of other gene products at that given moment, and not explicitly on past concentrations or their rate of change over time.

R4

f is such that g is always non-negative and bounded. This is because g is a vector of concentrations and concentrations cannot be negative. Similarly, cells are finite and, hence, cannot produce an infinite amount of gene products. Thus, f should be such that gene product concentrations never go to infinity.

R5

f is monotonously increasing or decreasing with respect to each gene product concentration. In other words, the sign of gene products interactions cannot change with the concentration of the interacting gene products (when everything else is kept constant). Thus, if a gene product k activates another gene product j, then k is an activator of j, even if its concentration changes. The same applies to inhibitory interactions (i.e. the sign of the regulation of k by j does not change with the concentration of j). In other words, the partial derivative of f with respect to each gene product j, should have the same sign for any value of g_j,

where f_k is the output of f regarding gene product k (i.e., the k-th component of f), and all gene product concentrations, except for g_j, are kept constant for the evaluation of J_kj. There are some known examples of particular gene product interactions that do not fulfill this requirement (Nelson et al., 2021) but, in general, the sign of the regulation of a gene product by another does not change in any complex way with the concentration of either gene product. Even in the reported cases where this happens, sign changes only occur once or twice (Nelson et al., 2021). Requirement R5 is simply a way to consider that pattern formation does not just arise from very simple gene networks with very complex interactions between each pair of genes (i.e. changes in the sign of the regulation depending on the concentration of the regulator) but from networks of relatively simple interactions between gene products.

Through this article when we refer to the “the parameters” we mean the parameters described in (1) (i.e. its parameters, M, D and J) and to any additional parameters a specific f may have. We consider that these parameters do not vary over time or space.

In this article we define the topology matrix T of a gene network as the matrix whose elements are −1, 1 or 0 depending on the sign of the corresponding entry in the Jacobian matrix (i.e., T =sgn(J)) (see Fig.2B). The topology of a given gene network is then a description of which gene products interact with which gene products and the sign of that interaction. In contrast, from now on, we will reserve the term gene network for a gene network topology with a specific f and all its parameter values so that a specific resulting pattern can arise from a specific initial pattern (i.e. a full implementation of equations (1)-(3) for specific f and specific values of the parameters). In this article a topological class is a set of topologies that fulfill some requirement definable at the level of topology, e.g. they have at least one positive loop. Studying gene network topology is important because for most developmental systems we only have information about T, and not so much about f (Gilbert and Barresi, 2023).

The results are organized into a set of logical arguments (for the least evident of them, we provide mathematical proofs in the Supporting Information). First we present some trivial requirements that the topology of all pattern-transforming gene networks need to fulfill. Secondly, we explain why all gene networks leading to non-trivial pattern transformations can be classified into three topological classes, and their combinations. Thirdly, we explain how the topological class to which a gene network belongs determines which non-trivial pattern transformations it can produce and which ones it cannot produce. In other words, belonging to one of these three classes is a necessary, but not sufficient, condition for a gene network to be able to lead to non-trivial pattern transformations. Among the types of pattern transformations possible from its topological class, a given gene network would produce one or another depending on the exact form of f and the value of the parameters. However, the topological class to which a gene network belongs indicates which non-trivial pattern transformations it cannot produce, regardless of its parameters and f.

Results

Any network with more than two gene products can be partitioned into a set of subnetworks but, as we will see, partitioning gene networks based on extracellular signals is especially useful. In this article, we call a signal subnetwork to the set of interactions between all gene products downstream of an extracellular signal (including the signal itself). Signal subnetworks, or from now on simply subnetworks, can partially overlap since a gene product can be downstream of several extracellular signals (see Fig.4). An extracellular signal can be upstream of other extracellular signals and, thus, a signal subnetwork can include other signal subnetworks (see Fig.4).

Figure 4:

Basic requirements on gene networks capable of pattern transformation

In this section we present some simple requirements on the topology of any gene network capable of pattern transformation. These are rather trivial but introducing them here facilitates explaining later results.

I1

Any gene network capable of non-trivial pattern transformation must include at least one signal subnetwork whose extracellular signal is either present in the initial pattern or positively downstream of some gene product that is. By ‘present in the initial pattern’, we mean that the corresponding gene product has a non-zero initial concentration somewhere in the domain (i.e., its gene is expressed in the initial pattern). Since we do not consider cell movement, mechano-transduction, and membrane-tethered signaling, non-trivial pattern transformations can only occur through the secretion of extracellular signals, which necessarily requires the activation of at least one signal subnetwork.

I2

All gene products are being degraded at all times (second term in (1)). Thus, for a gene product to be present in the resulting pattern, it should receive positive regulation from some gene products. This ultimately requires that all the gene products present in the resulting pattern are either within a self-activatory loop or downstream of it (otherwise, their concentration will decay to zero over time). This broad case includes genes that are constitutively expressed, meaning that they are always produced by cells regardless of other gene products (this can be viewed as these genes being within a self-activatory loop of their own or downstream of one).

Considering I1 and I2, it follows that all gene networks capable of non-trivial pattern transformation should contain a positive regulatory loop that is upstream of the patterned gene products and downstream of gene products present in the initial pattern. These loops can be extracellular or intracellular (i.e., include extracellular signals or not).

Gene network classification

In this section we explain that all conceivable signal subnetworks can be classified into just three classes based on how cells respond to extracellular signals in terms of secreting, or not, the same extracellular signals (see Fig.4). Let A be the most upstream extracellular signal of a given signal subnetwork. Then, there are trivially only two options, either A is downstream of itself (directly or indirectly), or it is not. In the latter case, we say that the subnetwork is hierarchical, or class H, while in the former case we say that the subnetwork is emergent (Salazar-Ciudad et al., 2000). Notice that this exhausts all possibilities: an extracellular signal is either downstream of itself or not, there are no other possibilities. Thus, all signal subnetworks can be classified as hierarchical, emergent or a combination of them (for subnetworks composed of other subnetworks).

The emergent class can be further divided into the L⁺ class, if A is positively downstream (and upstream) of itself, or into the class L^-, if A is negatively downstream (and upstream) of itself. This gives three topological classes: H, L^- and L⁺.

H subnetworks can contain regulatory loops as long as these do not include an extracellular signal (because in this case the loop would be extracellular, and the subnetwork, emergent). In fact, because of I2, any H network leading to pattern transformation must contain a positive intracellular loop (we label these as l⁺).

A gene network topology can be composed of several signal subnetworks. We classify whole gene networks topologies in the same way as their composing subnetworks. For example, a gene network topology with only H subnetworks is an H gene network, while a gene network topology with only L⁺ subnetworks is an L⁺ network. Composite gene networks topologies are labeled according to their composing subnetworks, independently of the number of subnetworks of each class they contain (see Fig.5). For example, an H L⁺ gene network can contain one, or more, of H subnetworks, and one, or more, L⁺ subnetworks.

Figure 5:

Our classification into classes H, L⁺ and L^- exhausts all possible topologies for gene networks. Within each topological class, however, there are gene networks that can lead to non-trivial pattern transformations and gene networks that cannot (e.g. depending on the parameters and the exact f). In the rest of this article we first identify further topological requirements that gene networks of each class must satisfy in order to lead to non-trivial pattern transformations. Then, we study the types of resulting patterns possible from each topological class. These types are described in very broad terms. In other words, we do not specify which exact resulting patterns arise from each gene network but some general features and commonalities among the resulting patterns possible from the gene networks within each class. We do the same for the gene networks combining different classes of subnetworks.

Linear stability analysis

In this section, we introduce the mathematical apparatus that allows us to link the topology of a gene network to its ability to produce non-trivial pattern transformation: the dispersion relation. We will treat each initial pattern as a perturbation of an otherwise spatially homogeneous steady state, and we will study the conditions under which such perturbation grows in time. If initial perturbations do not grow, one recovers the original homogeneous steady state and, thus, there is no pattern transformation. Even when these perturbations grow, one may recover a different homogeneous steady state as a resulting pattern, i.e. no non-trivial pattern transformation either. Thus, instability against perturbations is a necessary but no sufficient condition for pattern transformation.

To analyze whether a perturbation will grow we will use the standard linear stability analysis of the dynamical system defined by our model equations (1)-(3). In this study we focus in those f functions for which equation (1) to (3) admit a steady state homogeneous solution g(t, x) ≡ g* (i.e. a homogeneous steady state) given by f (g*)− M (g*)^m=0. Hence, let denote a small concentration perturbation from the homogeneous steady state g*. Then, plugging back into (1)-(3) shows that the evolution of such perturbation in space and time is given by the linearized reaction-diffusion equations,

where is the reaction-diffusion matrix in the absence of diffusion; and represents the initial perturbation. In the reaction diffusion matrix, is the Jacobian matrix of f (i.e. the matrix of derivatives of f (g) with respect to each gene product) evaluated at the steady state g*, and denotes component-by-component vector multiplication.

Each initial pattern can be described as a different small initial perturbation around an otherwise homogeneous steady state g*. In homogeneous-with-noise initial patterns, g*>0 and is small white noise. In spike initial patterns, g*=0 and only for positions x such that ||x−x_c||< L, where x_c is the center of the spike and L its width. In the combined spike-homogeneous initial pattern, g*>0 and, again, for positions x in the spike such that ||x−x_c||< L. In both cases, the initial concentration in the spike must be small enough for linearization (5) to hold (see S of the SI for the case of larger spikes).

Over time, the initial perturbation can either decay (i.e., as t →∞) or grow away from the original homogeneous steady state. In the former case, we say that is a stable perturbation; in the latter, that is an unstable perturbation. Since we are interested in the gene networks that can lead to pattern transformations through extracellular signaling, we are interested in the gene networks in which perturbations are unstable when there is extracellular signal diffusion (i.e., when some elements of D are non-zero).

Solutions to (5) can be found using the spectral decomposition of the Laplacian ∇² with zero-flux boundary conditions (Murray, 2003). In this sense, let W_k(x) be a zero-flux eigenmode of the Laplacian operator in our domain, that is, a non-trivial solution (i.e., W_k(x)≠0 for at least one x) of the following equations,

for some k∈ℝⁿ, with n=1,2 or 3. We denote by σ(∇²) the countable set of all k∈ℝⁿ for which equation (6) admits a non-trivial solution (Baker et al., 2008), and we say that k∈σ(∇²) are the wavenumbers of the corresponding eigenmode W_k(x). This is a slight abuse of notation based on the fact that, in simpler spatial domains (e.g., a 2D square), eigenmodes W_k(x) are trigonometric functions for which each entry in k counts the number of maxima and minima that W_k(x) displays along each spatial dimension.

Eigenmodes W_k(x) can be seen as a generalization of classical Fourier modes and hence, we can now consider the series expansion of in terms of W_k(x), that is,

where each is given by,

If we now plug series expansion (7) into system (5) and use the linearity of the equations to split the sum into terms for each k, we get,

and then, if we use (6) to solve the space derivatives, and eliminate alike terms from each side of the equations (i.e., W_k(x)), we get that the time evolution of each ĉ_k(t) is given by the equation,

where R(k²)=R(0)+k²D=J_f(g*)−m⊙M(g*)^m−1 is the so-called reaction-diffusion matrix. Solutions to equation (10) are given by,

where exp(R(k²)t) is a matrix exponential; and each is computed as in (8), changing by the corresponding initial perturbation .

At this point, we can plug (11) back into (7) and thus, we get that solutions to the linearized reaction-diffusion equations (5) can be written as a series expansion of the form,

However, for our linear stability analysis, we are only interested in the long-term behavior of (12), namely whether grows over time or eventually decays back to 0. In this sense, the long-term behavior of exp(R(k²)t) is dominated by the eigenvalue of R(k²) with the greatest real part and then, we can reduce expression (12) to,

where Λ_k∈ℂ is the eigenvalue of R(k²) with the greatest real part, and is the spectral projection of onto the corresponding eigenspace with eigenvalue Λ_k.

The different eigenvalues of R(k²) vary with k, and this variation is given by the so-called dispersion relation. By definition, the dispersion relation corresponds to the characteristic equation of the reaction-diffusion matrix R(k²) (Murray, 2003; Baker et al., 2008), that is,

where I denotes the N_g×N_g identity matrix. Then, Λ_k∈ℂ in the exponent of (13) corresponds to the solution of (14) with the greatest real part for each given k (i.e., Λ_k= λ_i such that Re(λ_i)≥Re(λ_j) for any other solution λ_j of (14)). Accordingly, we will refer to Λ_k as the principal branch of the dispersion relation.

The principal branch in of the dispersion relation, (13), determines whether a perturbation is stable or unstable. If Re(Λ_k)<0 for some k∈σ(∇²), then decays to 0 as t →∞ and we say that k (resp. W_k(x)) is a stable wavenumber (resp., eigenmode). Conversely, if Re(Λ_k)>0 for some k∈σ(∇²), then grows as t →∞ and we say that k is an unstable wavenumber (resp., eigenmode). Hence, perturbation is stable if all wavenumbers k∈σ(∇²) are stable, and is unstable if there exists at least one unstable wavenumber k∈σ(∇²) for which the principal branch has a positive real part. Note that, in accordance with requirement R4, the unbounded exponential growth of unstable eigenmodes in (13) will be halted, (13), eventually, by the higher order terms in (1) that we disregard in the linear approximation (6).

Each gene network has a , M and D matrices, and each such set of matrices defines a different reaction-diffusion matrix R(k²). This means that each gene network has and associated dispersion relation that is given by equation (14). In this sense, a gene network can only produce pattern transformations if its associated dispersion relation has at least one unstable wavenumber k∈σ(∇²) for which Re(Λ_k)>0. Moreover, each initial perturbation can be described as a series expansion in terms of wavemodes W_k(x) (e.g., a homogeneous-with-noise initial pattern will result in a series with many eigenmodes W_k(x) and many of them with large wavenumbers k). In this sense, each gene network can be seen as leading, or not, to pattern transformation by specifically destabilizing some of these wavenumbers while stabilizing others. The later means that some perturbations, the ones with the destabilized wavenumbers, grow over time until some non-linearities in f (R2 and R4) may preclude their further growth. This means that the resulting patterns, can be characterized, to a large extent, by these destabilized eigenmodes (Murray, 2003), even though some other wavenumbers may be destabilized by nonlinear effects (e.g., resonant eigenmodes (Castelino et al., 2020)).

The principal branch (i.e. Λ_k) allows to classify all the pattern transformations arising from gene networks into two mutually exclusive types: those with a finite amount of unstable wavenumbers and those with an infinite amount of unstable wavenumbers.

We say that a gene network is RD-unstable of the first kind if, for some choice of the parameters, its associated dispersion relation yields a finite amount of unstable wavenumbers (see Fig.6A-B). In this case, the resulting patterns, if heterogeneous, are periodic since series expansion (13) reduces to a sum of finitely many unstable eigenmodes, and eigenmodes are trigonometric functions for the vast majority of simple 1, 2 and 3-dimensional domains considered for pattern transformations in development (Baker et al., 2008). Notice that, since none of the initial patterns we consider are periodic, the emergence of these periodic resulting patterns constitutes a non-trivial pattern transformation (i.e. new concentration peaks and valleys will arise in some places).

Figure 6:

We say that a gene network is RD-unstable of the second kind if, for some choice of the parameters, its associated dispersion relation yields an infinite amount of unstable wavenumbers (see Fig.6). Previous research (Klika et al., 2012) has shown that the dispersion relation never goes to positive infinity (see Fig.6E) and that gene networks that are RD-unstable of the second kind have a dispersion relation that saturates to a positive value and, thus, have an infinite queue of large unstable wavenumbers (see Fig.6C-D). Since there is an infinite number of unstable wavenumbers, the resulting patterns, if heterogeneous, can have an infinite sum of eigenmodes in (9). There exist both periodic and non periodic functions that admit such an infinite series representations (Stein & Shakarchi, 2003). Nevertheless, as we later explain, these gene networks lead to three broad types of heterogeneous resulting patterns: periodic (but not trigonometric) patterns, noisy patterns, and radially symmetric multi-peaked patterns. Which of these can arise depends on the initial pattern and the topological class of the gene network.

Positive regulatory loops determine the kind of RD-instability

For a gene network to be RD-unstable it needs to have one or more positive regulatory loops. If such loops are all extracellular, then the gene network can only be RD-unstable of the first kind or RD-stable ; and if such loops are all intracellular, then the gene network can only be RD-unstable of the second kind or RD-stable. In here, we provide a proof of these two statements for the particular case of gene networks with only one positive loop and linear degradation, while in S2 of the SI we give a general prove for gene networks with any number of loops.

In the square matrix in (13), each off-diagonal entry corresponds to the interaction between two different gene products, while each diagonal entry contains the variable λ, the degradation rate of a given gene product, its diffusion coefficient and, perhaps, some self-regulation term (i.e., J_ii≠0). By the Leibniz formula for determinants (Strang, 2016), each term in det (i.e., the characteristic polynomial of the reaction-diffusion matrix) is the product of N_g entries of different rows and columns (i.e. different permutations). This means that each term of the characteristic polynomial of R(k²) contains only one entry per gene product (i.e., row or column) and hence, only the terms that correspond to regulatory loops are non-zero (since all other terms contain at least one zero element in the product), except for the product of all diagonal entries. The characteristic polynomial reduces, thus, to a sum of loop terms. This implies that, if a gene network contains no loops, the characteristic polynomial of R(k²) has just one term: the product of the elements of its diagonal. Namely, the dispersion relation of the network reads,

and, thus, its eigenvalues are given by,

for all i=1, …, N_g. In this case, all eigenvalues are negative and, thus, all wavenumbers are stable: no pattern transformation is therefore possible. This proves that, without regulatory loops in the network, pattern transformations are not possible.

Let us now consider the case of a gene network with only one loop, which we assume to contain N_l>1 gene products. Without loss of generality, gene products can be relabeled so that the submatrix of of the loop contains all non-zero elements under the diagonal and in the top-right entry. Namely,

Similarly, given that there exists no other loop in the network, we can relabel the remaining gene products so that the rest of the matrix is upper triangular (Bang-Jensen, 2008). Thus, the characteristic polynomial has two terms: the product of all diagonal entries, and the product of the non-diagonal entries in the loop times the remaining N_g−N_l diagonal entries. The corresponding dispersion relation reads,

where the minus sign before each J_ij comes from (13) and the factor comes from the Leibniz formula for the case of a loop like (16). A (−1) factor can be extracted from each Jacobian entry to obtain a factor that then simplifies (17) to,

Similarly, the factor (λ+μ_i+k²d_i) is found N_g times in the left term and N_g−N_l times in the right term and so, (18) can be further simplified into,

All terms in the first product of (19) are positive and, thus, its expansion is an N_g-polynomial in λ with positive coefficients,

and then, the characteristic polynomial has a negative constant term if:

Notice condition (21) only holds if the regulatory loop is positive (i.e. all regulations are positive or there is an even number of negative ones). If condition (21) holds for some k, it follows that the characteristic polynomial takes a negative value for λ=0. Moreover, since the term with the largest degree is positive, the characteristic polynomial becomes positive as λ →∞. This implies that the polynomial has to cross the λ-axis at least at one positive λ and such positive root is a positive eigenvalue of R(k²) for any k satisfying (21). In other words, we have proven that the dispersion relation of the gene network with a positive regulatory loop can yield unstable wavenumbers and, thus, the network can be RD-unstable depending on the parameters.

If the positive regulatory loop is extracellular, then at least one of the diffusion coefficients of its gene products is non-zero (i.e., d_i≠0 for some i in (21)). Then, given any specific choice of the parameters (i.e., specific values for J_ij) condition (21) will not hold for large wavenumbers k because the right-hand side simply grows unbounded as k²→∞. This means that condition (21) can only hold, at most, for a finite amount of small wavenumbers k. Thus, gene networks that have a single positive loop and in which this loop is extracellular, can only be RD-unstable of the first kind, or not RD-unstable at all. A general proof for the case of gene networks with several positive extracellular loops, and in fact for any gene network topology, is given in S2 of the SI.

If the positive loop is intracellular, then the diffusion coefficients of all its gene products are zero (i.e., d_i=0 for all i in (21)). Consequently, condition (21) reduces to,

In this case, the existence of positive eigenvalues does not depend on any specific wavenumber k. This means that, whenever the parameters allow (22) to hold, all wavenumbers k are unstable (i.e., there are infinitely many unstable wavenumbers). In other words, gene network that have a single intracellular positive loop can only be RD-unstable of the second kind, or not RD-unstable. Again, a general proof for the case of gene networks with several positive intracellular loops is given in S2 of the SI.

As we have seen, at least one positive regulatory loop is a necessary for a gene network to lead to pattern transformation. As we will see, for non-trivial pattern transformations to occur, negative regulation is also required (otherwise, the resulting pattern can only be homogeneous). Given that positive loops are required, and that there are only three classes of signal subnetworks, it follows that the gene networks leading to pattern transformation can only be: gene networks with positive intracellular loops and intracellular negative regulation (i.e., H networks); gene networks with positive intracellular loops and negative extracellular loops (i.e., L^- gene networks with positive intracellular loops l⁺); gene networks with extracellular positive loops and extracellular negative loops (i.e., emergent L⁺L^- gene networks); or a combination of any of the above. As we will see, each of these classes of gene networks leads to qualitatively different types of resulting patterns and the dispersion relation, together with some topological considerations, can be used to shed some light into the types of resulting patterns possible.

Hierarchical gene networks

By definition, hierarchical networks do not have extracellular loops and, due to requirement I2, hierarchical gene networks leading to non-trivial pattern transformations include at least one positive intracellular loop. As we saw in the previous section, this means that these networks are RD-unstable of the second kind.

Pattern transformations from spike initial patterns in H networks with a single extracellular signal

Let us begin by considering H gene networks consisting of a single H subnetwork with a single extracellular signal. We call these H⁰ networks. From a spike initial pattern, it is only the cells in the spike that, by definition, express the extracellular signal. This implies that it is only these cells that are secreting the extracellular signal. Thus, the concentration of the extracellular signal in all other cells (i.e., cells outside the spike) is fully governed by signal diffusion from this source and its degradation. Indeed, if we consider the set of cells either on the right-hand side of the spike (computations are equivalent for the left-hand side), we can solve equation (1) assuming linear degradation in the signal to show that its concentration forms an exponential gradient around the spike,

where A is the extracellular signal in the H⁰ network, is its concentration in the resulting pattern, x is the distance to the border of the spike and L is half the width of the spike. In 2D, the signal concentration gradient decays as , where ||x||₂ denotes euclidean norm; and in 3D, as (Sommerfeld, 1949; Stakgold & Holst, 2011). In all cases signal concentration has a radially decaying spatial distribution centered on the spike. In S3 of the SI we calculate C.

In the spatial distribution arising from (23), each cell at each side of the spike is at a unique distance from the spike and, thus, experiences a unique concentration of A. In that sense, there is a unique correspondence between the concentration of A and the distance to the spike. This correspondence has been used in the literature (Wolpert, 1968) to propose that, with a single extracellular signal, any pattern transformation is possible (see Fig.S1 in the Supplementary Information). The only thing that would be required is that cells interpret the signal concentration in a different way for each different resulting pattern attained (Capek and Müller, 2019; Sharpe, 2019). However, what would that interpretation be, or which would be its underlying gene network, was not specified by this author (Horder, 2001). Allegedly, any interpretation should rely on specific gene networks. A common idea for such interpretation is that cells express specific genes if they receive A at concentrations beyond some threshold value that is different for different genes. Then, different genes may become expressed at different distances from the spike and, thus, a non-trivial pattern transformation would occur (Capek and Müller, 2019; Sharpe, 2019).

There are, however, many other ways in which H⁰ networks can lead to non-trivial pattern transformations. In the following, we explain some fundamental requirements that H⁰ networks need to fulfill in order to lead to non-trivial pattern transformations and which are their possible non-trivial pattern transformations. This new requirements should be added on top of requirements R1-R5 to enable non-trivial pattern transformations in H⁰ networks.

Requirement RH1

For non-trivial pattern transformations in H⁰ networks, at least one gene product that is positively downstream of the extracellular signal A (let us say gene product k) must inhibit at least one gene product that is also positively downstream of A (let us say gene product j).

Requirement (RH1) has been previously discussed by other authors (Munteanu et al., 2014), especially for networks with three gene products, but we include it and expanded it here for completeness (see Fig.7). This requirement can be split into two parts: 1) there needs to be at least two gene products that are positively downstream of A; and, 2) one of them has to inhibit the other. Part 1 is simply requirement I1 applied to H⁰ networks: any gene product undergoing a non-trivial pattern transformation must be positively downstream of an extracellular signal and in H⁰ networks, this is A because it is the only signal there is.

Figure 7:

To understand part 2, let us first consider H⁰ networks in which part 1 holds but part 2 does not i.e., j and k are both positively downstream of A, but j is not inhibited by k. In that case, according to requirement R5 on f, both k and j will increase their concentration when that of A increases, and decrease their concentration when that of A decreases. Since the concentration of A decreases away from the spike (23), it follows that the concentrations of k and j also decrease with the distance to the spike and thus, k and j have their concentration peak in the same place as A and, thus, no non-trivial pattern transformation occurs. In fact, the same will occur for any gene products that, in a H⁰ network, are positively, and only positively, downstream of A.

Let us now consider the case in which part 1 and part 2 of RH1 hold. In that case, there is at least one gene product, k, that is downstream from A and inhibits another gene product, j, that is also downstream from A. In this case, the concentration of j can become very small where that of k is large (i.e., in the spike) and thus, form a concentration valley centered in the spike. This valley will in turn be flanked by two concentration peaks of j at each side (see Fig. 8). In the 2D and 3D cases there would be a central basin and a ridge surrounding it (see Fig.8 and Fig.9). However, the emergence of these concentration peaks and valleys is only possible if function f fulfills at least one of the two following additional requirements (see Fig.7):

Figure 8:

Figure 9.

Requirement RH2a

function f is such that the proportion between the activations of k and j by A changes with the concentration of A (i.e., f_k(g_A, …)/ f_j(g_A, …) is not a constant function of g_A). In other words, the functional form of the activation of k by A (i.e., f_k(g_A, …)) has a different non-linear part than the functional form of the activation of j by A (i.e., f_j(g_A, …)).

Requirement RH2b

The inhibition of j by k is non-linear in the sense that it is disproportionately larger for high concentrations of k than for low concentrations of k.

Both requirements imply an inhibition of j that is disproportionately different for different values of g_A and, thus, for different distances to the spike. In RH2a, this disproportionality arises from the different non-linear parts of the activations of j and k by A, while in RH2b, it arises from an intrinsic non-linearity in the inhibition of j by k (see Fig.7).

To better understand these requirements lets consider what happens if neither RH2a nor RH2b hold and all relevant interactions are linear. In this case, gene products k and j have concentrations that are, in each cell, proportional to those of A (and, thus, proportional to each other). In this case, if k inhibits j linearly, then g_j decreases proportionally to g_k in each cell. Consequently, g_j remains proportional to A over space but it is simply expressed at a lower overall level (i.e., no non-trivial pattern transformation can occur).

If the inhibition of j by k is disproportionately higher when g_k is large i.e. non-linear inhibition, then g_j would be disproportionately smaller where g_k is large (i.e. around the spike) than where g_k is small (i.e. away from the spike). This can lead to the formation of a valley of g_j in and around the spike and, consequently, to the formation of two peaks of concentration around it (in two or higher dimensions that would be a basin and a ridge surrounding it, see Fig.9).

To illustrate these arguments lets us consider an example gene network and f in which signal A directly activates gene products j and k, and k directly inhibits j in such a way that j and k are regulated by no other gene product in the network. In H⁰ networks the concentration of A is given by equation (23) and the dynamics of j and k are given by (24).

Where we assume, for simplicity, that both j and k suffer linear degradation (i.e., m_j=m_k=1). Now, assume that regulations f_k(g_A) and f_j(g_A, g_k) are both linear. Then, system (24) reads,

where B_kA, B_jA, B_jk >0 are real positive constants. Given that we are only interested in stationary resulting patterns, we can take time derivatives in (25) to zero and then we find that the concentrations of k and j in the resulting pattern are given by,

where we have used the gradient shape of from (23). Then, as we can see in (26), both concentrations have negative exponential distributions in space with the same exponential decay but different heights (i.e., different leading coefficients). Moreover, equation (26) also shows that there is a single concentration peak for A, j and k, and that these peaks are all centered at the same position (i.e. the center of the spike) and, thus, there is no non-trivial pattern transformation. This proves that no non-trivial pattern transformations are possible if all regulations are linear (requirement R2).

Now, let us consider the case in which all regulations remain linear except for the inhibition of j by k (as requested by RH2b). Then, system (24) becomes,

where f_jk is a non-linear inhibition of j by k. Hence, taking time derivatives in (27) equal to zero, the resulting pattern concentrations of j by k are,

where we have used the gradient shape of from (23). For gene product k, we have the same result as before: a single concentration peak centered in the spike. Things are, however, different for the resulting pattern of gene product j. From requirement R5, we know that f_jk is monotonically increasing with g_A and hence, it decreases with x. The two terms in have a maximum at x=0 and then, since the second term is negative, the overall value of in (28) can be lower at x=0 than elsewhere. If this happens, a concentration valley of j can form at x=0, flanked by two concentration peaks of j at each side of it. This implies the formation of new concentration peaks and valleys, and thus, that non-trivial pattern transformations are possible (see Fig.7) when requirement (RH2b) holds.

Next, let us consider the case in which RH2a holds and RH2b does not. In this case, the activation of k by A at high concentrations of A is disproportionately larger than that of j by A. As a result, g_k would be disproportionately larger than g_j in, and close to, the spike and not away from it (i.e. g_k and g_j would decay with a different slope from the spike). Then the inhibition of j by k, even when linear, can lead to the formation of a valley of g_j in and around the spike and, consequently, to the formation of two concentration peaks away from the spike (in two and higher dimensions that would be a concentration ridge and basin, see Fig.9).

Let us illustrate these arguments with a specific example in which regulations f_k(g_A) and f_j(g_A, g_k) are non-linear, but the activations of k and j by A have a different non-linear part in their functional forms (as required by RH2a) In this case, system (24) reads,

where f_jA(g_A) is a non-linear activation of j by A that is different from that of f_k(g_A) (e.g., one is a quadratic monomial and the other, a cubic monomial). We maintain linear inhibition of j by k because we are considering the case in which RH2a holds but RH2b does not. As in previous cases, if we equal time derivatives in (29) to zero, we get that the stationary resulting patterns of j and k are given by.

where we have used (23). Requirement R5 ensures that f_k(g_A) and f_jA(g_A) are both monotonically increasing functions of g_A and consequently, monotonically decreasing functions of x. For gene product k, this means that its resulting pattern will have the same monotony as (including its unique maximum at x=0), but with a different slope. Similarly, the two terms in share that same monotony with a maximum at x=0 but different slopes (e.g., if and , then , with . Consequently, since the second term is negative, the overall value of in (30) can be lower at x=0 than it is in nearby regions and, if this happens, a concentration valley of j can form at x=0, flanked by two concentration peaks of j at each side of it. Like before, when and where these concentration peaks and valleys form depends on the specific choice of the parameters (again, the parameters choice must be such that for all x). However, if they form, then gene product j would undergo a non-trivial pattern transformation as two new concentration peaks emerge at positions where there was no previous concentration peak (see Fig.7). This proves that requirement (RH2b), together with the proper choice of the parameters, enables non-trivial pattern transformations in H⁰ networks.

It is worth mentioning that, although we used a very simple H⁰ network for the illustration of requirements RH2a and RH2b, both requirements apply even if the regulations involved are not direct (e.g., if k activates some intermediate gene product i and j activates j). In this case, the different non-linear part in the functional form of regulations can be in any intermediate regulation (e.g., in the activation that k exerts on i, or the one that i exerts on j).

Pattern transformations from spike initial patterns in general H networks

In general, H networks may have different extracellular signals and these may affect patterning in downstream gene products. These networks can also lead to non-trivial pattern transformations and they do it with a set of less restrictive requirements than H⁰ networks (see Fig.7C). In this case, inhibition does not need to occur between gene products downstream of the same extracellular signal A, but can occur between gene products downstream from different extracellular signals secreted from the spike, e.g. signals A₁ and A₂. Then, requirements RH1, RH2a and RH2b are no longer necessary as long as these different extracellular signals have different diffusion coefficients (i.e., ), or different degradation coefficients (i.e., ). Hence, former additional requirements shall be substituted by the following one:

Requirement RH2c

For non-trivial pattern transformations H networks with several extracellular signals, there must exists some gene product j that is positively downstream of one extracellular signal A₁ and negatively downstream of another extracellular signal A₂, both secreted from the spike and with different diffusion coefficients (or different degradation rates).

If this occurs, then according to (23), these two signals have concentration gradients with different slopes along x. In other words, some of the extracellular signals have a higher proportion of their concentration close to the spike than others. As in the previous case then, concentration valleys can form around the spike for those gene products that are inhibited by one signal and activated by the other, directly or indirectly (see Fig.7). To illustrate this possibility, let us consider a simple H network where extracellular signal A₁ directly activates an intracellular gene product j, while signal A₂ directly inhibits it. Thus, if we assume that regulation is linear, the concentration dynamics of j is given by,

where B_jA1,B_jA2>0 are positive real constants. Then, taking the time derivative in (31) equal to zero, we have that the emerging patter of j is given by,

Equation (32) is essentially the difference of two negative exponential functions. Then, if the slope of the first exponential is smaller than that of the second exponential, this subtraction can lead to the formation of concentration valley in at x=0 flanked by two concentration peaks at each side of it (thus leading to non-trivial pattern transformation).

Pattern transformations arising from the combination of multiple H subnetworks

Gene products that undergo a non-trivial pattern transformation through an H network can lead to further non-trivial pattern transformations in downstream gene products without the need for activating additional extracellular signals (Scalise & Schulman, 2014). This can occur in gene networks in which gene products with a heterogeneous resulting pattern activate a gene product that is, in turn, inhibited by gene products with another different pattern (see Fig. 10). The resulting pattern in the downstream gene product can then have concentration valleys where its inhibitory gene products have concentration peaks, and concentration peaks where its activating gene products have concentration valleys (we call this peak subtraction).

Figure 10.

Similarly, a gene product that is positively downstream of gene products with different patterns can also develop a new pattern in which concentrations peaks would form wherever any of the upstream gene products have concentration peaks, as long as the concentration peaks of upstream activators are sufficiently distant from one another, (see Fig.10). We call this peak addition.

In both cases, the downstream gene product may acquire new concentration peaks and valleys and, thus, peak addition and subtraction can lead to non-trivial pattern transformations. The total number of peaks in the resulting pattern is less than or equal to the sum of the number of peaks in the resulting pattern of upstream gene products (as explained in section S4 and S5 of SI).

The ensemble of possible pattern transformations from spike initial patterns in H networks

As we have seen in the previous subsection, H subnetworks can be combined in simple ways to lead to resulting patterns with many concentration peaks and valleys. In this subsection we argue that, with some broad restrictions, for any radially symmetric resulting pattern there is a H network producing it from a spike initial pattern.

Let us first acknowledge that any resulting pattern can be seen as a sequence of concentration peaks of different heights and steepness at different positions (or basins and ridges in higher dimensions). Now, let us consider a simple gene network that we call the diamond network (see Fig.11A). This gene network has an intracellular gene product in the spike (say gene product 1) that self-activates and activates two extracellular signals (say signals 2 and 3) that, in turn, regulate a downstream intracellular gene product (say gene product 4). Gene product 4 is activated by signal 2 and inhibited by signal 3. Finally, the gene network includes an intracellular gene product 5 that is activated by gene product 4.

Figure 11:

As we know from previous subsections, if extracellular signals 2 and 3 have different diffusion coefficients, or degradation rates, the resulting pattern of 4 can exhibit a pair of concentration peaks at a distance from the spike. With some very mild restrictions (specified in S4 of SI), this distance can be made larger or smaller by changes in these parameters and the rest of parameters (see S4 of the SI for a more formal discussion) (see Fig.11). The resulting pattern of gene product 5 will exhibit concentration peaks in the same place as gene product 4. However, by choosing an adequate f₅(g₄) (while keeping R1 to R5) the concentration peaks of gene product 5 can have different heights and steepness (i.e., different decrease rates from its peak). In the most extreme case (e.g., using a sigmoidal function for f₅(g₄)), the concentration of gene product 5 can be made to be only appreciable in the cell where it attains its peak and its most immediate neighboring cells (i.e. very narrow concentration peak).

This simple diamond gene network can be combined any number of times to produce any number of peaks in any combination of positions and heights: simply, gene product 1 would activate different diamond gene networks and each such networks would have different parameters that lead to concentration peaks in different locations and with different heights and steepness. Then, the peak addition described above (see Fig.10) can put together all such peaks in the resulting pattern of a gene product downstream of all such diamond networks. Hence, the resulting pattern can produce any array of concentration peaks and valleys (i.e. target resulting pattern) as long as it is symmetric around the spike, it is continuous over space, and peaks and valleys are sufficiently distant from each other (peaks that are very close together may merge; see S4 of SI). Moreover, it is worth noting that if the resulting pattern consists of very narrow concentration peaks, function f and its parameters may need to be tuned very precisely and thus, patterns made of extremely narrow peaks may not be biologically feasible in practice.

The hierarchical gene networks studied in the previous sections are just simple examples. Since these simple H gene networks can lead to nearly all resulting patterns, one can expect that there are other, more complex, H networks that can also do it.

Finally, it is worth noting that hierarchical gene networks acting on a spike initial pattern tend to have similar variational properties (as studied in Salazar-Ciudad et al., 2000, 2001). Thus, in many H networks, the properties of the different peaks and valleys in a pattern (i.e. their height, position and steepness) are able to vary in respect to those of many other peaks in a pattern (while keeping radial symmetry). In the combination of diamond networks just described, for example, each pair of peaks arises from a specific diamond subnetwork with specific extracellular signals and its own parameters. This implies that it is possible, as we have seen, to regulate the height, steepness, and position of each pair of peaks in the resulting patterns (or, of each ridge in higher dimensions) independently from each other through changes in the parameters of each subnetwork (see Fig. S2 for other simple examples). However, not in all H gene networks all peaks can vary independently from each other. Peak addition and subtraction, for example, can lead to some peaks varying together with parameter changes (Salazar-Ciudad et al., 2000).

Pattern transformations from combined spike-homogeneous initial patterns in H networks

The pattern transformations possible from the combined spike-homogeneous initial patterns are similar to those possible from the spike initial conditions and have the same properties (see Fig.8-9). It suffices to add steady-state concentration to the signal gradient in (23) and all subsequent arguments in the previous subsections follow similarly.

Pattern formations from homogeneous-with-noise initial patterns in H networks

As we have seen in previous section, for pattern transformations to be possible, the homogeneous steady state needs to be unstable to perturbations. This means that the gene network, f and the parameters must be such that equations (1) admit at least three equilibrium points (the initial steady state and some others equilibrium points).

Gene networks that admit only two such equilibrium points cannot lead to non-trivial pattern transformations. The reason is that the initial equilibrium point needs to be unstable for pattern transformation to occur and that noise is applied to all cells. If the initial equilibrium is indeed unstable, then noise will lead all cells to change from the initial unstable equilibrium point to the second equilibrium point. In other words, the homogeneous-with-noise initial pattern will transforms into another homogeneous pattern with just different concentrations of the gene products.

Gene networks that admit three or more equilibrium points can lead to pattern transformations. In this case, cells may move to one equilibrium point or another depending on the magnitude of the random perturbation they initially receive. For example, in the case of gene networks that admit three different equilibrium points such that for all gene product i, the gene product concentrations in each cell where the initial concentration is below may change to g^**. Similarly, any cell where the initial concentration is above may change to g^***. This way cells randomly settle in different equilibrium points (g^** and g^*** in the example).

If the above gene networks do not include any extracellular signal, the resulting pattern will simply be a random distribution of these equilibrium points among cells. This will be, in other words, a random distribution of concentration peaks and valleys. Each of these peaks and valleys will form where the initial pattern also had, due to noise, tiny peaks and valleys. In that sense, no new concentration peaks and valleys appear and, thus, without signaling, there are no non-trivial pattern transformation (not for H networks nor in any other gene networks).

In the case of H gene networks and homogeneous-with-noise initial patterns there are, by definition, extracellular signals and, thus, concentration peaks and valleys can form where the initial pattern had tiny concentration peaks and valleys (as explained in the previous paragraph) but also at a distance from them, as explained for H gene networks under other initial patterns. In this case, thus, new peaks and valleys can arise and non-trivial patterns transformations are possible. The resulting patterns, however, still consist of random distributions of peaks and valleys (basins and ridges in higher dimensions)(see Fig.8-9). This is, in fact, what the dispersion relation indicates. Hierarchical networks are RD-unstable of the second kind and, thus, have an infinite number of unstable wavenumbers. Except for sampling limitations, the initial pattern with white noise contains all wavenumbers. This means that many, if not most, wavenumbers present in the initial noise can be amplified and, thus, the resulting patterns most likely consist of random distributions of concentration peaks and valleys (see Fig.8-9).

Emergent gene networks

L⁺ gene networks

These are gene networks with one or several L⁺ subnetworks (and no other types of subnetworks). From a spike initial pattern, the signal, or signals, in the L⁺ subnetworks diffuse from the spike across the extracellular space and activate their own production in the cells that receive it. Since these signals promote their own production, the cells receiving them also end up producing and secreting them. As a result, each signal progressively spreads further until it ends up being expressed by all cells (see Fig.8-9). Due to the averaging effect of diffusion and the self-activating nature of L⁺, all cells end up having the same concentration of the signal (even noise is suppressed) and, thus, the only possible resulting patterns are homogeneous (i.e. no non-trivial pattern transformation occurs). The same happens for homogeneous-with-noise or combined spike-homogeneous initial patterns.

L^- gene networks

By definition, these networks only contain subnetworks with extracellular signals that inhibit their own production, directly or indirectly, in the cells that receive them. According to requirement I2, gene networks leading to non-trivial pattern transformations must include one or several positive loops (i.e. l⁺) positively upstream of the gene products present in the resulting pattern. For L^- gene networks, these loops have to be intracellular, otherwise the gene network would be L⁺L^- not L^-. Accordingly, we must only consider the L^- networks that include at least one L^- subnetwork downstream of at least one l⁺.

Pattern transformations in L^- subnetworks from spike initial patterns

In the spike initial pattern, only the cells in the spike initially express the genes of the L^- subnetwork. Since the secreted signal inhibits itself, the expression of this signal cannot spread beyond the spike (i.e. it will only be secreted from the spike) and, thus, its peak concentration will be in the spike itself (i.e. no non-trivial pattern transformation).

Pattern transformations in L^- subnetworks from homogeneous-with-noise initial patterns

Lets us first consider the gene networks with a single L^- subnetwork and a single l⁺. Among these, let us first consider the case of L^- gene networks with a single L^- subnetwork positively downstream, but not upstream, of a single l⁺. Since l⁺ is not affected by any extracellular signal, its gene products can only have a homogeneous resulting pattern, or a noisy resulting pattern like the one described for gene networks without extracellular signals (in the “Pattern formations from homogeneous-with-noise initial patterns in H networks” section). Neither case constitutes a non-trivial pattern transformation. By definition, the most upstream extracellular signal, A, inhibits itself. Then it can only be expressed in the same place than its positive regulators, l⁺, and thus, these gene networks cannot lead to non-trivial pattern transformations either. The case of a L^- subnetwork upstream of an l⁺ but not downstream of it does not fulfill requirement I2 and, thus, cannot lead to pattern transformations either.

The only case left is that of a L^- subnetwork that is both downstream and upstream of the l⁺. In the simplest case, l⁺ activates A and A inhibits some gene product in l⁺, and thus indirectly all of them and itself. We call these networks over-Turing gene networks since, as we will see, they are capable of producing periodic patterns that resemble those of classical Turing gene networks (see Fig. 8-9). Since the only positive loop is intracellular, these gene networks are, depending on the parameters, either RD-unstable of the second kind or RD-stable (see Fig.6C-D) (Klika et al., 2012).

In each cell, l⁺ activates A but A diffuses extracellularly and inhibits l⁺ in the cells that receive it and, thus indirectly itself. This means that, indirectly, the l⁺ in the different cells are competing with each other. In other words, by inhibiting the l⁺ of the other cells, the l⁺ in a cell is indirectly activating itself. If the system is made of just two cells, the concentration of the gene products in l⁺ (as well as A) can only increase in one cell (the one experiencing a weaker effective inhibition of A by l⁺ at the initial time) and decrease in the other (the one subject to a stronger initial inhibition of A by l⁺). This implies that the initial difference in the concentration of A between cells can only but grow until some steady state may be reached (in accordance with requirement R4). In other words, the initial concentration difference between cells (i.e. the initial noise in the pattern) becomes amplified.

The two-cell example can be used to understand larger systems with homogeneous-with-noise initial patterns. These systems can be understood as being composed of multiple overlapping pairs of cells that, due to noise, have small initial differences in the concentration of the gene products in l⁺ (or any gene product upstream of it). From the arguments explained in the previous paragraph, one can expect that the initial concentration differences between each pair of cells can only but grow. In the case of systems made of many cells, however, each cell is contiguous to more than one cell (e.g., two in 1D) and, thus, it can occur that two contiguous cells both increase their concentration of A because they are both contiguous to other cells that, by chance, have smaller initial concentrations of A. Which cells end up in concentration peaks or valleys depends on which cells had higher signal concentration in the initial pattern and the inhibition arising from the extracellular signals, as well as on the parameters of the network.

The resulting pattern consists of randomly distributed concentration peaks and valleys of A and the gene products in l⁺ (see Fig.8-9). These patterns are very similar to the ones that can arise from H networks acting on homogeneous-with-noise initial patterns (even if the underlying dynamics are quite different). This is, in fact, what the dispersion relation indicates in the case of RD-unstable networks of the second kind: high wavenumbers are unstable and thus, they all can be amplified to form a noisy resulting pattern. Indeed, the resulting pattern can be understood, to some extent, as an amplification of the initial noise in the homogeneous-with-noise pattern (Marcon et al., 2016; Diego et al., 2018). Because of this, each resulting pattern will be different, with different peaks and valleys in different positions and different distribution of heights between them (see Fig.S3A).

Let us now consider the ensemble of over-Turing gene networks. This is any gene network combining any number of l⁺ loops and L^- subnetworks downstream and upstream of each other (but no L⁺ subnetwork). In all of them it is still the case, by the definition of L^- subnetworks, that extracellular signals inhibit their own production in the cells that receive them and, that thus, cells are inhibiting each others’ production of each extracellular signal, directly or indirectly. This means that, as in the simple examples above, any initial difference in the initial concentration of extracellular signals between cells can only increase. In fact, over-Turing networks contain positive intracellular loops and, thus, irrespectively of their complexity, they can only be, depending on the parameters, RD-unstable of the second kind or RD-stable. In the former case, most wavenumbers would be unstable and, since there is an infinite number of them, the resulting patterns should include many of them, i.e. the resulting pattern will likely be a random distribution of peaks and valleys as in the case of the simplest over-Turing networks.

Pattern transformations in L^- subnetworks from spike-homogeneous initial patterns

As in the previous initial patterns, and for the same reasons, only the over-Turing gene networks can lead to non-trivial pattern transformations. Over-Turing gene networks are known to be able to lead to relatively complex and periodic non-trivial pattern transformations from the spike-homogeneous initial pattern (Marcon et al., 2016; Wang et al., 2022). By definition, the concentration of A starts being higher in the spike than elsewhere. As A diffuses, it inhibits itself in the cells around the spike and this leads to the formation of concentration valleys of A. However, A activates gene products in l⁺ and thus, similar concentration valleys of gene products in l⁺ will form around the spike (see Fig.S3B-C). Since l⁺ activates signal A, the cells in these valleys will produce less A and, in turn, cells beyond the valleys will receive a smaller amount of it. As a consequence, there is less inhibition of the gene products in l⁺ in such cells, and the concentration of A and the gene products in l⁺ will increase in those cells. This increase in the concentration of A leads to the formation of a concentration valley further away, and these valleys to further peaks, and so on. This process continues until the whole domain is occupied by a sequence of concentration peaks and valleys of A and the gene products in l⁺. These have very similar heights, shapes and distances between them (see Fig.S3), with some variation in those closer to the spike or the boundary of the domain. In 2D there are no concentration peaks but concentric rings of high and low gene product concentration centered around the spike. In 3D there are concentric spherical shells of low and high concentration. We refer to these patterns as frozen-wave patterns as they result from the action of a patterning wavefront that moves away of the spike leaving a stationary wave-like pattern behind. As in the case of H networks, these frozen-wave patterns are symmetric around the initial spike but, in contrast to H networks, they extend over the whole system.

The formation of these frozen-wave resulting patterns can also be understood as a consequence of the fact that, in over-Turing networks, any initial difference in the concentration of A between contiguous cells can only increase (as we explained in the previous subsection). In this case, the initial concentration differences are between the cells at the margin of the spike and the cells just outside of it in every radial direction. These differences can only but grow over time and thus, the concentration of A and gene products in l⁺ will decrease in cells just outside of the spike. This, in turn, leads to a concentration difference between the cells just outside the spike and cells further away, in which the concentration of A and gene products in l⁺ will increase again. This process continues leading to an alternation over space of concentration peaks and valleys (or concentric concentration rings and shells in higher dimensions) or, in other words, to the frozen-wave pattern.

The frozen-wave patterns cannot arise from spike initial patterns. This is because concentration valleys cannot form where gene product concentrations are zero (i.e. around the spike), since gene product concentration cannot be negative.

Changes in the parameters of over-Turing gene networks (e.g., diffusion coefficients) can lead to changes in the number, height or width of concentration peaks and valleys in the resulting pattern. However, contrary to what happens for H gene networks, concentration peaks and valleys in over-Turing networks cannot change independently from each other (see Fig.S4). This lack of variational independence in the properties of peaks and valleys is easy to understand from the fact that over-Turing gene networks produce many different concentration peaks from a small number of interactions in the gene network and thus, no specific part of the network is responsible for any specific concentration peak (as it can happen in H gene networks). This implies that the ensemble of pattern transformations attainable by these networks is less diverse than the ensemble of pattern transformations attainable by H gene networks. In other words, the over-Turing gene networks can produce patterns with equally sized and shaped peaks, while the H networks can also produce patterns with concentration peaks and valleys with different characteristics (different heights and shapes, Salazar-Ciudad et al., 2000, 2001).

Gene networks combining different classes of subnetworks

Subnetworks can be combined in parallel, in series, or both. By parallel combinations we mean that different subnetworks are upstream of the same set of gene products (while being downstream of some gene product in the initial pattern, requirement I1). Since in this case the different subnetworks do not regulate each other, the resulting patterns in downstream gene products can have concentration peaks and valleys wherever their upstream gene products have them (just as described for H gene networks above, see Fig.10).

By in series combinations we mean that some subnetworks are upstream of others (i.e., the signal of the first subnetwork is upstream of the signal of the other). In this case we are interested in the resulting patterns of the downstream subnetwork (since the pattern transformations in the upstream subnetworks are unaffected by the downstream gene network and have already been described in previous sections).

Subnetworks can be combined in series positively or negatively (i.e., the signal of the upstream subnetwork is negatively upstream of the signal of the other subnetwork). In the case of negative regulation, however, some gene products in the initial pattern must positively regulate the downstream subnetworks independently of the upstream subnetwork, otherwise the genes in the downstream subnetwork only receive negative regulation and, thus, cannot be expressed (i.e. requirement I1). In this case, the resulting pattern of the downstream network would be the one it would produce on its own, but with that of the other subnetwork subtracted from it, as we described for H subnetworks (see Fig. 10). The case of positive in series combinations of subnetworks is described in the next section. Naturally, several subnetworks can be combined in series to form a chain in which each subnetwork is upstream of the next one.

Subnetworks can also be combined in series but forming a loop (i.e., the signal of each subnetwork is upstream and downstream of the signals of other subnetworks). However, if an H subnetwork is combined forming a loop with any other class of subnetwork, then the signal of the H subnetwork is indirectly regulating itself and is, thus, an emergent subnetwork. In this sense, combining subnetworks in a loop leads to emergent subnetworks.

Pattern transformations in the combination of H and L⁺ and H and L^- subnetworks

As we have seen, all gene networks capable of pattern transformation can be RD-unstable of the first kind or RD-unstable of the second kind. Thus, regardless of how sub-networks are combined, there are only two main types of resulting patterns. H and L⁺ subnetworks can only be RD-stable or RD-unstable of the second kind (since the only positive loop is extracellular). Beyond that we can acquire a more detailed intuition of the possible pattern transformations through some simple considerations on the pattern transformations we have described for the H and L⁺ subnetworks individually.

Gene networks composed of L⁺ subnetworks upstream of H subnetworks are equivalent to H subnetworks acting on a homogeneous initial pattern without noise. This is because, as explained in the L⁺ section, these subnetworks lead to homogeneous resulting patterns whichever the initial pattern and, in the process, eliminate noise. Thus, the downstream H subnetworks effectively receive and input pattern that is homogeneous (i.e. without noise) and, thus, cannot lead to pattern transformations.

In gene networks in which H subnetworks are upstream of L⁺ subnetworks no pattern transformations are possible downstream of the H subnetworks because the L⁺ subnetworks homogenize any input patterns that the H subnetworks may convey.

Gene networks composed of L^- subnetworks upstream of H subnetworks can only lead to non-trivial pattern transformations in the case the L^- subnetworks are over-Turing, for the reasons explained in the L^- subnetworks section. The same applies to gene networks composed of H subnetworks upstream of L^- subnetworks. In this case the H subnetwork can be seen as providing the L^- subnetworks with a pattern but as we have seen L^- subnetworks cannot lead to non-trivial pattern transformations unless they are over-Turing, regardless of the initial pattern they act on.

In the case of H subnetworks upstream of over-Turing subnetworks and combined spike-homogeneous initial patterns, the frozen-wave pattern that would normally spread from the spike perturbation can spread, instead, from the peaks and valleys produced by the H subnetwork. Thus, the resulting is still a frozen-wave pattern but locally modified by the H subnetwork so that the peaks close to the initial spike can vary in size (this effect keeping radial symmetry around the spike) (see Fig.12). If the over-Turing subnetworks are upstream of H subnetworks, then the radially symmetric patterns that are attainable by the latter can be repeated around each concentration peak of the amplified noise pattern. When acting on a homogeneous-with-noise initial pattern, the only possible non-trivial pattern transformations are the noisy resulting patterns reported for both over-Turing and H networks (since both kinds of networks produce similar patterns on their own anyway).

Figure 12:

Pattern transformations in the combination of L⁺ and L^- subnetworks

In this section we first consider combinations between single L⁺ and single L^- subnetworks. If an L⁺ subnetwork is upstream of an L^- subnetwork, and not downstream of it (i.e., they form no loop), the L⁺ subnetwork is unaffected by the L^- subnetwork and consequently, its genes will be homogeneously expressed regardless of the initial pattern (see Fig. S5-6). This reduces to what we have seen is possible for L^- gene networks acting on homogeneous-with-noise initial patterns, except that the L⁺ subnetwork can eliminate the noise and preclude any pattern transformation. If the L^- subnetwork is upstream of the L⁺ subnetwork, then the resulting pattern of the latter subnetwork can only be homogeneous because, as we explained in previous subsections, L⁺ subnetworks produce these patterns whenever any of its signals is expressed in at least one cell (see Fig S5-6).

L⁺ and L^- subnetworks can be combined into loops. Previous research has extensively studied pattern transformations occurring from such networks (Cotterell & Sharpe, 2010), for the case of gene networks with up to three gene products. From these, the only resulting patterns and developmental dynamics that we have not discussed so far are those of Turing gene networks, also called Turing mechanisms (Turing, 1952; Murray, 2002; Maini et al., 2006; Meinhardt, 2008). These occur, for example, when a L⁺ subnetwork is positively upstream of an L^- subnetwork, and the latter is negatively upstream of the former (e.g., the signal of the latter inhibits, directly or indirectly, the signal of the former).

In Turing’s seminal work (Turing, 1952), each subnetwork is in fact a single extracellular signal, but it has been widely reported that these mechanisms also apply to gene networks where extracellular signals regulate each other indirectly, through an intracellular part of the gene network (Salazar-Ciudad, et al. 2000; Satnoianu, et al., 2000; Maini et al., 2006; Meinhardt, 2008; Cotterell & Sharpe, 2010).

Given that Turing gene networks contain an extracellular positive loop (i.e. an L⁺ subnetwork), they are, depending on the parameters, either RD-stable or RD-unstable of the first kind (see Fig. 6A-B). In the former case the resulting patterns are periodic. For homogeneous-with-noise initial patterns, a good approximation for the wavenumber of the resulting pattern is given by the unstable wavenumber closest to the maximum of the real part of the principal branch of the dispersion relation (Murray, 2002) (see Fig. 6A-B). For combined spike-homogeneous initial pattern a more detailed analysis is required to understand the wavenumbers in the resulting pattern (Tarumi & Mueller, 1989; Klika et al., 2024) but the resulting pattern is also periodic.

In Turing gene networks, the larger wavenumbers are stable (i.e. do not grow) because of the L⁺ subnetwork. The signal in the L⁺ subnetwork diffuses between cells and activates its own production in nearby cells. This produces a homogenizing effect in which, contrarily to what happens in over-Turing networks, the initial differences in concentration between contiguous cells do not necessarily grow over time. In fact, the positive extracellular loop precludes large concentration differences between nearby cells and thus, precludes large unstable wavenumbers (i.e. having many concentrations peaks per unit of space and so, large relative differences in concentration between contiguous cells).

In 1D, the possible resulting patterns of Turing networks are similar to the periodic patterns arising from over-Turing gene networks (see Fig. 8). The difference is that, while in Turing networks these periodic patterns can arise from both homogeneous-with-noise and combined spike-homogeneous initial patterns, in over-Turing networks they can only arise from combined spike-homogeneous initial patterns. In 2D and 3D, Turing networks acting on homogeneous-with-noise initial patterns can lead to resulting patterns consisting of several concentration peaks (e.g., dots) or ridges (e.g., stripes) with a constant spacing and regularly distributed over space, or forming labyrinths (Meinhardt, 1982; Murray, 2002), while over-Turing networks acting on those initial patterns can only produce noisy resulting patterns (see Fig.9). From combined spike-homogeneous initial patterns in 2D and 3D domains, over-Turing networks can lead to resulting patterns consisting of concentrically arranged concentration ridges and basins due to the radial symmetry of the problem (see Fig.9). This is also possible from Turing gene networks (Meinhardt, 1982).

As in the case of over-Turing gene networks, changes in the parameters of a Turing gene network (e.g., in diffusion coefficients) can change the global number, height or width of concentration peaks and valleys in the resulting pattern (Turing, 1952; Murray, 2002; Maini et al., 2006; Meinhardt, 2008). However, in contrast to what we see in H networks (Salazar-Ciudad et a., 2000), these parameter changes lead to all concentration peaks and valleys to change in the same way. The reasons for this difference are the same we discussed for over-Turing gene networks (see Fig.S4 and S5) (i.e. a relatively small number of gene products is responsible for the resulting pattern and no part of the network is responsible for any specific subset of concentration peaks or valleys). For this same reason the ensemble of distinct resulting patterns attainable by Turing gene networks is smaller and less diverse than the ensemble of distinct resulting patterns attainable by H gene networks. Simply, the ensemble of possible Turing networks can only produce resulting patterns where all concentration peaks (and valleys) are identical while the ensemble of possible H networks can produce resulting patterns where different concentrations peaks and valleys can have different heights, weights and spacings (see Salazar-Ciudad, et al. 2000, 2001 for a more detailed exploration).

Turing subnetworks can be combined into composite gene networks (see Fig.12). By combining Turing subnetworks in series, the resulting pattern of one subnetwork can serve as the initial pattern of the other (Fujita & Kawaguchi, 2013; Moustakas-Verho et al., 2014). As long as these combinations the only positive loops are extracellular, the system can only be, depending on the parameters, RD-stable or RD-unstable of the first kind (just like single Turing networks). Consequently, regardless of how complex the composite network is, the possible resulting patterns are still periodic patterns, that may combine several wavelengths (Fujita & Kawaguchi, 2013; Moustakas-Verho et al., 2014). Combining Turing gene networks with over-Turing gene networks in series leads to composite gene networks that behave either as RD-unstable networks of the first or second kind depending of model parameters (see section S2 in SI).

In general, gene networks that combine Turing gene networks and H gene networks in series have variational properties with features of both the Turing and H gene networks (see Fig.12). There are, in fact, many articles studying this (Salazar-Ciudad et al., 2001; Miura, 2013; Green & Sharpe, 2015; Glim et al., 2021; Tzika et al., 2023). If the Turing networks are upstream of the H networks, then each concentration peak resulting from the Turing network acts as a different spike in the input pattern received by the H network. Thus, the peaks and valleys resulting from the H network can be repeated around each of the concentration peaks produced by the Turing network. If the peaks produced by the H network are wider than those produced by the Turing, however, these latter peaks may fuse and the resulting pattern simply resembles that of the Turing network (see Fig.12).

In gene networks in which H networks are upstream of the Turing networks, the resulting pattern can be periodic, but the H networks can make the periodic pattern slightly different in different spatial regions (by specific gene product in the H networks regulating specific gene products in the Turing networks). Thus, the height or spacing of concentration peaks in the resulting pattern of the Turing network can be larger, or smaller, where patterned genes of the H network are expressed (see Fig.12).

Besides being combined into Turing or over-Turing subnetworks, or their combinations, L⁺ and L^- can be combined into any arrangement (e.g. into multiple nested loops). As we have shown in previous sections, however, a gene network with such complex combinations can only be, nonetheless, RD-unstable of the first or second kind. In the former case the resulting pattern will be periodic, as when combining Turing gene networks. In the second case, we do not know for sure which are the possible resulting patterns but since the number of unstable wavenumbers is infinite, we can expect that, from homogeneous patterns with noise, the resulting patterns are amplified noise patterns (just as combining several over-Turing networks). Similarly, from combined spike-homogeneous patterns we expect froze-wave patterns, possibly combining different periods.

Discussion

The main conclusions of this article are that all gene networks capable of non-trivial pattern transformations can be classified into three distinct topological classes (and their combinations) and that each of these classes allow for qualitatively different types of resulting patterns. For gene networks within each class we found some additional topological requirements.

All gene networks can be decomposed into subnetworks. Because we study pattern transformation by extracellular signaling, we decompose gene networks based on the subnetworks downstream of each of its extracellular signal. These subnetworks are classified based on whether their extracellular signals are positively downstream of themselves (L⁺ subnetworks), negatively downstream of themselves (L^- subnetworks) or not downstream of themselves (H subnetworks). This classification exhaust all possibilities, except for subnetworks that arise as a combination of several of these simpler subnetworks.

We have shown that pattern transformation requires gene networks with positive regulatory loops and that, if these loops are extracellular, then the resulting patterns, when heterogeneous, are periodic. We have also shown that, if the positive regulatory loops are intracellular, then the resulting patterns, when heterogeneous, are characterized by an infinite amount of unstable wavenumbers. However, negative interactions are also required for pattern transformation. The combination of the former findings with this latter requirement implies that there are only three types of gene networks leading to pattern transformations: 1) Hierarchical gene networks with positive intracellular loops and negative regulations 2) Over-Turing gene networks (i.e. gene networks with negative extracellular loops and positive intracellular loops forming a loop between them) 3) and Turing gene networks (i.e. gene networks with positive and negative extracellular loops forming a loop between themselves). These three fundamental classes, and their combinations, exhaust all possibilities for gene networks capable of non-trivial pattern transformations through extracellular signaling.

In this article we also show that these three classes of gene networks, and their combinations, each lead to specific types of resulting patterns (see Fig.8-9). Hierarchical networks lead to resulting patterns that are radially symmetric and made of any combination of concentration peaks and valleys (ridges and basins in higher dimensions). Depending on the parameters of the network, these resulting patterns can contain more or less peaks and valleys with different positions, heights and different shapes. In over-Turing gene networks infinitely-many wavenumbers are linearly unstable and thus, many get amplified in the resulting pattern. From homogeneous-with-noise initial patterns, this leads to amplified-noise patterns, while from combined spike-homogeneous patterns, this leads to radially symmetric periodic patterns extending all over the domain. Similarly, we also prove that Turing networks only amplify a finite number of wavenumbers, if any, and then, the only non-trivial pattern transformations they can lead to result in periodic patterns. In the last part of the article, we discuss how the further combination of hierarchical with over-Turing, hierarchical and Turing or Turing and over-Turing networks leads to resulting patterns that can be seen, to a large extent, as the combination of the types of patterns possible from each of the three classes of networks.

Throughout our article, we have considered that cells are arranged in simple regular ways over space. In our arguments, the only element that depends on the geometry of the cellular arrangement is the explicit form of the eigenmodes in (6), whereas the construction of the dispersion relation itself is geometry-independent. The existence of a limited number of topological classes capable of generating patterns does not depend on the geometry of space (i.e., on how cells are distributed within the embryo). The type of pattern transformations for each class does not depend on the geometry of space either (e.g., Turing networks lead to periodic patterns, H networks do not, etc). The geometry of space, however, can change specific aspects of the specific resulting patterns possible from each specific gene network. Thus, for example, the height, spacing and shape of concentration peaks and valleys can change close to the boundaries of the system (Diambra & Costa, 2006; Glimm et al., 2014; Castelino et al., 2020). In Turing networks it is even possible that the resulting patterns changes substantially depending on the shape of the system (e.g., from producing stripes to producing spots), but in all cases the type of pattern, in the coarse way we have defined, remains periodic (Murray, 2002).

All the gene network topologies capable of pattern transformations that we report have been individually reported before in one way or another (Turing, 1952; Salazar-Ciudad et al., 2001; Cotterell et al., 2015; Wang, 2022). What has not been shown before is that these are the only possible gene network topologies in non-trivial pattern transformation, given the biologically reasonable restrictions on f that we describe. We think this is important both for theory and empirical studies. For theory, it is important because we provide a general description of what is possible at the level of non-trivial pattern transformations and underlying gene network topologies (as long as there is no cell movement). Experimentally, this is also important because it provides a simple guideline for experimental developmental biologists trying to understand pattern transformations in specific organs over a set of developmental stages. Thus, if no cell movements occur during these stages, the underlying gene network should have a network topology and variational properties among the ones described in this article (see Fig.13). This significantly narrows down the range of possible underlying gene networks and patterning mechanisms to consider when trying to understand the development of a multicellular system. Similarly, the results of this article allow to improve the study of the evolution of gene networks since they show that the gene networks capable of non-trivial pattern transformations can only change within a finite set of classes (Salazar-Ciudad et al., 2001).

Figure 13:

Of the three classes of gene network topologies we identify, two have been widely studied before: hierarchical and Turing networks. Most of the research in Turing gene networks has been theoretical (Turing, 1952; Maini et al., 2006; Meinhardt, 2008), but there is also direct experimental research on the involvement of these gene networks in the development of many organs (e.g., Glover et al., 2017; Johnson et al., 2023; Tzika et al., 2023; Tseng et al., 2024, just to name some recent ones), although the underlying networks are usually more complex than the ones studied theoretically. Although hierarchical networks are usually not called this way, they constitute the bulk of the gene networks experimentally studied (Salazar-Ciudad et al., 2001; Gilbert and Barresi, 2023) and they are, by far, the easier to understand. This, and the fact that many networks seem hierarchical when only some of their interactions are known (Salazar-Ciudad, 2009), may have biased research to focus on them. Although over-Turing gene networks seem simpler and easier to understand than Turing networks, they have only been discovered very recently and only based on theoretical work (Cotterell et al., 2015; Marcon et al., 2016; Wang et al., 2022). It is thus, an open and suggestive question whether this class of gene networks is widely used in pattern transformations in multicellular systems. In fact, it is perfectly possible that organs that are though to use Turing networks may actually be using over-Turing networks instead, since the latter produce resulting patterns that are similar to the former, especially when the organ is effectively 1D.

Data availability

All data is provided within the article itself (this is a purely theoretical article).

Acknowledgements

We thank H. Cano-Fernández, S. Greendal, C. Hedges and A. Loreto-Velázquez for their comments. We are also grateful to the Spanish Ministry of Science and Innovation for funding PID2021-122930NB-100 and CNS2022-135397 to I.S.-C. (MCIN/AEI/10.13039/501100011033/FEDER, UE) ; and to the Catalan Office of Universities and Research for 2022FI_B 00610 to K.M.-A. This work was likewise supported by the Spanish State Research Agency, through the Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence in R&D (CEX2020-001084-M). We thank CERCA Programme/Generalitat de Catalunya for institutional support.

Additional files

SI

Additional information

Funding

MEC | Spanish National Plan for Scientific and Technical Research and Innovation (Plan Estatal de Investigación Científica y Técnica y de Innovación) (PID2021-122930NB-100)

• Isaac Salazar-Ciudad

MEC | Spanish National Plan for Scientific and Technical Research and Innovation (Plan Estatal de Investigación Científica y Técnica y de Innovación) (CNS2022-135397)

• Isaac Salazar-Ciudad