Deciphering and modelling the TGF-β signalling interplays specifying the dorsal-ventral axis of the sea urchin embryo

During embryonic development, cell fate is specified by transcription factors activated in response to instructive signals. Regulatory interactions between signalling molecules and their target genes form networks, called gene regulatory networks (GRN) (Arnone and Davidson, 1997). Deciphering such GRNs is a key for developmental biologists to understand how information encoded in the genome is translated into cell fates, then into tissues and organs, and how morphological form and body plan can emerge from the linear DNA sequence of the chromosomes (Levine and Davidson, 2005). Noteworthy, the GRN orchestrating the morphogenesis of the ectoderm along the dorsal-ventral (D/V) axis of the embryo of the model sea urchin Paracentrotus lividus has started to be uncovered in great detail (Haillot et al., 2015; Lapraz et al., 2015; Li et al., 2012, 2013, 2014; Molina et al., 2018; Range et al., 2007; Saudemont et al., 2010; Su et al., 2009 and reviewed by Molina et al., 2013).

The ectoderm of the sea urchin larva comprises two opposite ventral and dorsal territories, separated by a central ciliary band (Fig. 1A): the ventral ectoderm is the territory at the centre of which the mouth will be formed. Specification of the ventral ectoderm relies on signalling by Nodal, a secreted growth factor of the TGF-β family (Duboc et al., 2004; reviewed by Molina et al., 2013). The expression of nodal is turned on by maternal factors, whereas Nodal stimulates and maintains its own expression through a positive feedback circuit (Bolouri and Davidson, 2010; Range et al., 2007; Range and Lepage, 2011). Nodal is zygotically expressed and is thought to dimerise with a maternally expressed TGF-β ligand called Univin (Range et al., 2007); the Nodal-Univin heterodimer promotes Alk4/5/7 signalling and the activation of Smad2/3 together with Smad4. The ventral ectoderm boundary is thought to be positioned by the activity of the product of the Nodal target gene lefty, which prevents the expansion of nodal expression beyond the ventral ectoderm region via a diffusion-repression mechanism (Chen et al., 2004; Cheng et al., 2004; Duboc et al., 2008; Molina et al., 2013; Sakuma et al., 2002).

The ciliary band ectoderm is a proneural territory located between the ventral and dorsal ectoderm (Angerer et al., 2011). The ciliary band is made of prototypical cuboidal epithelial cells and runs along the arms of the pluteus larva. Unlike the specification of the ventral and the dorsal ectoderm, which actively requires TGF-β signalling, the specification of the ciliary band tissue does not rely on Nodal or BMP signalling, and this tissue develops as a ‘default’ state of the ectoderm in the absence of these signals (Saudemont et al., 2010).

The dorsal ectoderm is the territory that will differentiate into the apex of the pluteus larva. Its specification relies on the diffusion of ventrally synthesised BMP2/4, which promotes dorsal fates by activating phosphorylation of Smad1/5/8 via the activation of the BMP type I receptors Alk1/2 and Alk3/6. The inhibition of BMP signalling on the ventral side and the translocation of BMP2/4 to the dorsal side require the product of the chordin gene, which is activated in the ventral ectoderm downstream of Nodal signalling (Lapraz et al., 2009). Glypican 5 is expressed downstream of BMP2/4 signalling and contributes to stabilise BMP signalling on the dorsal side via a positive feedback circuit (Lapraz et al., 2009). In addition, the BMP ligands Admp1 and Admp2 dampen BMP signalling fluctuations through an expansion-repression mechanism (Ben-Zvi et al., 2008; De Robertis, 2009; Joubin and Stern, 2001; Kimelman and Pyati, 2005; Lapraz et al., 2015; Lele et al., 2001; Reversade and De Robertis, 2005; Willot et al., 2002). This mechanism, which relies on the transcriptional repression of admp1 expression by BMP2/4 signalling, allows admp1 expression to increase and ADMP1 protein to be shuttled to the dorsal side by Chordin when BMP signalling decreases. Thus, an increase in admp1 expression compensates for the reduction of the intensity of BMP signalling.

One prominent feature of the D/V onset specification is that it relies extensively on the maternal inputs Panda and Univin, which repress and promote, respectively, ventral fate (Range et al., 2007; Haillot et al., 2015; reviewed by Molina and Lepage, 2020). Univin is a TGF-β related to Vg1 and GDF1/3. Univin is essential for Nodal signalling through the Nodal/Activin receptors; in the absence of Univin, Nodal expression stops and D/V axis specification fails. Panda is a member of the TGF-β superfamily presumed to repress ventral fate by a still unidentified mechanism (Haillot et al., 2015). Finally, the transcriptional repressor Yan/Tel acts as a negative regulator of nodal expression, the expression of which is required downstream of Panda (Molina et al., 2018).

Previous studies have shown that Nodal produced by the ventral ectoderm is a strong ventralising signal (Duboc et al., 2004; Lapraz et al., 2015; Saudemont et al., 2010). Overexpression of nodal causes all ectodermal cells to adopt a ventral fate (Fig. 1B), whereas a loss of Nodal function prevents specification of both ventral and dorsal fates and causes the ectoderm to differentiate as a ciliary band (Fig. 1C). Conversely, the activity of BMP2/4 protein promotes dorsalisation, and overexpression of BMP2/4 forces all ectoderm cells to adopt a dorsal fate (Haillot et al., 2015; Lapraz et al., 2009, 2015) (Fig. 1D). In contrast, removing the function of the BMP2/4 ligand from fertilisation onwards prevents specification of dorsal fates, leading to the formation of an ectopic ciliary band territory in the dorsal region (Fig. 1E). In addition, knocking down panda expression on top of BMP2/4 loss-of-function enables nodal to be expressed through the dorsal side of the ectoderm and promotes ventral fates in all ectodermal cells (Fig. 1F). Conversely, a local panda overexpression promotes dorsal fates, suggesting that panda is sufficient to break the radial symmetry of the embryo and necessary to specify the D/V axis (Haillot et al., 2015). The BMP and Nodal ligands thus show strongly antagonistic activities. However, the mechanism underlying this antagonism and the resulting cell fate decision still awaits clarification.

Owing to the largely non-cell autonomous nature of the D/V GRN, to the many events of protein diffusion and to the intertwined feedback circuits involved, an intuitive understanding of the logic of the network is hard to obtain. For example, Nodal and BMP2/4 are co-expressed in the ventral territory, but the corresponding signalling pathways are operating at opposite poles of the D/V axis. In this context, a model of the D/V GRN is very useful to formalise the complex regulatory interactions at stake (Wilczynski and Furlong, 2010). A GRN can be modelled as a static regulatory graph using standardised graphical conventions to represent the molecular interactions between relevant regulator components (Arnone and Davidson, 1997). Such regulatory graphs can be supplemented with threshold levels and Boolean rules to obtain a predictive dynamical model (Mbodj et al., 2016; Peter et al., 2012; Thomas and D'Ari, 1990).

In the present study, we aim at integrating the existing knowledge on D/V axis specification mentioned above into a common predictive framework. In this respect, we built a logical model (i.e. using Boolean algebra for the regulatory rules) of the sea urchin D/V specification GRN to: (1) assess its accuracy, (2) compare the simulations of different perturbations with the observed gene expression patterns, (3) explore the dynamics of the system and (4) develop a multicellular framework to test the ability of the model to generate the spatial patterns observed in wild-type and perturbed embryos. Although the model does not address the totality of the underlying mechanisms mentioned above (e.g. the ventral fate repression driven by Panda or the translocation of BMP2/4), we gain strong insights into several key regulations, such as the antagonism between Nodal and BMP2/4.

In order to build a model of the GRN driving the D/V patterning of the sea urchin ectoderm, we started by compiling experimental data to identify the key genes and regulatory interactions (Fig. 2). The raw data that provided the spatial and temporal expression information to build the model were derived from the analysis of high resolution in situ hybridisations, northern blot and systematic perturbations experiments (Table S1). Loss-of-function experiments via morpholino injections are particularly important for GRN reconstruction as they allow us to test whether a gene is required for the activation of another gene. Gain-of-function experiments via mRNA injection were also used in many instances to test for the ability of a given gene to induce another gene when overexpressed.

Based on these data, we first built a regulatory graph representing the D/V GRN in a single cell, using the software GINsim (https://ginsim.org; Naldi et al., 2018a). The model has been built through an iterative process, alternating simulations and model refinements, to optimise our simulation results compared with experimental observations (Fig. 2). An extensive analysis of published data enabled us to delineate relevant regulatory interactions and to formalise them in terms of a signed, directed graph, complemented with Boolean rules (see Materials and Methods). The resulting model was then subjected to computational analyses, which were confronted to experimental data. The analysis of discrepancies led us to revise the network and the Boolean rules (e.g. by incorporating a new interaction), until the simulations qualitatively recapitulated all relevant experimental observations (Fig. 2). The main challenge encountered during this process consisted in the integration of multiple interactions targeting the same component. Indeed, the number of possible Boolean rules exponentially increases with the number of incoming interactions. These rules must correctly reflect the relative importance of each incoming signal, taking into account different qualitative ranges. These rules are encoded using the classical Boolean operators (AND/OR/NOT), using multilevel variables when justified (see Materials and Methods). At the end of this iterative process, the model was further refined based on the results of novel experiments designed to disentangle the interplay between Nodal and BMP2/4 pathways (Fig. 2).

Regarding the computational analyses performed, we combined three main approaches: (1) a direct computation of stable states for relevant input combinations, (2) stochastic simulations to estimate the probabilities of alternative cell fates (chordin morpholino), and (3) logical multicellular simulations to cover intercellular signalling events (see below: Multicellular simulations emphasise the crucial role of long-range signal diffusion).

We present hereafter the final model and selected simulations recapitulating the patterns observed experimentally. The model encompasses a total of 30 nodes, linked by 23 activations and 15 inhibitions (Fig. 3). The nodes included into the model correspond to signalling and regulatory components (representing proteins or genes), whereas the signed arcs denote regulatory interactions between these components (representing either protein-protein or protein-DNA interactions). Signalling proteins are modelled as input nodes (in yellow in Fig. 3), providing activatory or inhibitory signals through the corresponding membrane receptors.

Each (non-input) component represents a gene (e.g. nodal) or a protein (e.g. Alk receptors), and is classified as ventral (eleven nodes shown in blue in Fig. 3), ciliary band (two nodes shown in pink in Fig. 3) or dorsal (eight nodes shown in green in Fig. 3), according to the reported pattern of expression in the presumptive ectodermal regions. For example, goosecoid is activated by the Nodal cascade in the ventral region in wild-type condition, and is thus considered a ventral gene. The model encompasses the main regulatory components of TGF-β signalling pathways, including the ligands, negative regulators such as the proteins that trigger receptor degradation, downstream transcription factors, and antagonists. Each component of the model is annotated with textual explanations and database links (in particular to PubMed) documenting our modelling assumptions (main references include the GRN diagrams published for P. lividus in Haillot et al., 2015; Lapraz et al., 2009; Lapraz et al., 2015; Molina et al., 2013; Range et al., 2007; Saudemont et al., 2010).

Among the 30 components of the model, 20 are associated with Boolean variables (ellipsoid nodes in Fig. 3, taking the values 0 or 1 depending on their activation state), whereas the remaining components are associated with multilevel variables (rectangular nodes in Fig. 3, associated with three or four integer levels, from zero to 2 or 3, see below). The nine input nodes (shown in yellow in Fig. 3) define 576 possible configurations of input values. Using the Java library bioLQM (Naldi, 2018), we identified 654 stable states, which can be split into three main patterns based on the active nodes. More precisely, 288 of these stable states correspond to a ventral fate, 126 to a ciliary fate and 240 to a dorsal fate (see the Jupyter notebook, https://github.com/ComputationalSystemsBiology/echinodal_notebook).

A key feature of our model consists in the antagonistic activities of BMP2/4 and Nodal. In this respect, additional experiments were conducted to probe the underlying mechanisms. We first tested whether difference in dosage between the corresponding signalling molecules could favour the establishment of a specific cell fate. Treatment with an intermediate dose of BMP2/4 protein resulted in embryos developing with a straight archenteron, no mouth and covered with a ciliary band-like ectoderm (Fig. 4A), which is a prototypical nodal loss-of-function phenotype. In contrast, for intermediate doses of Nodal, the embryos developed with a reduced apex, reminiscent of the BMP2/4 loss-of-function phenotype. These observations suggest that ectodermal cells receive both antagonistic ventralising Nodal and dorsalising BMP2/4 signals, and integrate them at the level of the cis-regulatory sequences of their target genes.

However, as these treatments were performed soon after fertilisation, it was not clear whether the outcome was reflecting an antagonism between Nodal and BMP occurring during cell fates allocation, or whether it resulted from one pathway being activated early and dominantly in all cells of the embryo following the injection of mRNA into the egg. To address this issue, we repeated the Nodal and BMP2/4 treatment at late blastula and at early mesenchyme blastula stage. Treatments with Nodal or with BMP2/4 proteins at late blastula radialised the embryos by respectively inducing ventral or dorsal fate in all ectodermal cells (Fig. 4B). These results confirm that Nodal and BMP2/4 pathways also act antagonistically during fate specification and that a dosage competition plays an important role in cell fate specification.

These results led us to pay particular attention to the encoding of this antagonism in our model. First, to account for dosage effects, we associated ten nodes of the model with multilevel variables (rectangular nodes in Fig. 3, taking values from 0 to 2 or 3). These multivalued nodes allow for a more fine-grained encoding of the activation states of key morphogens and downstream components, the effects of which are dose-sensitive (nodes Nodal, Chordin, BMP2/4, Alk4/5/7, Alk2/3/6, Smad2/3/4, Smad1/4/5/8 in Fig. 3). For example, Nodal activity can either overcome (level 2) or be counteracted by (level 1) the inhibition of Lefty. In addition, Nodal input node is represented by a quaternary variable, in order to distinguish wild-type activity levels (levels 1-2) versus still higher gain-of-function (level 3) activity. Second, the antagonism between the Nodal and BMP pathways is encoded as a reciprocal inhibition between Smad2/3/4 and Smad1/4/5/8 (Fig. 3), which implements the competition of these signalling complexes for the shared molecule Smad4. Noteworthy, each of these two inhibitory interactions can be counteracted by an increased activity of the other antagonistic pathway, according to our dose-dependent competition hypothesis (this is encoded in the corresponding Boolean rules, see Table 1).

To test our model, we ran simulations and compared the resulting stable states with the wild-type pattern observed experimentally. We applied different sets of values for the input nodes, each set corresponding to a specific territory of the ectoderm (Table 1; Materials and Methods). For the relevant combinations of active inputs, the resulting stable states are compared with the list of marker genes expected to be expressed in the corresponding territory, based on in-situ hybridisation experiments. We first considered the combinations of inputs corresponding to early 32-cell stage embryo signalling (Fig. 5), which corresponds to the onset of specification of the ventral organiser, forming at the opposite side of the gradient of Panda mRNA (Haillot et al., 2015). This pattern is correctly recapitulated by the stable states obtained for the wild-type situation (Fig. 5). The resulting stable states were then used to specify the input conditions reflecting a later blastula stage, after diffusion and shuttling of zygotic factors have taken place (Fig. 5). Indeed, as multiple diffusion events occur, some model inputs expressed in one territory are active in a broader region for blastula simulations (Fig. 6). After some iterative refinements of the rules, model simulations qualitatively recapitulated the expression patterns expected for each individual territory (ventral, ciliary, dorsal) (Fig. 6). We can therefore conclude that the regulatory graph shown in Fig. 3 supplemented by the logical rules of Table 1 are sufficient to account for the specification of the three main ectodermal D/V patterns of the sea urchin embryo.

To further validate and explore the properties of our model, we simulated loss- and gain-of-function experiments (mRNA or morpholino injection) by restricting the range of reachable levels for one or multiple node(s), e.g. to zero for a loss-of-function, or to a higher value for a gain-of-function (see Materials and Methods). Similarly to the wild-type conditions, we challenged our model by comparing the resulting stable states with the in-situ patterns observed following mRNA or morpholino injection in the embryo at an early stage. For seven of the eight morphants simulated, the model returned a unique single stable state in each region, which qualitatively matched experimental observations (Fig. 6).

In the cases of nodal morphants and of lefty mRNA overexpression, the ventral cascade fails to be established, leading to the absence of both Nodal and BMP2/4 pathway activities, and to the presence of a default ciliary state in all the ectodermal cells. Following nodal mRNA overexpression, the competition between Smad2/3 and Smad1/5/8 for Smad4 turns in the advantage of the ventral cascade, producing a fully ventralised embryo. The same pattern is obtained for the lefty morpholino, because this perturbation impacts the diffusion-repression mechanism controlled by Lefty (Juan and Hamada, 2001), enabling nodal expression to propagate without restriction (Duboc et al., 2008). In the case of overexpression of the dorsal fate repressor chordin, or in the case of BMP2/4 morphants, the absence of BMP2/4 signalling fosters a ciliary band state in the presumptive dorsal territory. However, as BMP2/4 is not necessary for the expression of nodal, genes belonging to the ventral cascade maintain a wild-type expression pattern in these morphants. Finally, following bmp2/4 overexpression, the competition for the common Smad is driven toward the activation of the dorsal cascade, giving a fully dorsalised ectoderm state.

Interestingly, in the case of the chordin morpholino, the model returned two stable states (denoted by the blue and green half circles at the bottom of Fig. 6) in the presumptive ventral region, corresponding to ventral and dorsal fates, respectively. This situation is further investigated using a probabilistic framework hereafter.

Using a probabilistic extension of the logical framework, one can unfold the temporal dynamics of the regulatory network for given initial conditions, and estimate the prevalence of a given cellular fate when the model predicts alternative stable states. In the case of the ventral region in chordin morphant conditions, we have seen that our model predicts two different stable states, corresponding to ventral and dorsal expression patterns. Using the software MaBoSS (https://maboss.curie.fr; Stoll et al., 2017), we performed stochastic temporal simulations of our model to generate mean time plots and estimate the probability to reach each of these stable states. In the absence of precise kinetic information, we first used equal rates for all (up or down) transitions.

In the wild-type ventral region, as expected, all stochastic simulations gave rise to a ventral expression pattern (Fig. 7A). In contrast, for the chordin morphant, the dorsal state is reached about twice as often as the ventral state (Fig. 7B). In other words, the dorsal pathway is more likely to win the competition for Smad4. This partial dominance of the dorsal pathway matches the weak dorsal patterns observed experimentally for chordin morphants (Fig. 6), which presumably result from the co-activation of the two antagonistic pathways. In addition, the prediction of a transient expression of ciliary genes in the ventral region (Fig. 7) is further supported by experimental evidence from Saudemont et al. (2010).

In the chordin morphants, BMP signalling diffuses unrestricted in the ventral ectoderm, promoting dorsal fates and repressing the ventral cascade. However, as Nodal is essential for bmp2/4 activation, nodal downregulation in turn leads to the repression of BMP2/4 signalling. Therefore, in the absence of Chordin, both the ventral and dorsal cascades are activated. This conclusion is supported by the experimental observation of transient Smad1/5/8 signalling and tbx2/3 expression in the ventral ectoderm (Lapraz et al., 2009).

To further assess whether this imbalance in favour of the dorsal pathway activation is sensitive to kinetic (transition) rates, we ran stochastic simulations with different Smad activation rates for the two pathways. An imbalance in favour of the dorsal Smad activation further increases the final proportion of dorsal fate compared with wild type (Fig. 7C). Conversely, an imbalance in favour of the ventral Smad activation inverses the relationship, with a higher fraction of ventrally specified states compared with dorsally specified states, almost mirroring the ratios obtained with equiprobable transition rates (Fig. 7D).

In conclusion, in the condition of a simultaneous activation of both ventral and dorsal cascades without inhibition (i.e. chordin morphants), the outcome of the competition between the BMP2/4 and Nodal pathways is at least partly driven by the kinetic rates of Smad activation. In such a situation, the pathway-specific Smad firstly activated more quickly reaches a concentration level sufficient to repress the other pathway by pre-empting Smad4 and thereby fosters the corresponding state stable.

In the preceding section, we simulated the behaviour of cells for each of the three different presumptive territories by selecting appropriate combinations of signalling input levels, which were considered as fixed during the whole duration of the simulations. To model more precisely the production and diffusion of signalling molecules across the ectoderm, we used the software EpiLog (https://epilog-tool.org) (Varela et al., 2019), which supports simulations of an epithelium encompassing multiple cells connected by diffusing signals. The behaviour of each cell of the epithelium is modelled by the same cellular logical model, but levels of input signals directly depend on the output signal values from neighbouring cells. The input signals perceived by a given cell are integrated into logical diffusion rules, which are updated synchronously during simulations (Materials and Methods; Table 2; Fig. 8A).

To simulate the wild-type situation, we initialised the model with a broad moderate (level 1) expression of Nodal and Panda in the presumptive ciliary and dorsal territories (corresponding to the six right-most range of cells, Fig. 8B), as observed experimentally (Haillot et al., 2015). The simulation correctly recapitulates the expected contiguous ventral, ciliary and dorsal ectoderm territories (Fig. 8C). This result suggests that the relatively simple diffusion rules properly account for the dynamics of the proteins governing the D/V patterning. In addition, this result highlights the crucial role of Nodal to direct the specification of the ectoderm along the whole dorsal-ventral axis. Interestingly, in the course of wild-type simulation, we can clearly distinguish the main steps of the D/V patterning, with the initial restriction of Nodal expression on the ventral side, followed by the diffusion of BMP2/4 signalling toward the dorsal side, until it activates the dorsal marker genes (Fig. S1).

As in the case of the unicellular model, we can apply specific perturbations to assess their impact at the tissue level. As shown in Fig. 8C, our multicellular simulations accurately recapitulate the phenotypes of the different morphant patterns observed experimentally. Note that the chordin morpholino has been discarded from these simulations, as it gave rise to two different stable states, which cannot be properly simulated with the EpiLog deterministic updating approach.

Using EpiLog, it is further possible to perform local perturbations by modifying the initial levels of one or several signalling molecules at specific epithelium locations. Using this feature, we could recapitulate in silico the results of an experiment reported by Lapraz et al. (2015), who injected nodal mRNA into two opposite cells of a four-cell-stage nodal knockdown embryo (i.e. following a nodal morpholino injection in the egg). This experiment triggered the formation of an ectopic inverted D/V axis and resulted in the development of Siamese pluteus larvae with two ventral sides, two ciliary bands and a central dorsal territory. Using EpiLog and imposing nodal and smad2/3/4 activity at the initial state in both the ventral and the dorsal side of the epithelium, our spatial logical simulation qualitatively recapitulates the Siamese pattern observed experimentally (Fig. 8D).

GRNs integrate documented interactions between transcription factors, signalling components and their target genes, which ultimately translate the information encoded into the genome into morphological form and body plan. However, as our delineation of developmental systems progresses, we are facing increasingly large and complex networks, which cannot be fully and rigorously understood without proper formalisation. This is clearly the case for the GRN governing D/V patterning of the sea urchin embryo, which relies on numerous signalling and regulatory factors, involved in multiple positive and negative feedback circuits.

In our modelling study, several key choices had to be made. As little is known regarding detailed mechanisms and kinetic parameters, we opted for a qualitative logical formalism. However, to properly model morphogen diffusion and dose-dependent effects, we considered a multilevel extension of the classical Boolean framework. Importantly, in the course of its conception, the model was systematically tested through extensive simulations of wild-type and perturbed conditions. In wild-type conditions, our unicellular model fully recapitulated each territory pattern independently. We further took advantage of a recent multicellular extension of the logical framework to explicitly simulate spatial pattern formation, the results of which can be more easily compared directly with the phenotypes of wild-type and morphant embryos.

In the process of model building, we identified several parameters that were crucial for proper simulation of the D/V patterning. For the unicellular model, the relative action of the multilevel nodes was a key component to properly model weak versus strong repressions (e.g. Nodal overcomes Lefty repression at level 2 but not level 1, Chordin overcomes Tolloid repression at level 2 but not level 1). Considering intermediate activity levels enables the writing of refined rules for the effect of interaction combinations on downstream nodes. For the multicellular model, the logical rule governing Chordin maximal diffusion distance was the key parameter modulating the size of the dorsal ectoderm region, hence suggesting that it may be a tightly controlled mechanism in the embryos.

A key step in our study was to model the interplay between the Nodal and BMP pathways. In this respect, we were guided by our experiments dealing with the treatment of embryos with recombinant Nodal or BMP2/4 proteins at blastula stage (i.e. after the initial specification of the ventral and dorsal territories). These experiments demonstrated that overactivation of one of these pathways is sufficient to abrogate signalling from the other pathway. Consequently, our model highlights the strong antagonism between Nodal and BMP2/4 signalling and suggests that the outcome of this competition relies at least in part on the relative doses of the two TGF-βs, rather than or in addition to their timing of activation. We further suggest that the antagonism stems from a competition for the recruitment of Smad protein complexes. Although evidence for binding sites of Smad complexes on the cis-regulatory regions of BMP2/4 and Nodal target genes has been described (Hill, 2016), the mechanism of competition between the two cascades has not been fully clarified yet to our knowledge, calling for further experiments to support and validate this hypothesis.

This competition between Nodal and BMP2/4 plays a key role in understanding the regulatory dynamics within the chordin knockdown experiments, which was the only morphant not fully recapitulated by our model for all three territories. In the case of the chordin knockdown, our logical model predicted that both the ventral and dorsal steady states were possible in the presumptive ventral region. Accordingly, in the chordin morphant, both Nodal and BMP2/4 pathways are activated, antagonising each other. Following this ectopic activation of BMP signalling, the ventral territory in chordin morphants displays a transient dorsalisation, before reversing towards a ventral ectoderm fate during gastrulation, as shown by the presence of a mouth opening. To further explore the underlying regulatory mechanism of this dorsal fate dominance, we performed a stochastic logical simulation of the unicellular model in the chordin knockdown condition. This analysis resulted in a higher proportion of active dorsal fate over ventral fate, in agreement with experimental observations. This result suggests that the transient dorsal dominance is encoded in the structure of the GRN. Indeed, even in the case of the chordin morphants, the model accurately recapitulates the conflict caused by the coactivation of the Nodal and BMP2/4 pathways in the ventral ectoderm. However, by modulating the rates associated with the different Smads and performing additional simulation, we showed that the outcome of the competition between the two pathways is sensitive to these rates in the context of co-activation.

At this point, our model remains limited to the major early D/V patterning events occurring in sea urchin. However, this model could be tentatively extended in the future to integrate novel data and explore more refined specification and differentiation events. For example, it could be extended to investigate the specification of the boundary ectoderm region, located at the interface between the ectoderm and endomesoderm, which plays a central role in positioning the clusters of primary mesenchyme cells and in patterning of the spicules (Armstrong and McClay, 1994; Armstrong et al., 1993; Duloquin et al., 2007; Hardin et al., 1992; Röttinger et al., 2008). This process is known to depend on Wnt signalling, presumably in conjunction with Nodal, BMP2/4 and ADMP2 signalling (Lapraz et al., 2015; McIntyre et al., 2013; Röttinger et al., 2008; Saudemont et al., 2010). With the current unicellular model, the simulation with the input levels corresponding to the boundary ectoderm (i.e. Admp2 and Wnt active) results in a dorsal stable state (Fig. S2).

Another possible extension to the model would be the integration of the negative feedback of Smad6 on the Nodal and BMP2/4 pathways. Indeed, Smad6 is activated by the dorsal signalling downstream of BMP2/4 (Massagué and Chen, 2000). Such a negative feedback circuit typically generates an oscillatory behaviour (Thomas and D'Ari, 1990). In the case of our model of the D/V GRN, adding this negative feedback circuit results in a cyclic activation of BMP and Smad6 nodes. We chose to not include this element in our current model as we still lack a sufficient understanding of the underlying mechanisms.

Similarly, noggin and follistatin (two inhibitors of BMP signalling) have been additionally identified in the sea urchin genome (Lapraz et al., 2006). Preliminary data regarding the regulation of their expression suggests that these genes are expressed at a very low level during early development. Hence, we have decided to leave these genes aside at this stage. As we cannot exclude that these genes nevertheless may play a role in the repression of the dorsal patterning of the ectoderm, we believe that their integration to the model will have to be considered in future versions to refine our understanding of the Nodal/BMP antagonism.

Our logical model focuses on the blastula and gastrula stages of sea urchin embryogenesis. One possible extension would be to further explore the regulatory interactions taking place at earlier stages. In the case of the 32-cell stage, our model correctly recapitulates wild-type pattern mainly driven by Panda expression. Furthermore, the simulation results of Panda loss-of-function in 32-cell-stage conditions mirror the fully ventralised phenotype obtained experimentally (Fig. S3A,B;Haillot et al., 2015). However, the simulations for Panda overexpression show discrepancies relative to the experimental observations. Indeed, our model predicts a loss of ventral fate specification, whereas global injection of Panda mRNA does not impact the wild-type pattern. Current models indicate that an asymmetry of Panda mRNA provides the spatial cue that in turn controls the polarised activation of downstream genes. Indeed, rather than a response to local absolute concentration of panda mRNA, the symmetry breaking is initiated by the detection of a relative change in panda mRNA concentration along the D-V axis. Such signalling based on a multicellular gradient cannot be currently recapitulated by our unicellular model, as it requires the integration of inputs from multiple surrounding cells and also the reliance on relative differences in concentration instead of absolute levels. As relative concentration differences are difficult to model with the logical formalism, a continuous framework (e.g. ODE) would be better suited to further explore the specificities of such multicellular gradient signalling.

To conclude, using a qualitative logical model we could capture several salient dynamical features of the GRN governing the early D/V patterning of sea urchin embryos, including the key role played by intercellular interactions. More precisely, our model integrates three key mechanisms responsible for the proper ectoderm patterning. First, the early initiation of Nodal activity is required to robustly maintain a ventral organiser through a positive feedback circuit. Second, the diffusion of dorsalising signal from ventral to dorsal regions is crucial for the onset of the dorsal patterning. Third, the antagonism between ventral and dorsal signalling pathways is the key mechanism shaping the territory boundaries along the D/V axis of the sea urchin embryo. In addition, the kinetic analysis of the model structure suggests the dominance of the BMP pathway over the Nodal pathway, in the context of simultaneous activation.

More generally, we have shown that qualitative models are useful to explore the interplay between maternal factors and zygotic genes, which orchestrates patterning of the ectoderm of the sea urchin embryo downstream of intercellular signals. To ease further model analyses and extensions, we provide our models together with a Jupyter notebook implementing all the simulations performed with GINsim and MaBoSS (see Data Availability section), and a Docker image containing the necessary software (https://github.com/colomoto/colomoto-docker).

Adult sea urchins (P. lividus) were collected in the bay of Villefranche-sur-Mer, France. Embryos were cultured as described in Lepage and Gache (1989, 1990). Fertilisation envelopes were removed by adding 1 mM 3-amino-1,2,4 triazole (ATA) 1 min before insemination to prevent hardening of this envelope followed by filtration through a 75 μm nylon net. Treatments with recombinant BMP4 (recombinant mouse BMP-4 protein, CF 5020-BP/CF, R&D Systems) or Nodal (recombinant human Nodal protein, 3218-ND-025/CF, R&D Systems) proteins were performed at the time indicated in the schemes by adding the recombinant protein diluted from stocks in 1 mM HCl in 24-well plates containing ∼1000 embryos in 2 ml of artificial sea water (Saudemont et al., 2010). All the experiments described in this study have been repeated two or three times. At least 200 wild-type and 50 injected embryos were analysed for each condition or experiment. In the case of treatments with recombinant proteins, more than 500 embryos were scored for morphological phenotypes for each condition and ∼250 embryos were used for in situ hybridisation. Only phenotypes observed in more than 90% of the embryos are shown.

For overexpression studies, capped mRNAs were synthesised from NotI-linearised templates using the mMessage mMachine kit (Ambion). After synthesis, capped RNAs were purified on Sephadex G50 columns and quantitated by spectrophotometry. The nodal, lefty, chordin and bmp2/4 pCS2 constructs have been described in Duboc et al. (2004, 2008) and Lapraz et al. (2009). In this study, nodal mRNA was injected at 400 µg/ml, lefty mRNA at 200 µg/ml, chordin mRNA at 1000 µg/ml and bmp2/4 mRNA at 500 µg/ml. Approximately 2-4 pl of capped mRNA was injected mixed with Tetramethylrhodamine Dextran (10,000 MW, Molecular Probes, #D1817) at 5 mg/ml (Duboc et al., 2004). Morpholino oligonucleotides were dissolved in sterile water together with Tetramethylrhodamine Dextran (10,000 MW) at 5 mg/ml and approximately 2-4 pl of the resulting solution was injected at the one-cell stage. nodal morpholino was used at 1 mM (Duboc et al., 2004), chordin and bmp2/4 morpholinos were injected at 1 mM and 0.2 mM, respectively (Lapraz et al., 2009), lefty morpholino at 1.5 mM (Duboc et al., 2008) and panda morpholino at 1.2 mM (Haillot et al., 2015).

Embryos were fixed with 4% formaldehyde for 15 min at swimming blastula stage (3 h after adding BMP2/4 protein) then briefly permeabilised using methanol. Anti-Phospho-Smad1 (Ser463/465)/Smad5 (Ser463/465)/Smad9 (Ser465/467) (D5B10, 13820, Cell Signaling Technology) was used at 1/400. Embryos were imaged using an Axio Imager.M2 microscope (Zeiss).

In situ hybridisation was performed using standard methods (Harland, 1991) with DIG-labelled RNA probes and developed with NBT/BCIP reagent. The nodal, chordin and tbx2/3 probes have been described previously (Duboc et al., 2004; Lapraz et al., 2009). Control and experimental embryos were developed for the same time in the same experiments. Embryos were imaged with an Axio Imager M2 microscope.

We built our model using the multilevel logical formalism introduced by R. Thomas (Thomas and D'Ari, 1990). This qualitative approach relies on graph-based representations of the network and of its dynamics. The network is formalised as a regulatory graph, in which nodes denote molecular species (e.g. proteins or genes), whereas signed arcs denote regulatory interactions, positive or negative. The nodes can take a limited number of integer values, only two (0 or 1) in the simplest Boolean case, but potentially more when biologically justified, for example in the case of morphogens with clearly distinct activity ranges. Hence, each regulatory arc is associated with an integer threshold, always 1 in the Boolean case, but potentially higher in the case of a multilevel regulator. Boolean rules further specify how each node reacts to the presence or absence of the corresponding incoming interactions. Specific (non-overlapping) Boolean rules are defined for each value of each node. In the case of nodes with unknown activator (e.g. Univin), a default basal level is fixed before the simulation (Table 1).

In this example, the regulatory interactions from Brachyury to FoxA and from Repressor_R1 to FoxA correspond to an activation and to an inhibition, respectively. The levels of the input (unregulated) nodes are defined at the start of simulations.

Using the Boolean rules of Table 1, we can simulate the behaviour of the system for different combinations of input values. In this respect, we use the asynchronous updating approach proposed by R. Thomas (Thomas, 1991), which consists in following all the different possible single unitary changes of values induced by the rules.

The dynamical behaviour of the model is generally represented in terms of a state transition graph, in which each node represents a discrete state of the model (i.e. a vector listing the values of the different components of the model), whereas each arc represents a state transition (i.e. the unitary change of the value of the component). In this work, we took advantage of the implementation of this logical formalism into the user-friendly Java software suite GINsim (version 3.0, http://ginsim.org; Naldi et al., 2018b). In our analyses, we particularly focused on stable states (see e.g. Fig. 5B), which typically denote cell differentiation states. These can be directly computed (i.e. without unfolding the state transition graph) using a very efficient algorithm implemented in GINsim (Naldi et al., 2007).

We simulated the behaviour of each D/V region independently, considering different sets of fixed values for the input nodes in the ventral, ciliary and dorsal presumptive territories. These sets of input values were defined based on previously published results (see Results section, Model building). To account for their transient activity, the input nodes for Panda and Admp2 were defined as having a basal value of 0, enabling them to turn off during the simulation if active in the initial state conditions.

As we simulate each territory individually, the unicellular model cannot directly take into account the diffusion of morphogens, which are therefore specified as input levels (e.g. the presence of Lefty is considered as an active input in the ciliary regions, although it is known that it diffuses from the ventral region). For each simulation, we extract the resulting stable state(s) and classify them as ventral, ciliary or dorsal pattern depending on output node levels. For example, the initial conditions for the simulation of the ventral ectoderm territory amount to considering the inputs Nodal (level 2), Lefty, Admp1, Tolloid, Chordin (level 2) and BMP2/4 as active (whereas the other inputs are inactive). This combination results in a stable state with all the ventral nodes active and the dorsal nodes inactive. In contrast, when the initial conditions are set with Nodal (level 1), Lefty, BMP2/4, Chordin (level 1) and Tolloid being active (and the other inputs inactive), the resulting stable state corresponds to the dorsal fate.

Genetic perturbations are defined in GINsim by constraining the values of selected nodes of the model. To simulate a knockdown morphant (e.g. injection of a morpholino), the level of the corresponding node is set and maintained at 0. In the case of an ectopic expression (e.g. injection of an mRNA), the level of the corresponding node is set and maintained at its maximal value, which can be 1 or higher, in the case of a multilevel node. Morphogen diffusion is taken into account through the specification of proper input values, which thus need to be adjusted for each morphant. For example, the ectopic activation of Nodal is known to induce the activation of its downstream target BMP2/4 very early on; hence, the corresponding input variables must be set at their highest levels for the simulation of ectopic nodal expression.

When several stable states can be reached (as in the case of chordin morpholino), we have performed probabilistic simulations to evaluate the probability of reaching each of these stable states from the specified initial conditions. In this respect we used the software MaBoSS (https://maboss.curie.fr), a C++ software enabling the simulation of continuous/discrete time Markov processes, applied to Boolean networks. The original unicellular model is converted into the MaBoSS-compliant format using a specific export functionality of GINsim, which involves the replacement of multilevel nodes by sets of Boolean variables, without affecting the model dynamic (Stoll et al., 2017). Per default, all up and down rates are considered equal, but these can be modified at will.

In this study, we used MaBoSS to simulate the chordin morpholino perturbation (comparing it with the wild-type situation), which resulted in two possible stable states in the unicellular model. The inputs were fixed as for the ventral configuration (Nodal, Lefty, BMP2/4 and Admp1 active) in the presence or inactivation of Chordin. We then modified the propensity to activate the ventral or the dorsal cascade by adjusting the ratios of the rates assigned to the Alk receptors corresponding to each of the two cascades: 0.5/0.5 (equiprobable rates), 0.75/0.25 (ratio favouring the dorsal Alk), 0.25/0.75 (ratio favouring the ventral Alk).

We took advantage of the recent software EpiLog (https://epilog-tool.org, v1.1.1; Varela et al., 2019) to perform multicellular simulations. Our epithelium model is nine cells wide and six cells long, made of hexagonal shaped cells, each one being in direct contact with at most six different neighbouring cells. The top and bottom parts of the epithelium are wrapped together to allow diffusion of signalling molecules through these two sides. Each hexagonal cell encompasses one copy of our unicellular model and thus behave according to the same logical rules during simulations. In contrast with our previous unicellular simulations, the cell inputs are dynamically updated based on the signals perceived from the corresponding output nodes of neighbouring cells (e.g. Nodal input level will be based on the quantity of neighbouring cells expressing Nodal as an output). The rules integrating the extracellular signals are identical for all cells of the epithelium. In our epithelium simulations, input nodes of all cells are updated in a synchronous manner. Hence, each epithelium simulation gives rise to a deterministic trajectory ending in a single attractor at the level of the whole tissue (a stable state for the simulations reported here).

The dynamical update of the input node levels implies the definition of additional logical rules for the diffusion of signals, e.g. of the values of input nodes depending on the output nodes active in neighbouring cells, taking into account their distance from the target cell. For example, the rule ‘{Nodal:2[1:], min=1}’ states that a cell will have its Nodal input node value converging toward the value 1 if at least one cell is expressing Nodal at a value of 2 at a minimum distance of one cell (i.e. all cells except the target cell itself).

For our epithelium simulations, we define the initial state by selecting the nodes that will be active in a specific set of cells at the start of the simulation. During simulations, the values of these nodes can change depending on the model state and on paracrine signalling. To simulate a wild-type embryo, we set the model to an initial state where nodal is broadly expressed at a moderate level (level 1), and where panda is expressed in the presumptive ciliary and dorsal territories (six right-most cell columns of the epithelium; Fig 8B), with the initial and transient activation of Panda and Nodal (output) nodes. Univin is also ubiquitously present at initial state. As in the unicellular simulation, to simulate loss- or gain-of-function perturbations, the value of the corresponding node is set and maintained at a fixed value.

To simulate the double injection of nodal mRNA leading to Siamese embryos, we use the wild-type logical model, with an initial state accounting for a ventral expression of Smad1/4/5/8 and Nodal on the ventral side, but also on the dorsal side of the epithelium (right-most cell column of the epithelium), i.e. a symmetrical activity pattern (Fig. 8D).

We thank Aline Chessel for excellent technical help. We are indebted to Guillaume Lavisse for insightful comments. We thank Aurélie Martres for taking care of the sea urchins. We thank Aurelien Naldi for continuous help and support on the development of the Colomoto Jupyter Notebook. We thank Claudine Chaouiya and Pedro Monteiro for their support regarding the use of EpiLog. We thank Mathurin Dorel for his help in defining a preliminary version of the cellular model.

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