## Abstract

Boundary formation in the developing neuroepithelium decides on the position and size of compartments in the adult nervous system. In this study, we start from the French Flag model proposed by Lewis Wolpert, in which boundaries are formed through the combination of morphogen diffusion and of thresholds in cell responses. In contemporary terms, a response is characterized by the expression of cell-autonomous transcription factors, very often of the homeoprotein family. Theoretical studies suggest that this sole mechanism results in the formation of boundaries of imprecise shapes and positions. Alan Turing, on the other hand, proposed a model whereby two morphogens that exhibit self-activation and reciprocal inhibition, and are uniformly distributed and diffuse at different rates lead to the formation of territories of unpredictable shapes and positions but with sharp boundaries (the ‘leopard spots’). Here, we have combined the two models and compared the stability of boundaries when the hypothesis of local homeoprotein intercellular diffusion is, or is not, introduced in the equations. We find that the addition of homeoprotein local diffusion leads to a dramatic stabilization of the positioning of the boundary, even when other parameters are significantly modified. This novel Turing/Wolpert combined model has thus important theoretical consequences for our understanding of the role of the intercellular diffusion of homeoproteins in the developmental robustness of and the changes that take place in the course of evolution.

## INTRODUCTION

The specification of territories in the nervous system relies on the precise positioning of boundaries between different functional areas (Flanagan, 2006; Kicheva et al., 2012; Kiecker and Lumsden, 2005). Each territory is characterized by the expression of a specific combination of molecular marks, including transcription factors (TFs), before developing into areas endowed with specific functions (O'Leary et al., 2007; Zilles and Amunts, 2010). The emergence of compartments in the cerebral cortex or in the spinal cord is a paradigmatic example of this process. From a theoretical perspective, the specification of territories in the nervous system is a particular instance of the general phenomenon of patterning. Alan Turing provided the first theoretical model of how patterns form. In his 1952 seminal article ‘The chemical basis of morphogenesis’, Turing explains how reaction-diffusion properties of two morphogens, in the presence of a catalyst, can lead to the emergence of heterogeneities, even if the tissue is initially homogeneous (Turing, 1952; see Fig. 1 for a schematic description of the model). This universal pattern formation mechanism through Turing instabilities has become increasingly popular in the developmental biology community (Kang et al., 2012; Marcon and Sharpe, 2012; Sheth et al., 2012; Raspopovic et al., 2014; Xu et al., 2014). In Turing's model and its enriched versions, in particular those proposed by Meinhardt and Gierer (2000), the interaction of a limited number of molecular species can create regular spatial patterns, provided that they exhibit different diffusion constants and have auto-activating and reciprocal inhibitory properties. In all cases, Turing-like mechanisms alone do not lead to the emergence of predictable shapes.

Another popular patterning mechanism has been proposed by Lewis Wolpert (1969) with the concept of Positional Information (PI). This model, also known as the French Flag model (FFM), requires a continuous morphogen gradient and the existence of thresholds, see Fig. 1. A typical abstract example is the differentiation of cells into blue, white and red populations when exposed to high, intermediate or low morphogen levels (i.e. the FFM), each territory corresponding to the expression of specific genes, in many cases transcription factors (TFs) that define specific areas within the neuroepithelium. This model has since evolved considerably to take into account the complexity of the cellular environment (Hornung et al., 2005; Kerszberg and Wolpert, 2007; Lander, 2007; Xiong et al., 2013).

If one compares the two models, Turing's model allows the formation of precise and neat boundaries but suffers from the absence of a historical pre-patterning that leads to a lack of reproducibility in their positioning. By contrast, the PI model provides a pre-pattern that constrains the positioning, but suffers from fuzziness, owing to an uncertainty in the morphogen concentration at which a threshold appears (especially when the morphogen slope is shallow). This represents a serious difficulty, as discussed by Gregor and colleagues (2007a). In addition to the positioning of boundaries, one has to consider the fate of misplaced cells not expressing a TF combination corresponding to their territory. Because, in PI models, each cell ‘works for itself’, cells close to thresholds may differentiate into different types, leading to a salt-and-pepper pattern in the region of the boundary. In the most parsimonious version of the model (no other mechanism added), the only solutions are migration or death of misplaced cells (Kiecker and Lumsden, 2005; Xiong et al., 2013), requiring additional mechanisms and information to regulate cell migration/guidance and/or death.

It might thus be useful to verify whether recent findings in developmental biology may permit the reconciliation of the advantages of the two models. In vertebrates, the most popular illustration of the PI theory is provided by the compartmentalization of the neural tube in response to the diffusion of the ventral and dorsal morphogens sonic hedgehog (Shh) and bone morphogenetic protein (BMP), respectively (Ribes et al., 2010; Dessaud et al., 2010). A continuous gradient activates ventral and dorsal genes, and territories are formed that express distinct TF subsets (Ashe and Briscoe, 2006; Dessaud et al., 2010, 2007; Kiecker and Lumsden, 2005). In this model, differentiation is based on the almost general rule that within a developing neuroepithelium, each side of a boundary expresses a TF, in most cases a homeoprotein (HP) transcription factor, which amplifies its own expression and represses that of its counterpart (on the other side). This is illustrated, among many other examples, by the Pax6/Nkx2.2 dorsoventral boundary and the Otx2/Gbx2 anteroposterior boundary in the neural tube, or the Emx2/Pax6 boundary in the cortex (Briscoe et al., 2000; Brodski et al., 2003; Joyner et al., 2000; O'Leary et al., 2007). An important novelty of this study is to introduce in the calculations the intercellular transfer of HPs allowed by two short peptidic sequences present in their DNA-binding domain (Spatazza et al., 2013a,b; Joliot and Prochiantz, 2004; Sugiyama et al., 2008; Wizenmann et al., 2009; Kim et al., 2014; Miyata et al., 2012; Yoon et al., 2012).

Direct communication between nearby nuclei in the context of cell assemblies is reminiscent of the studies where direct morphogenetic functions were attributed to transcription factors diffusing in the fly embryo at the syncytial stage (Driever and Nüsslein-Volhard, 1988a,b). The parallel is made even more striking by recent studies suggesting that such local diffusion between nearby nuclei represses developmental noise, allowing the precise positioning of transcriptional domains (Gregor et al., 2007a,b). It is not usual to think of a transcription factor as a morphogen, and if Bicoid was easily labeled ‘morphogen’ in spite of being a HP transcription factor, it is rather because of its graded expression and of the fact that the *Drosophila* embryo is a syncytium that allows Bicoid direct transfer from nucleus to nucleus. Therefore, the similarity between the Bicoid model and our own hypothesis is limited to the fact that HP diffusion is involved. Indeed, Bicoid in the fly is a morphogen, as defined by Wolpert, whereas, in our model, HPs are morphogens by Turing's definition.

Indeed, their intercellular transfer added to self-amplification and reciprocal inhibition properties may convey to HPs the quality of local Turing's morphogens. In that sense, nature may have combined Turing's morphogen diffusion (HPs) with PI provided by classical morphogen gradients (e.g. Shh). This reasoning is at the basis of the parsimonious model presented in this study that takes into account the presence of morphogen gradients, as in the PI theory, but also incorporates a Turing-like mechanism based on the local diffusion of HP transcription factors (see Fig. 1). A major and counterintuitive finding of our study is that, even in the limit of infinitesimal diffusion, HP transfer across cells is sufficient to ensure precise boundaries with reliable location. Beyond the case explored here in the context of neural development, this study has led to us to discover an important mathematical property, universal in systems with competing species subject to diffusion, as shown in another study (Perthame et al., 2014). This paper does not present these formal mathematical details, but illustrates this theory with one minimalistic example that can be precisely analyzed mathematically and simulated. In conclusion, it is demonstrated that the addition of the simple property of HP transfer integrates a local Turing's mechanism within the PI model first proposed by Wolpert and provides a very parsimonious model for the formation of precise and stable boundaries.

## RESULTS

We propose a model that takes into account the basic mechanisms at play during neuroepithelium development when different combinations of genes are expressed in abutting differentiating domains, including HPs, that dictate the morphological and functional fate of territories (Kiecker and Lumsden, 2005).

The simple and parsimonious model that we propose is schematically represented in Fig. 1. It considers that the differentiation between two areas, *A* and *B*, is driven by the dynamic competition between the expression of two homeogenes associated with distinct HPs: *T _{A}* and

*T*. Three important processes propel this mechanism:

_{B}The presence of one or several morphogens that form gradients along the developmental axis.

The competition between the different HPs through autocatalytic activation and reciprocal inhibition.

The activity of non cell-autonomous HPs captured from the closest neighboring cells (up to three cell ranks) through extracellular diffusion.

As the neuroepithelium develops, epigenetic phenomena take place and modify the homeogene expression repertoire by favoring those that are the expressed the most. Eventually, a classical self-limiting process, such as saturation within the cell, imposes a plateau to gene expression.

All these phenomena provide a well-defined equation for the evolution in time of the HPs in each cell. We provide the detailed mathematical model in the Materials and Methods. Overall, the model qualitatively depends on only three effective parameters that are the ratio of (1) the autocatalytic activation rates, (2) the saturation/inhibition rates and (3) the extracellular diffusion rates.

The problem of boundary formation and stability consists in determining: (1) whether the piece of neural tissue clearly splits into separate regions in which cells either express *T _{A}* or

*T*; and (2) the site where this partition takes place as a function of initial conditions and the stability of the boundary position upon random parameter variations.

_{B}### Ambiguous boundary in the absence of non cell-autonomous processes

In the absence HP diffusion, the behavior of each cell is governed by an autonomous equation (independent of the behavior of the other cells) that depends on the local concentration of morphogen. Within each cell, homeogenes compete for expression, and the outcome of this process is that the ‘winner-takes-all’: one TF will be expressed at the expense of the other that eventually disappears. The differentiation of a cell into A or B depends on their position within the morphogen gradient. We demonstrate that in the regions where the expression of one HP (e.g. *T _{A}*) is highly promoted by the morphogen gradients, the cells can only differentiate into type A: morphogens ‘select’, in these regions, the winner (see Model in the supplementary material). However, in the regions of intermediate concentrations of morphogen, the cells can differentiate into A or B, and the fate of one cell is governed by initial concentrations of HPs and the transcriptional noise. In other words, there exists a non-trivial set of morphogen concentration levels in which the system has a stochastic patterning. In a differentiating tissue, the region corresponding to these morphogen concentrations is ambiguous: the system displays an exponential number of possible stable differentiated states. (If there are

*k*cells in the ambiguous region, each cell can be of type

*A*or

*B*independently, therefore the total number of possible solutions is equal to 2

*.)*

^{k}The sensitivity of the differentiation process within the ambiguous region leads, in physiological noisy conditions, to an unpredictable patterning, and the vast majority of the solutions display an alternation between the two cell types, precluding the definition of smooth boundaries between cells but rather leading to a salt-and-pepper pattern. In the absence of additional processes leading to cell reprogrammation, migration or death (Kiecker and Lumsden, 2005), this salt-and-pepper regime is ubiquitous (see Model in the supplementary material). This is a property of a wide class of abstract models of cell differentiation, where systems of competing species yield two winner-takes-all states that, except in extremely fine-tuned situations, do not change stability exactly at the same points in space, and hence are generally both stable in a region of space defined as ambiguous.

This is clearly visible in the numerical simulations of a two-dimensional tissue in Fig. 2: in the two cases, the initial randomness persists in the final state, and the overall shape of the domain dramatically depends on the choice of the initial condition. In all cases, we indeed observe the salt-and-pepper type of boundary due to the randomness in the initial condition. In order to illustrate this sensitivity to noise and to initial concentrations, we present three examples with no specific initial pattern, a partially differentiated initial pattern or a fully differentiated initial pattern. The end state of the differentiation process always show an imprecise boundary, and we can see the dramatic dependence on the initial condition with, in these three cases, a patterning that is globally very distinct.

It can be concluded that in the absence of HP diffusion, different steady-state solutions appear and remain stable. Moreover, the differentiated domains are highly irregular and subject to fluctuations upon variation of the initial conditions and parameters.

### Unpredictable patterns in the absence of morphogen gradients

We now discuss the behavior of a differentiating tissue within which molecular species diffuse but in the absence of positional information given by morphogen gradients. Turing was the first to suggest that the diffusion of self-activating and reciprocal inhibitor elements is at the basis of boundary formation (Turing, 1952). In order to support pattern formation, the original Turing model makes the assumption that an additional molecular species plays a catalytic role on the expression of both of A and B. This molecular species contrasts with the graded expression of the different morphogens of the PI model in at least two aspects: it has a no spatial source, and therefore does not define any preferred place in space for one species to be expressed; and it promotes the expression of both A and B.

A major mathematical finding in this model is the now-called Turing's instability: when the rates of diffusion of the two species are very different, several homogenous ‘winner-takes-all’ abutting territories emerge at random places (the ‘leopard spots’, see Fig. 1). The patterns so generated are unpredictable: they are highly sensitive to noise and to initial conditions.

In our model, one can consider HPs as Turing's self-activating and reciprocal inhibitor species, and the morphogen showing a graded expression along the differentiating pluricellular tissue (central in Wolpert's FFM) plays the role of Turing's catalytic species. But it no longer has a spatially homogeneous concentration. Its graded monotonic expression will stabilize the Turing patterns, leading to regular, predictable and highly reproducible boundaries between distinct ‘winner-takes-all’ abutting territories, as we now show. Moreover, in contrast to Turing's instability, our model does not necessitate that the competing species have different diffusion rates: precise and smooth patterning arises as soon as both species have the ability to diffuse.

### Precise patterning for competitive systems with spatial cues and HP diffusion

From the two above sections, we conclude that HP diffusion in the absence of morphogen gradients (Turing) leads to unpredictable patterning with clearly defined boundaries, while the presence of spatial cues (positional information) in the absence of HP diffusion (Wolpert) yields to a patterning predictable ‘at large’ but with imprecise boundaries. Our model combines both spatial cues (external morphogen gradient) and HP diffusion across cell membranes. The classical morphogen in the Wolpert's definition (e.g. Shh) creates zones of expression of distinct HPs (the ‘French Flag’) with blurred and unstable boundaries, and HPs are now locally diffusing secondary morphogens in the Turing's definition (self-activation and reciprocal inhibition). These two processes, when combined, lead to smooth and predictable pattern formation and the location of the boundary is very robust to random fluctuations of the initial conditions and parameters, even at very low diffusion levels. This is a surprising property of the equations. Indeed, this stabilization takes place for arbitrarily small values of the diffusion constants, meaning that most solutions present in the case σ* _{A}*= σ

*=0 disappear in favor of a unique solution with precise front location. This stabilization property is mathematically demonstrated in our study on general models of competitive-type systems driven by monotonic gradients (Perthame et al., 2014).*

_{B}In order to illustrate this phenomenon, simulations of the system are provided in Fig. 3. For the sake of consistency with the biological problem, we performed the simulations adopting the topology of a neural tube. Two sources (representing Shh and BMP, for example) are fixed at the floor plate and roof plate, respectively, and free diffusion was simulated to form the gradients. The BMP source was arbitrarily chosen as being stronger than the Shh source (ratio 3:2), and initial HP concentrations were chosen close to zero, with small fluctuations across different cells. In the absence of diffusion emerges a noisy boundary that is consistent with the previous analysis. But even a very small diffusion leads to a dramatic stabilization and regularization of the boundary, at a location that depends only on the parameters of the system (strength of the gradients and intensities of the reactions) but not on the choice of the initial conditions (see Model in the supplementary material).

This dramatic regularization and stabilization of the boundary position is a direct consequence of HP local diffusion (see Fig. 3, right). First, in contrast with the cell-autonomous situation, diffusion prevents the persistence of small isolates of one cell type, e.g. *B*, within a large domain of the other cell type, e.g. *A*. If such an isolate were to appear, diffusions of *T _{B}* and

*T*(out and into the isolate, respectively) would rapidly translate into a ‘

_{A}*T*-takes-all’ situation. In addition to forcing isolated cells to adopt the identity of their dominant neighbors, HP diffusion also contributes to the determination of a highly conserved boundary position between territories

_{A}*A*and

*B*, even for a large range of initial conditions. Indeed, as in the cell-autonomous situation, the regions of high morphogen concentration rapidly differentiate into the

*A*or

*B*type, thus anchoring the differentiation of the field at both of its extremities. Closer to the future boundary, HP diffusion extends the competing domains until the two fronts meet, resulting in continuous and monotonous

*T*and

_{A}*T*gradients. Then local competition based on HP local diffusion and the ability of the two HPs to self amplify and to repress each other, will settle a smooth boundary along the level sets of the morphogen gradients.

_{B}### Stability of the front

In physiological conditions, several phenomena may occur and perturb the position of the front. An important source of variability comes from the heterogeneity of the cell population, and in particular from the fact that the characteristics of gene expression vary from cell to cell. Moreover, noise can arise from cell division, cell death and random movements of the cells that modify the sensed value of the morphogens, which may join their effects to perturb the position of the boundary. In fact, the boundary location predicted in the idealized model proves surprisingly resilient in all these situations, as we now illustrate.

In order to quantify the sensitivity of the boundary location to the heterogeneity of the cell population, we investigated the effect of having heterogeneous rates of self-activation and inhibition between TFs (i.e. varying from cell to cell). These two parameters completely characterize gene expression in a given cell in our model. We considered, for example, these rates randomly chosen according to a normal distribution, with mean *g=*1 and different standard deviations λ (see Fig. 4). The end-state for λ=0.05 is displayed in Fig. 4 (left) superimposed with the end-state in the homogenous case λ=0 (dashed line). We can observe that even if the precise concentration levels in the different cells are modified compared with the homogeneous predicted solution, the position of the front barely changes. This is due to the very sharp drop of concentration across the boundary. We quantified this stability by looking at the distribution of the front location for 500 independent realizations. The histograms of the front location are displayed in Fig. 3 (right) for different values of the heterogeneity level. Even for large values of the heterogeneity, the position of the front is conserved relatively precisely. For example, for a noise on the coefficients of λ=0.05, the front position is barely modified (maximal displacements of two cells); for λ=0.2, although the solution appears relatively different from realization to realization, the front location remains relatively stable, with maximal errors of 10 cell ranks (on a total of 100 cells)⇓.

Cell division occurring during development may also result in variations in the position of the boundary. In order to investigate this effect, we simulate the system with a variable *N* that randomly depends on time. *N* is set to 100 at initial time, and we consider that one new cell appears as a Poisson process (i.e. cell division occurs at independent exponentially distributed times). When a cell divides, it shares its contents (number of TFs *T _{A}* and

*T*) between the two new cells that conserve the same epigenetic marks as the mother cell (here, transcription intensities

_{B}*D*and

_{A}*D*). A typical trajectory of the front is depicted in supplementary material Fig. S4. Numerical results show that the position of the boundary is barely modified by this process: transient displacements of the boundary that may arise when divisions occur close to the boundary are rapidly overcome, as visible in supplementary material Movie 1. The stability of the boundary location upon variations of morphogen activity (parameters

_{B}*F*and

_{A}*F*) was analyzed in order to account for possible random movements of the cells and fluctuations of the environment [e.g. random arrival of morphogen molecules at their target and readout noise, see Gregor (2007a)]. Again, the front remains stable with time, varying by, at most, a few cell ranks even for large values of the noise (see supplementary material Fig. S4).

_{B}## DISCUSSION

In this paper, we describe a parsimonious model for the formation of boundaries within an epithelium. It is in the spirit of the seminal paper where Lewis Wolpert proposed, almost 50 years ago, the FFM to explain boundary formation and, in many ways, it extends this model (Wolpert, 1969). We started with the idea that the compartments created by a diffusing morphogen, as in the FFM, are marked by the expression of secondary morphogens (not morphogens in the presently most accepted term, but in the sense coined by Turing) of the HP transcription factor family and introduced two hypotheses: first, that HPs diffuse locally between cells; second, that HPs on either side of a boundary activate themselves and are reciprocal inhibitors (at the transcription level). HP diffusion was indeed demonstrated in a number of biological systems and situations (Brunet et al., 2007, 2005; Di Lullo et al., 2011; Spatazza et al., 2013a,b; Stettler et al., 2012; Sugiyama et al., 2008; Wizenmann et al., 2009; Kim et al., 2014; Miyata et al., 2012; Yoon et al., 2012). In addition, the sequences responsible for HP secretion and internalization are highly conserved between HPs (Joliot and Prochiantz, 2004), supporting the idea that most HPs are local ‘Turing’ morphogens. The second hypothesis is also supported by a large number of experiments and is illustrated by the fact that genetic gain or loss of function of one of the two ‘abutting’ HPs results in a shift in boundary position (Millet et al., 1999; O'Leary et al., 2007; O'Leary and Sahara, 2008; Puelles et al., 2004; Toresson et al., 2000; Yun et al., 2001; Dessaud et al., 2008). From a mathematical standpoint, the phenomenon of disambiguation and stability of the boundary is relatively surprising, because for arbitrarily small values of the diffusion constants, most solutions present in the case *σ _{A}*=

*σ*=0 (no HP diffusion) disappear in favor of a unique solution with precise boundary location. The characterization of similar phenomena in partial differential equations (PDE) in the small diffusion limit is a very interesting mathematical problem and constitutes an active field of research (Bages et al., 2012). It is actually possible to prove that, in the continuous limit, the viscosity solutions of this equation (i.e. asymptotic solutions in the limit where the diffusion tends towards zero) present a unique and perfectly defined boundary.

_{B}Our model requires only three molecules to form a boundary (one graded morphogen and two HPs). It is, thus, as parsimonious as the FFM, while avoiding the introduction of explicit thresholds. Its main advantage is that the robustness of the positioning of boundaries is highly increased by the diffusion and reciprocal inhibition HP properties. Our model can also be compared with that proposed by Turing (Turing, 1952). Indeed, HPs can be considered as morphogens in the sense of Turing because they amplify their own expression, are reciprocal inhibitors and have non-local properties. However, in the reaction-diffusion Turing's model, boundaries appear in a morphogenetic field due to dynamic instabilities arising when the rate of diffusion of the two species in competition are sufficiently different. The mechanisms by which a biological system could be composed of species with very different diffusion constants are still largely unknown. Moreover, when Turing instability forms a pattern, the boundary location is unpredictable. In sharp contrast, our model forms regular and predictable patterns regardless of the respective value of the diffusion of the two species. In other words, the diffusivity of the species in competition do not need to be different to form a boundary; moreover, the boundary forms at a precise position and remains stable under variations of the initial conditions and fluctuations of parameters.

By putting aside the ability of HPs themselves to form a gradient through their iterative induction across a large territory that was considered recently (Holcman et al., 2007; Kasatkin et al., 2008), we have been able to base our developments only on solidly established data and to neglect several parameters, thus giving direct access to the comparison with the models proposed by Turing and Wolpert. If we think of other models, such as those proposed by Hans Meinhardt, by concentrating on HP local diffusion we could also make the economy of the long-range inhibition hypothesis (Meinhardt, 1978, 1983, 1994; Meinhardt and Gierer, 2000). Indeed, our model does not preclude that such long-range inhibitions take place, but does not need it in the first place. Other studies have proposed that bistable dynamics could be the source of reliable patternings (Lewis et al., 1977). Their model is somewhat simpler in that it only considers auto-activation (and ignores cross-inhibition) and the presence of a long-range gradient. But the cells no more respond monotonically to gradients: they have a more complex nonlinear dynamics that, over a certain range of values of the morphogen gradient, can differentiate into different populations. This bistability is naturally built into our model and emerges from the competition between the two species. Yet, in the absence of diffusion, any bistable system bears ambiguity on the patterning: the boundary will, in particular, depend on the initial condition (as in our system in the absence of diffusion). However, similar to what we show here, adding a diffusion term in bistable models such as proposed by Lewis et al. (1977) would allow stabilization of the boundary.

This is actually a deep mathematical property. From a mathematical viewpoint, the problem of neurodevelopment in the presence of diffusing HPs is one of the few examples in which biology led to the discovery of a universal mathematical property. Motivated initially by the mechanism of gene expression described here, we demonstrate that all competitive systems in the presence of monotonic cues yield the formation of a stable and regular boundary between two abutting domains; this property is valid even at arbitrarily low levels of diffusion (Perthame et al., 2014). This mathematical result, beyond applications to other domains, has major implications from a biological viewpoint. Indeed, it ensures that the phenomenon of reliable pattern formation does not depend on the details of the model under consideration, but only on a few qualitative properties that are very natural in the context of neurodevelopment.

Because HPs are very ancient molecules present in all phyla (Derelle et al., 2007) and as transduction takes place in plants and animals, it is speculated that this mode of signaling was operating in the first multicellular organisms. In that sense, it may have preceded other signaling mechanisms based on classical signaling entities (e.g. growth factors and their receptors) and pathways. Indeed, in the case of Bicoid (Dubnau and Struhl, 1996; Mayfield, 1996; Rivera-Pomar et al., 1996), it was shown that internalized HPs could regulate local translation (Alvarez-Fischer et al., 2011; Stettler et al., 2012; Yoon et al., 2012). The recruitment, later in the evolution of multicellular organisms, of classical signaling pathway is likely to have added robustness to the formation of territories and to other functions involving HP transduction. For example, it was shown that the patterning of terminals from the retinal ganglion cells within the tectum/superior colliculus depends on an interaction between Engrailed HP and Ephrin/Eph signaling (Stettler et al., 2012; Wizenmann et al., 2009). How HP and classical signaling pathways have evolved in parallel and by interaction is of the highest importance to understand the morphogenesis of multicellular organisms and its evolution. In this context, proposing a parsimonious mechanism is a first step in the further analysis of these complex phenomena.

## MATERIALS AND METHODS

We now present the mathematical details of the model we have been developing and using in our analyses and simulations.

### Theoretical description

The model describes the time evolution of the quantities *T _{A}* and

*T*in a spatially extended neural tissue composed of

_{B}*N*differentiating cells. Their dynamics are the result of cell-autonomous mechanisms and non cell-autonomous diffusion. Specifically, they satisfy the equations:

### Cell-autonomous HP competition

The expression of the genes is the result of the competition between the expression of the two combinations of genes modulated in our system by gene expression capacities *D _{A}* and

*D*that evolve according to epigenetic mechanisms that we discuss below. We take into account the following phenomena (described for one combination of gene A; the same phenomena is considered for B):

_{B}Morphogens stimulates TF expression:
The quantity *F _{A}*(

*x*) denotes the rate of production of

*T*induced by the morphogen on cells at location

_{A}*x*. It is a monotonic function along a developmental axis (gradient direction of the morphogen).

The auto-inducer properties of TFs are taken into account by considering that *T _{A}* stimulates its own expression with a positive rate

*g*. This intensity is modulated also by the gene expression capacity

_{A}*D*

_{A}: The cross-inhibition properties imply that the presence of

*T*inhibits the expression of

_{B}*T*causing, in the cell, a decrease of the production rate of

_{A}*T*at a certain rate

_{A}*S*. The simplest way to express this competition is to write: Finally, saturation of the number of proteins inside the cell is taken into account by considering that the rate of production of the species decreases when

_{A}*T*exceeds a certain level. We choose here the logistic saturation law that is classical to ecologists: These equations characterize the expression dynamics within a cell. All phenomena requiring gene expression occur at a rate that is scaled by a coefficient

_{A}*D*, taking into account the epigenetic phenomena. This coefficient accounts for the fact that the more one combination of genes is expressed, the more likely it is to be expressed. This facilitation-inhibition of the transcriptional activity results from the fact that

_{A}*D*is an increasing function of

_{A}*T*and decreasing function of

_{A}*T*: where the map

_{B}*G*is such that, by convention: In this scaling,

*D*=1 corresponds to a maximal expression activity and

_{A}*D*=0 to no gene expression at all.

_{A}### Non-cell-autonomous transfers

In addition to the cell-autonomous mechanisms, and given that homeoproteins are endowed with direct non-cell-autonomous properties, we include in the set of equations what we call a diffusion operator Δ. From a modeling viewpoint, we incorporate in the dynamics of *T _{A}* and

*T*the ability to be transferred to neighboring cells. To emphasize this very local mechanism (Layalle et al., 2011), we limit this diffusion to one cell in all directions. In detail, the time evolution of the transcription factor level

_{B}*T*(

_{A}*x*) within the cell at location

*x*is added a nonlocal term corresponding to the exchange of transcription factors from and towards the set of neighboring cells

*υ*(

*x*) (the number of neighbors is denoted |

*υ*(

*x*)|):

In other words, TFs have the ability to be transferred to all neighboring cells at a rate σ (the intensity in time of the transfer), creating outward inward fluxes.

## Acknowledgements

The authors warmly acknowledge Benoît Perthame for his help in the mathematical analysis of the model and insightful discussions on the mathematical properties of competitive systems that lead to the writing of Perthame et al. (2014).

## Footnotes

**Competing interests**The authors declare no competing or financial interests.

**Authors contributions**C.Q. developed the model, performed the numerical simulations and mathematical analysis, analyzed the data, and wrote the paper; A.P. conceived the study, developed the model, analyzed the data and wrote the paper; J.T. conceived the study, developed the model, performed the mathematical analysis and numerical simulations, and wrote the paper.

**Funding**This study was funded by a grant from the Agence Nationale pour la Recherche [ANR-11-BLAN-069467], by the Global Research Laboratory Program [2009-00424] from the Korean Ministry of Educations, Science and Technology, and by an ERC Advanced Grant HOMEOSIGN [339379].

**Supplementary material**Supplementary material available online at http://dev.biologists.org/lookup/suppl/doi:10.1242/dev.113688/-/DC1

- Received June 23, 2014.
- Accepted March 23, 2015.

- © 2015. Published by The Company of Biologists Ltd