 Research article
 Open Access
 Published:
Landscape reveals critical network structures for sharpening gene expression boundaries
BMC Systems Biology volumeÂ 12, ArticleÂ number:Â 67 (2018)
Abstract
Background
Spatial pattern formation is a critical issue in developmental biology. Gene expression boundary sharpening has been observed from both experiments and modeling simulations. However, the mechanism to determine the sharpness of the boundary is not fully elucidated.
Results
We investigated the boundary sharpening resulted by three biological motifs, interacting with morphogens, and uncovered their probabilistic landscapes. The landscape view, along with calculated average switching time between attractors, provides a natural explanation for the boundary sharpening behavior relying on the noise induced gene state switchings. To possess boundary sharpening potential, a gene network needs to generate an asymmetric bistable state, i.e. one of the two stable states is less stable than the other. We found that the mutual repressed selfactivation model displays more robust boundary sharpening ability against parameter perturbation, compared to the mutual repression or the selfactivation model. This is supported by the results of switching time calculated from the landscape, which indicate that the mutual repressed selfactivation model has shortest switching time, among three models. Additionally, introducing cross gradients of morphogens provides a more stable mechanism for the boundary sharpening of gene expression, due to a twoway switching mechanism.
Conclusions
Our results reveal the underlying principle for the gene expression boundary sharpening, and pave the way for the mechanistic understanding of cell fate decisions in the pattern formation processes of development.
Background
A persisting focal issue in developmental biology is how the spatial pattern of gene expression is formed. During embryonic development, cells with different expression patterns occupy different domains separated by sharp borders. It is suggested that cell exploits the diffusive molecules and morphogens to receive positional information [1,2,3]. Morphogens, by forming concentration gradients, specify distinct cell fates and regulate patterns. So, cell fate decisions are vital for the pattern formation in the developmental process [4,5,6]. The challenge is to understand how cells form sharp gene expression boundaries from light morphogen gradients. Some mechanisms involve cell sorting resulting from cell movements and differential adhesion [7, 8]. Also, some cellintrinsic boundaryforming mechanisms, which transfer the differences in morphogen gradients into sharp change in downstream gene expression, have been proposed [9]. However, the mechanisms for boundary formation and sharpening remain not fully described.
In cells, there are intrinsic fluctuations from limited number of molecules and external fluctuations from inhomogeneous environments [10,11,12,13,14]. Therefore, the gene expression fluctuations need to be considered to simulate the real cellular environments. Some noise attenuation mechanisms have been suggested in cell fate decision processes during developmental patterning [15, 16]. Recently, Zhang et al. proposed a noise induced cell fate switching mechanism to explain gene expression boundary sharpening in the zebrafish hindbrain [15]. In their work, noise has been shown to play some critical roles in developmental pattern formation processes by facilitating gene expression state switchings. However, the underlying mechanism for how the noise promotes the sharpening of gene expression borders remains elusive.
In this work, we constructed a mutual repressed selfactivation (MRSA) model, which interacts with morphogens. The involved two genes (X and Y) are both activated by a morphogen M. In particular, we introduce the fluctuations to both morphogen gradients and gene expression levels. To explore the boundary sharpening mechanisms related to gene state transitions, we resort to a probabilistic landscape approach [17,18,19], which has been employed to study the stability of attractors for the gene networks and the switchings between the attractors. The landscape was proposed by Waddington to characterize development and differentiation of cells, as a metaphor [20]. From the landscape theory, different phenotypes can be depicted as the basins of attraction on a potential landscape surface, and the cell fate decision process can be viewed as a ball rolling down from one basin to another on the landscape surface.
We mapped out the landscape for the mutual repressed selfactivation model at different spatial locations (corresponding to different Morphogen levels). At the position close to the gene expression boundary, the landscape displays an unbalanced bistable state, where one attractor is much less stable than the other. The landscape view for the system transforms the boundary sharpening problem to the understanding of cell fate transitions between basins, and thus provides a natural explanation for the noise induced boundary sharpening behavior. We compared three common motifs  mutual repressed selfactivation, mutual repression, and selfactivation  to see the effects of network topology on the boundary sharpening ability. We found that mutual repressed selfactivation model has the better boundary sharpening performance against fluctuations, which is consistent with the results that mutual repressed selfactivation model possesses the shorter switching time calculated from the landscape among three models. Adding a second morphogen to the models provides more stable boundary sharpening ability due to a twoway switching mechanism. Our results reveal the underlying mechanisms for noise induced boundary sharpening of gene expression, and shed light on the deeper understanding of spatial pattern formation in embryonic development.
Results
Landscape explains the gene expression boundary sharpening from a mutual repressed selfactivation model
We first investigate a mutual repressed selfactivation model (MRSA), interacting with a morphogen M (Fig. 1a). In this circuit, gene X and gene Y mutually repress each other and selfactivate themselves. The morphogen M activates the expression of both gene X and gene Y. Figure 1bd show the boundary sharpening effects over time at different gene expression noise level d (coefficient of variance) from twodimensional simulations (A for dâ€‰=â€‰0, B for dâ€‰=â€‰0.01, C for dâ€‰=â€‰0.02). Here each cell is represented by a grid square. The blue grids represent the cells with gene X expressed, and the red grids represent cells with gene Y expressed. For the system without gene expression noise (Fig. 1b), the expressions for gene X and Y generate the rough boundary, due to the fluctuations in the morphogen M. If the noise level is too high (Fig. 1d), the gene expression boundary is also rough because the large fluctuations in the gene expression level dominate the expression states and blur the boundary. Surprisingly, when the noise level is in a certain medium range, the boundary becomes sharpening over time (Fig. 1c). The gene expression noise induced the state transition from X expressed state (blue) to Y expressed state (red). Therefore, around the boundary area most of cells are switched to Y state, which makes the boundary sharpened.
To elucidate the mechanism of how noise induces gene state switchings, we acquired the probabilistic landscape at different spatial positions (corresponding to different M level). Landscape is defined as Uâ€‰=â€‰â€‰âˆ’â€‰â€‰logâ€‰(P_{ ss }) [17, 18, 21,22,23], and P_{ ss } is the steady state probability distribution (see methods for how to obtain landscape). Figure 2ac show the probabilistic landscape for the expression level of gene X and Y for the MRSA model at different M levels. Here, the blue region represents the higher probability or lower potential state (attractors), and the red region represents the lower probability or higher potential state. Two stable states (basins) appear on the landscape (bistability). One of them is X expressed state (labeled by X) and the other one is Y expressed state (labelled by Y). The shape of landscape determines the transition difficulty between attractors (X and Y). To quantify the transition difficulty, we define the barrier height (Fig. 2a) as the potential difference from the local minimum (say X attractor) to the saddle point between X attractor and Y attractor. The larger the barrier height, the more difficult the system switches from X state to Y state. Figure 2ac show that the landscape shape changes as the morphogen M level increases. When M level is 0.5 (Fig. 2a), the landscape exhibits a bistable shape (two basins). As the M increases, the left basin (X attractor) becomes shallower gradually relative to the right basin (Y attractor), and the barrier for X attractor becomes lower, as shown in Fig. 2b. When M increases further (Fig. 2c), the landscape eventually becomes a monostable state (only the Y attractor exists).
The landscape change provides an explanation for the boundary formation and the boundary sharpening effects of MRSA circuit (Fig. 1). When M level is high (Fig. 2c), only Y attractor is stable as indicated on the landscape, so cells will end in the Y expression state (red grids in Fig. 1c). When M level is low (Fig. 2a), gene expression of X and Y form a bistable state, and cells will end in X expression state when assuming gene X is initially expressed [15]. For the medium M level (the central position), landscape also exhibits a bistable state, but X attractor is much less stable than the Y attractor. Induced by certain level of fluctuations, cells will switch from X attractor to Y attractor, which explains why the gene expression boundary will be sharpened. As shown in Fig. 1c, around the boundary position the blue grids change to red grids as T increases from 10 to 1000, which sharpens the boundary. The landscapes obtained here also explain why the noise level has to be in certain range for inducing the boundary sharpening. That is because if the noise level is low it is hard to trigger the switchings from X attractor to Y attractor (cell can not escape from X basin). But if the noise level is too high, the switchings can happen in either direction, meaning that cells can change state from X to Y (jump from X basin to Y basin) but can also change from Y to X (jump from Y basin to X basin). In that case, the boundary sharpening will not take place.
During the development of rhombomeres in the zebrafish hindbrain, the morphogen retinoic acid (RA) induces expression of gene hoxb1a in rhombomere 4 (r4) and gene krox20 in r3 and r5. Fluorescent in situ hybridization reveals rough edges around these gene expression domains, in which cells coexpress hoxb1a and krox20 on either side of the boundary, and these sharpen within a few hours [15]. It is found that the boundary sharpening effects are not clear in the r3/r4 boundary from a similar mathematical model [15]. The landscape picture we acquired here provides an explanation for this phenomenon. The r3/r4 boundary is corresponding to the low level M position in our model, which is corresponding to a bistable landscape (Fig. 2a). Here, both two basins are stable, not like the case in the Fig. 2b in which one basin is much less stable than the other. The noise induced transition will not help the sharpening of the boundary because the switching is bidirectional. Cell can switch from X basin to Y basin or from Y basin to X basin. But from population level, the cell type switchings are averaged out, and will not lead to sharpened boundaries.
To further explore the gene state switching mechanisms, we calculated the mean first passage time (MFPT) from state X to state Y, characterizing the average switching time. Now that the boundary sharpening problem has been transformed to a cell fate transition problem, we suppose that the boundary sharpening time should be correlated to the gene state switching time. We calculated the MFPT from the trajectories based on the solution of stochastic ordinary differential equations (see methods). The parameter A represents the ratio of the synthesis rate for the two genes, which measures the asymmetry level of the system. The parameter R represents the strength of the mutual inhibition between gene X and Y. Figure 2df show the average switching time from X state to Y state at different noise level, separately for different M, different A, and different R. Here, every curve (color) represents one set of parameter value. Each curve shows that given a noise level what is the average switching time or given a time interval how large noise is needed to trigger the switchings. As the noise level increases, the average switching time declines (Fig. 2df), because larger noise promotes the transition between attractors. Also, with M level increased the average switching time decreases (Fig. 2d), i.e., the switching from X state to Y state becomes faster. This is consistent with the landscape analysis because, as M increases, X (Y) attractor becomes less and less (more and more) stable (Fig. 2c), making the transition from X to Y easier. Figure 2e and f show that as the asymmetry level A increases or as the mutual repression strength R decreases, the switching from X state to Y state becomes faster. This indicates that a stable gene state switching for the MRSA model requires larger asymmetry between gene X and Y, and smaller mutual repression strength between gene X and Y. Therefore a natural prediction from these results is that the boundary sharpening will be faster for a circuit with more asymmetry or with weaker mutual repression.
In Fig. 2g and h we compared the MFPT and the barrier heights at different M level. The barrier height quantifies the change of the landscape topography (Fig. 2ac). The barrier height increases as the switching time increases (Fig. 2h), because a larger barrier makes the switching from X attractor to Y attractor more difficult and thus a longer switching time.
We also calculated the results for different amplitude ratio (Lâ€‰=â€‰Ïµ/d) between morphogen noise level Ïµ and gene expression noise level d for MRSA model. As we can see from Figs. 1 and 3, the gene expression noise level d has to be in an appropriate range to generate boundary sharpening effects. If the gene expression noise level d is too big, the boundary sharpening effects will not occur no matter how we choose the value of morphogen noise level or the amplitude ratio L (Fig. 3df). Given an appropriate value of d (Fig. 3ac, dâ€‰=â€‰0.01), the morphogen noise level or the amplitude ratio L also needs to be tuned. If the morphogen noise level or the amplitude ratio L is too big, boundary is not sharpened (Fig. 3c). If the morphogen noise level or the amplitude ratio L is too small, the boundary keeps sharpened most of the time except for a quick transition period from initial uneven distributions (Fig. 3a). We also provide the quantitative comparisons of boundary sharpening effects over time from sharpening index (SI), for MRSA models at different morphogen noise level Ïµ and gene expression noise level d value (Additional file 1: Figure S1), which provide quantitative support to our above conclusions.
To see how the vertical resolution influence the effects of the boundary sharpening, we also double the number of the vertical grids and made the simulations for the MRSA model (Additional file 1: Figure S2). Our results show that increasing the grid number in the vertical direction does not influence the major conclusions (Additional file 1: Figure S2).
Mutual repressed selfactivation (MRSA) model is more robust for generating gene expression boundary sharpening than selfactivation (SA) model or mutualrepression (MR) model
To explore how the boundary sharpening effects depend on the topology of the circuit, we also investigated the other two models, selfactivation (SA) model (Fig. 4j) and mutualrepression (MR) model (Fig. 5j). Simulations show the boundary sharpening effects over time for the selfactivation model at different noise level (Fig. 4ac). When noise level is too low or too high, the system tends to form rough boundary (Fig. 4a and c), whereas for the medium level of noise the boundary is sharpened with time (Fig. 4b). Similar to the case for the MRSA model, the boundary sharpening effects of SA model can be explained from the probabilistic landscape perspective. With M level increased (Fig. 4df), the landscape changes from a relative balanced bistable state (X and Y attractor coexist) to a Y inclined bistable state, and finally to a monostable Y state. For the medium M level, the noiseinduced transitions make cells switch from X state to Y state, and thus lead to the boundary sharpening. For selfactivation model, a critical parameter is the selfactivation strength k. We calculated the change of average switching time with noise level at different k. When the selfactivation strength k is larger, the switching is faster, indicating a stronger selfactivation promotes gene switching from X state to Y state, and therefore leads to a faster boundary sharpening effect.
Figure 5ac show the boundary sharpening effects over time for the mutualrepression model (Fig. 5j) at different noise level. We also observed the boundary sharpening effects at middle level of noise from simulations (Fig. 5b). Similarly, the landscape (Fig. 5df) provides corresponding explanations for the boundary sharpening mechanisms.
To make comparisons for three models, we define a sharpening index (SI) to quantify the boundary sharpening effects of the systems. The sharpening index is defined based on the number of the grid columns with mixed expression of X and Y around the boundary. A smaller SI means a sharper boundary. For each of the three models, we gave parameters a perturbation amplitude (Ïƒâ€‰=â€‰0.1), and obtained an average SI from multiple simulations. Figure 6a shows the comparison results for three models. For all three models (three curves in Fig. 6a), SI decreases with time and almost reaches a steady value. This indicates all three models have the potential with the boundary sharpening ability at certain parameter regions. However, we found that the SI for the MRSA model reaches a lower value than the other two models. To validate this result, we also calculated the switching time from the landscape at different noise level and different parameter regions, and made comparisons for three models (Fig. 6b). In Fig. 6b, the red, blue, green curves represent MRSA model, MR model and SA model, separately. Three different style of curves represent three typical parameter choices for each of the three models. Although it is not easy to make comparisons directly for three models due to the limit of specific parameter choices, it appears that the MRSA model in general has the shortest switching time (the red lines are inclined to be lower than blue and green lines) among three models.
From the landscape description of the system, the boundary sharpening problem has been transformed to the understanding of state switching from X basin to Y basin, i.e. the boundary sharpening time should be roughly correlated with the state switching time between X basin and Y basin, and shorter switching time means faster boundary sharpening. Our results show that the MRSA model provides a better boundary sharpening performance than the other two models (smaller sharpening index as shown in Fig. 6a, and shorter switching time as shown in Fig. 6b). From the topological perspective, this is because MRSA motif provides a stable structure for bistability, hysteresis, and ultrasensitivity, which many biological systems employ to achieve different functions. For example, MRSA motif has been exploited in stem cell developmental system [24, 25], cancer stem cell system [23], and EMT system [26], to generate multistable states.
Gene expression boundary sharpening for cross morphogen gradients
It was proposed that cells can also respond to multiple morphogen signals. For example, a BMPFGF interacting morphogen system was shown to be critical to the pattern formation in the development of forebrain [27]. Here we ask how such morphogen gradients influence the boundary sharpening effects. By adding a second gradient of morphogen to the MRSA model (Fig. 7a and b), we observed the boundary sharpening effects from simulations (Fig. 7c and d).
We found that the boundary sharpening can take place in two different situations (Fig. 7c and d). One of them is gene expression boundary is sharpened because of the switchings from state Y to state X (red to blue, Fig. 7c), and the other is gene expression boundary is sharpened because of the switchings from X state to Y state (blue to red, Fig. 7d). The difference of these two cases is that they have different initial conditions for X and Y expression. The mechanisms for this twoway boundary sharpening can be revealed from the landscape shown in Fig. 7eg. The landscape in the central position (Fig. 7f) displays a bistable state, and two attractors have the similar stability, whereas for the left position (Fig. 7e) the landscape exhibits a bistable state biased to X attractor, and for the right position (Fig. 7g) the landscape exhibits a bistable state biased to Y attractor. So, if the initially Y is chosen to be more expressed then X (initial condition is chosen closer to Y attractor), the boundary will appear in a position between Fig. 7e and g, generating the boundary sharpening behavior as shown in Fig. 7c. On the contrary, if the initial condition choice is X is more expressed than Y (initial condition is chosen closer to X attractor), the boundary will appear in between Fig. 7f and g, generating the boundary sharpening behavior as shown in Fig. 7d.
For the case with single morphogen gradient, a requirement for the appearance of boundary sharpening effects is the appropriate choice of initial conditions (X needs to be more expressed than Y, Fig. 2) [15]. One advantage for the system with cross morphogen gradients is that it can trigger boundary sharpening effects regardless of the choice of initial conditions, i.e., the requirements for the initial conditions can be relaxed. Therefore, the cross morphogen system provides a more stable mechanism for the boundary sharpening of gene expression domain.
Discussion
Cells employ morphogen gradients to control expression of different genes and form distinct spatial patterns. One key issue is how the boundary between gene expression domains is generated and sharpened. Studies suggest that morphogens interacting with downstream gene regulatory networks lead to ultrasensitivity and border formation [28,29,30]. We investigated such genemorphogen interaction networks under fluctuations. Previous works showed that fluctuations, commonly regarded as detrimental to the robustness of regulatory networks, may play a critical role in pattern formation process, and appropriate noise level may promote the sharpening of gene expression boundaries [15]. Here, we aim to disclose the underlying mechanisms for noise to promote the boundary sharpening. By uncovering the probabilistic landscape of a morphogengene interaction network, we found that as the morphogen level M increases, the landscape changes from a relative balanced bistable state to a biased bistable state, and finally to a monostable state. The boundary appears in the position corresponding to the biased bistable landscape. Starting from the less stable state of the bistable states, noise will induce the state transition from the less stable state to the relative stable state. This provides a natural explanation for the noiseinduced boundary sharpening mechanism.
From the analysis for the boundary sharpening ability of three circuits, we propose that to possess the boundary sharpening potential for a gene regulatory circuit with a single morphogen gradient, three conditions need to be met:

1.
The system needs to form a bistable state system.

2.
The bistable state needs to be asymmetric, meaning that one of the two stable states (say X state) is less stable than the other stable state (say Y state).

3.
Gene X (corresponding to the less stable state) needs to be expressed initially, i.e. the initial condition should be chosen at a position closer to X attractor rather than Y attractor.
As can be seen, some basic gene regulatory motifs such as MRSA, SA, and MR, are all able to generate bistable states with appropriate parameter choices. So, the condition 1 is easy to be satisfied. However, only bistability is not sufficient to trigger boundary sharpening for the system. Only if the bistability is inclined to one of the two states (parameters need to be tuned), i.e. one of the two states is very unstable, the noiseinduced switchings take place and lead to the boundary sharpening for the pattern formation system.
To explore the influence of the topology on the boundary sharpening ability of the network, we compared the boundary sharpening ability for three models: mutual repression with selfactivation model, selfactivation model, and mutual repression model. We define the sharpness index (SI) to quantify the boundary sharpening ability of the circuit. We found that MRSA model possesses the most stable boundary sharpening ability against parameter variations among the three models, which is supported by the results of average switching time indicating that the MRSA model has the shortest average switching time among three models. This might be because, compared to mutual repression or selfactivation, MRSA model is more robust to generate bistable states or multiple states under perturbations, and thus possesses better boundary sharpening abilities. We also introduced a second morphogen gradient to explore the effects of cross morphogen gradients on the boundary sharpening of the system. We found that in the cross morphogen gradients model the boundary sharpening can appear independent on the initial concentrations of gene X and gene Y, because of a twoway switching mechanism. For single gradient model, in order to induce the boundary sharpening, gene X has to be expressed more than gene Y. Therefore, adding second gradient can make the boundary sharpening behavior more universal, providing a more stable mechanism for the sharpening of gene expression boundary.
As suggested by recent studies, the gene expression boundaries are formed and sharpened in vivo [31,32,33]. Our results provide possible explanations for these observations from a perspective of noiseinduced cell fate switching, and insights into how noise utilizes simple gene regulatory network to perform meaningful biological functions. Our unpublished work on the formation of the cortexdorsal midline border in the developing telencephalon suggests that the mutually inhibitory BMP and FGF signals can lead to the boundary refinement. This study provides a simple explanation on the purpose of such twogradient system in boundary sharpening.
Our results help elucidate the critical network structures for noise attenuation and pattern formation in development. Our method is general, and can be applied to other pattern formation or spatial relevant biological systems. With more biological regulatory data available, some more realistic morphogengene regulatory networks can be constructed. We anticipate that, by investigating these more realistic networks, more intricate mechanisms for the boundary formation and sharpening can be discovered. This will further our understanding of boundary sharpening, pattern formation, and other spatial related issues.
Conclusions
In this study, we discovered the probabilistic landscape for three motifs (MRSA, SA, and MR), interacting with morphogens, and investigated their boundary sharpening effects. The landscape results, along with the average switching time between attractors, provide a natural explanation for the boundary sharpening behavior depending on the noiseinduced gene state switching. The MRSA model displays more robust boundary sharpening ability against parameter perturbation, compared to the MR or SA model. This is supported by the results of switching time between attractors, because the MRSA model has shortest switching time among three models. In addition, introducing cross gradients of morphogens provides a more stable mechanism for the boundary sharpening of gene expression, due to a twoway switching mechanism. Our results reveal the critical network structures for noiseinduced boundary sharpening effects of gene expression, and promote the understanding for the critical network structures of spatial pattern formation in embryonic development.
Methods
From the topology of the network for three models, we construct dynamic models to describe the temporal evolution for the expression level of different genes in the network. We first construct an ordinary differentiation equation (ODE) model based on Hill cooperativity form representing activation or repression [22]. The ODE model include three terms: basal synthesis rates, activation or repression regulations from other genes, and selfdegradations. The three models (parameters can be found in Additional file 1: Tables S1S4) are represented by the following equations:

1.
The mutual repressed selfactivation model (MRSA)

2.
The selfactivation model (SA)

3.
Mutualrepression model (MR)

4.
The mutual repressed selfactivation model with cross morphogen gradients
Here, the ODE systems describe the temporal evolution of expression levels of X and Y genes. S represents the threshold of the sigmoidal function, and n is the Hill coefficient, which determines the steepness of the sigmoidal function [22]. The parameter A represents the ratio of the synthesis rate for the two genes, which measures the asymmetry level of the system. The parameter R represents the strength of the mutual inhibition between gene X and Y. In addition, a is the basal synthesis rate and k is the degradation rate for gene X and Y (see Additional file 1 for the descriptions of parameters, and Table S1 for the values of parameters). Take Eq. (1) as an example, the first term represents the regulation effect from gene X and Y, the second term represents the regulations from morphogen level M, and the last term represents the degradation of gene X or Y.
To generate the morphogen gradient with fluctuations, we used the model from [15], and made the spatial simulations (Eq. 5). Here [RA]_{ out } and [RA]_{ in } denote separately extracellular and intracellular RA concentrations. \( {\left[ RA\right]}_{out}\frac{{\mathit{\partial}}^2{W}_{out}\left(t,x\right)}{\mathit{\partial t\partial x}} \) and \( {\left[ RA\right]}_{in}\frac{{\mathit{\partial}}^2{W}_{out}\left(t,x\right)}{\mathit{\partial t\partial x}} \) represent standard white noise for extracellular and intracellular RA concentrations. RA_{ in } is the morphogen M we used in the simulations, and Ïµ quantifies the noise level in morphogenÂ M. For spatial simulations, no flux boundary conditions are applied along the y axis (mediallateral direction). Along the x axis, a noflux boundary condition is applied at the anterior margin and a leaky boundary condition is used at the posterior margin on the assumption that extracellular RA only leaves the embryo by diffusing through cells. We have varied Ïµ to see itâ€™s influence on the simulation results.
For the simulations from stochastic ODEs, the initial conditions are specified as (X0, Y0)â€‰=â€‰(0.5+ sN(0,1), 0.01), as we require the system starts from a location close to X attractor, i.e. X is much larger than Y. Here N is a random variable obeying standard norm distribution, and sâ€‰=â€‰0.05 is the magnitude for the fluctuations. For each model, we ran multiple simulations (100 times) to get the average results (Figs. 1, 4, 5 and 7).
Probabilistic landscape
In the cells, there exist intrinsic noise from statistical fluctuations of the finite number of molecules, and external noise from highly dynamical and inhomogeneous environments. Both of them can be significant to the dynamics of the system [11,12,13]. Therefore, one needs to study the cellular network dynamics in fluctuating conditions in order to model the cellular inner and outer environments realistically. The dynamics of a gene network in fluctuating environments can be addressed by: \( \dot{x}=\mathbf{F}(x)+d\ast \mathbf{x}\ast \zeta \), where xâ€‰=â€‰(x_{1}(t),â€‰x_{2}(t),â€‰â€¦,â€‰x_{ n }(t)) represents the vector of gene expression levels. F(x) is the vector for the driving force of gene regulations. Î¶ is Gaussian noise term satisfied with: <Î¶_{ i }(x,â€‰t)â€‰>â€‰â€‰=â€‰0 and <Î¶_{ i }(x,â€‰t)Î¶_{ j }(x,â€‰t^{â€²})â€‰>â€‰â€‰=â€‰2D_{ ij }Î´_{ ij }Î´(tâ€‰âˆ’â€‰t^{â€²}) (Î´_{ ij }â€‰=â€‰1 for iâ€‰=â€‰j, and Î´_{ ij }â€‰=â€‰0 for iâ€‰â‰ â€‰j), where Î´(t) is Dirac delta function and D is diffusion coefficient matrix. Here, d is the amplitude for white noise. So, the above ODEs (Eqs. 1, 2 and 3) can be transformed to stochastic ordinary differential equations (SODE).
The probability evolution for a stochastic dynamical system can be captured by the diffusion equations [34, 35]. For a 2dimensional system, the diffusion equation has the form:
Here, F_{1}and F_{2} represent the driving force for the system from above ODEs. D is the diffusion coefficient matrix.
By solving the diffusion equations for the long time limit, we obtained the steady state probabilistic distribution of the system. We used COMSOL Multiphysics (version 4.3) to solve the diffusion equations. In this way, we mapped out the potential landscape for the system by Uâ€‰=â€‰â€‰âˆ’â€‰ln(P_{ ss }) [17, 18, 22, 23, 36]. Here, P_{ ss } represents the probability distribution of the steady state, and U is the dimensionless potential.
We calculated the mean first passage time (MFPT) from the temporal trajectories of the variable X and Y. Starting from a random initial state at X attractor, following the temporal evolution of the system trajectory, we are able to find the time when the system first switches to the Y attractor. The time difference between initial time and the final time is defined as the first passage time (FPT) for the transition process from X attractor to Y attractor. Repeating this process we can obtain the average of the FPT, defined as the MFPT for the transition process from X attractor to Y attractor.
Abbreviations
 FPT:

First passage time
 MFPT:

Mean first passage time
 MR:

Mutualrepression
 MRSA:

Mutual repressed selfactivation
 SA:

Selfactivation
References
Lander AD. Morpheus unbound: reimagining the morphogen gradient. Cell. 2007;128(2):245â€“56.
Lander AD. How cells know where they are. Science. 2013;339(6122):923â€“7.
Ferguson EL, Anderson KV. Decapentaplegic acts as a morphogen to organize dorsalventral pattern in the Drosophila embryo. Cell. 1992;71(3):451â€“61.
Bouchoucha YX, Reingruber J, Labalette C, Wassef MA, Thierion E, DesmarquetTrin DC, Holcman D, GilardiHebenstreit P, Charnay P. Dissection of a Krox20 positive feedback loop driving cell fate choices in hindbrain patterning. Mol Syst Biol. 2013;9(1):690.
Xue Y, Zheng X, Huang L, Xu P, Ma Y, Min Z, Tao Q, Tao Y, Meng A. Organizerderived Bmp2 is required for the formation of a correct bmp activity gradient during embryonic development. Nat Commun. 2014;5(4):3766.
Calzolari S, Terriente J, Pujades C. Cell segregation in the vertebrate hindbrain relies on actomyosin cables located at the interhombomeric boundaries. EMBO J. 2014;33(7):686.
Xu Q, Mellitzer G, Robinson V, Wilkinson DG. In vivo cell sorting in complementary segmental domains mediated by Eph receptors and ephrins. Nature. 1999;399(6733):267.
Dormann D, Weijer CJ. Chemotactic cell movement during Dictyostelium development and gastrulation. Curr Opin Genet Dev. 2006;16(4):367â€“73.
Melen GJ, Levy S, Barkai N, Shilo BZ. Threshold responses to morphogen gradients by zeroorder ultrasensitivity. Mol Syst Biol. 2005;1(1)
Hasty J, Pradines J, Dolnik M, Collins JJ. Noisebased switches and amplifiers for gene expression. Proc Natl Acad Sci U S A. 2000;97(5):2075â€“80.
Swain PS, Elowitz MB, Siggia ED. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Natl Acad Sci U S A. 2002;99(20):12795â€“800.
Kaern M, Elston TC, Blake WJ, Collins JJ. Stochasticity in gene expression: from theories to phenotypes. Nat Rev Genet. 2005;6(6):451â€“64.
Thattai M, Oudenaarden AV. Intrinsic noise in gene regulatory networks. Proc Natl Acad Sci U S A. 2001;98(15):8614â€“9.
BalÃ¡zsi G, Oudenaarden AV, Collins JJ. Cellular decision making and biological noise: from microbes to mammals. Cell. 2011;144(6):910â€“25.
Zhang L, Radtke K, Zheng L, Cai AQ, Schilling TF, Nie Q. Noise drives sharpening of gene expression boundaries in the zebrafish hindbrain. Mol Syst Biol. 2012;8:613.
Tian H, Fung ES, Lei Z, Grace H, Monuki ES, Qing N. Semiadaptive response and noise attenuation in bone morphogenetic protein signalling. J R Soc Interface. 2015;12(107):20150258.
Wang J. Landscape and flux theory of nonequilibrium dynamical systems with application to biology. Adv Phys. 2015;64(1):1â€“137.
Li C, Wang J. Landscape and flux reveal a new global view and physical quantification of mammalian cell cycle. Proc Natl Acad Sci U S A. 2014;111(39):14130â€“5.
Li C. Identifying the optimal anticancer targets from the landscape of a cancerimmunity interaction network. Phys Chem Chem Phys. 2017;19(11):7642â€“51.
Waddington CH. The strategy of the genes A discussion of some aspects of theoretical biology. London: Allen and Unwin; 1957.
Wang J, Xu L, Wang E. Potential landscape and flux framework of nonequilibrium networks: robustness, dissipation, and coherence of biochemical oscillations. Proc Natl Acad Sci U S A. 2008;105(34):12271â€“6.
Li C, Wang J. Quantifying cell fate decisions for differentiation and reprogramming of a human Stem Cell Network: landscape and biological paths. PLoS Comput Biol. 2013;9(8):e1003165.
Li C, Wang J. Quantifying the landscape for development and Cancer from a Core Cancer stem cell circuit. Cancer Res. 2015;75(13):2607â€“18.
Li C, Wang J. Quantifying Waddington landscapes and paths of nonadiabatic cell fate decisions for differentiation, reprogramming and transdifferentiation. J R Soc Interface. 2013;10(89):20130787.
Zhang B, Wolynes PG. Stem cell differentiation as a manybody problem. Proc Natl Acad Sci U S A. 2014;111(28):10185â€“90.
Lu M, Jolly MK, Levine H, Onuchic JN, Benjacob E. MicroRNAbased regulation of epithelialhybridmesenchymal fate determination. Proc Natl Acad Sci U S A. 2013;110(45):18144â€“9.
Srinivasan S, Hu JS, Currle DS, Fung ES, Hayes WB, Lander AD, Monuki ES. A BMPFGF morphogen toggle switch drives the ultrasensitive expression of multiple genes in the developing forebrain. PLoS Comput Biol. 2014;10(2):e1003463.
Rogers KW, Schier AF. Morphogen gradients: from generation to interpretation. Annu Rev Cell Dev Biol. 2011;27(1):377â€“407.
Balaskas N, Ribeiro A, Panovska J, Dessaud E, Sasai N, Page KM, Briscoe J, Ribes V. Gene regulatory logic for reading the sonic hedgehog signaling gradient in the vertebrate neural tube. Cell. 2012;148(1):273â€“84.
Manu SS, Spirov AV, Gursky VV, Hilde J, AhRam K, Ovidiu R, VanarioAlonso CE, Sharp DH, Maria S. Canalization of gene expression and domain shifts in theDrosophilaBlastoderm by dynamical attractors. PLoS Comput Biol. 2009;5(3):e1000303.
Holloway DM, Lopes FJ, Da FCL, TravenÃ§olo BA, Golyandina N, Usevich K, Spirov AV. Gene expression noise in spatial patterning: hunchback promoter structure affects noise amplitude and distribution in Drosophila segmentation. PLoS Comput Biol. 2011;7(2):e1001069.
Sosnik J, Zheng L, Rackauckas CV, Digman M, Gratton E, Nie Q, Schilling TF. Noise modulation in retinoic acid signaling sharpens segmental boundaries of gene expression in the embryonic zebrafish hindbrain. Elife. 2016;5:e14034.
Li CJ, Hong T, Tung YT, Yen YP, Hsu HC, Lu YL, Chang M, Nie Q, Chen JA. MicroRNA filters Hox temporal transcription noise to confer boundary formation in the spinal cord. Nat Commun. 2017;8:14685.
Van Kampen NG: Stochastic processes in physics and chemistry, 3rd edn. North Holland, Amsterdam: Elsevier; 2007.
Hu G. Stochastic Forces and Nonlinear Systems. Shanghai: Shanghai Scientific and Technological Education Press; 1994. p. 68â€“74.
Zhang F, Xu L, Zhang K, Wang E, Wang J. The potential and flux landscape theory of evolution. J Chem Phys. 2012;137(6):065102.
Acknowledgements
C.L. thanks Dr. Qixuan Wang for helpful discussions.
Funding
C.L. is supported by the National Natural Science Foundation of China (11771098), Natural Science Foundation of Shanghai, China (17ZR1444500), and Thousand Youth Talents Program. L.Z. is partially supported by the National Natural Science Foundation of China No. 11622102, 91430217, 11421110001. Q.N. is partially supported by NIH grants R01GM107264, R01ED023050, and R01NS095355, NSF grant DMS1562176, and a grant by Jayne Koskinas Ted Giovanis Foundation for Health and Policy joint with Breast Cancer Research Foundation.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Author information
Authors and Affiliations
Contributions
CL, LZ and QN conceived and designed the study, analyzed the data. CL performed the study. CL, LZ and QN wrote the paper. All authors read and approved the final manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Competing interests
The author Qing Nie is currently acting as an Associate Editor for BMC Systems Biology.
Publisherâ€™s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1:
Figure S1. Comparisons of boundary sharpening effects over time quantified by sharpening index (SI), for MRSA models at different morphogen noise level Ïµ and gene expression noise level d value. Figure S2. Two dimensional simulations show the boundary sharpening effects over time at different vertical resolution (A for 10 grids and B for 20 grids). Blue: X is expressed, red: Y is expressed. Table S1. Parameters of the mutual repressed selfactivation (MRSA) model. Table S2. Parameters of the selfactivation (SA) model. Table S3. Parameters of the mutual repression (MR) model. Table S4. Parameters of the cross morphogen gradients model. (PDF 3832 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Li, C., Zhang, L. & Nie, Q. Landscape reveals critical network structures for sharpening gene expression boundaries. BMC Syst Biol 12, 67 (2018). https://doi.org/10.1186/s1291801805955
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1291801805955
Keywords
 Gene expression boundaries
 Morphogen
 Landscape
 Switching time