Intrinsic Noise Profoundly Alters the Dynamics and Steady State of Morphogen-Controlled Bistable Genetic Switches

During tissue development, patterns of gene expression determine the spatial arrangement of cell types. In many cases, gradients of secreted signalling molecules—morphogens—guide this process by controlling downstream transcriptional networks. A mechanism commonly used in these networks to convert the continuous information provided by the gradient into discrete transitions between adjacent cell types is the genetic toggle switch, composed of cross-repressing transcriptional determinants. Previous analyses have emphasised the steady state output of these mechanisms. Here, we explore the dynamics of the toggle switch and use exact numerical simulations of the kinetic reactions, the corresponding Chemical Langevin Equation, and Minimum Action Path theory to establish a framework for studying the effect of gene expression noise on patterning time and boundary position. This provides insight into the time scale, gene expression trajectories and directionality of stochastic switching events between cell states. Taking gene expression noise into account predicts that the final boundary position of a morphogen-induced toggle switch, although robust to changes in the details of the noise, is distinct from that of the deterministic system. Moreover, the dramatic increase in patterning time close to the boundary predicted from the deterministic case is substantially reduced. The resulting stochastic switching introduces differences in patterning time along the morphogen gradient that result in a patterning wave propagating away from the morphogen source with a velocity determined by the intrinsic noise. The wave sharpens and slows as it advances and may never reach steady state in a biologically relevant time. This could explain experimentally observed dynamics of pattern formation. Together the analysis reveals the importance of dynamical transients for understanding morphogen-driven transcriptional networks and indicates that gene expression noise can qualitatively alter developmental patterning.

where the bursty expression of each gene is introduced in a first approximation as a stoichiometry in the production reactions. More complete descriptions of the kinetic reactions would involve mRNA dynamics, promoter state dynamics or protein-protein interactions. The Chemical Langevin Equation (5) is the limiting description of this reaction scheme [1]. The CLE approximation will break down if the reaction rates change significantly within the time scale of the slowest reaction, and hence full simulation of the reaction kinetics (ST.1) is required [1]. Comparison between full simulations of (ST.1) and the solution of the CLE are presented in Fig. 5 showing the validity of the approximation for the biological scenario studied. Additionally, it is also interesting to point out that the regulatory functions described in (4) are in fact the binding probabilities for a haploid system. As long as the number of molecules of the relevant transcription factors and polymerase are not very low (assumptions already necessary to the derivation of (4)), in diploid cells, the kinetic equations are the same, subject to a doubling of the values of α i . The interpretation of the terms in (4) is slightly different since the production rate is proportional to 2 · Prob(both alleles of promoter bound)+ Prob(one allele of promoter bound). This statement also extends to the stochastic formulation where the same CLE is obtained.

Patterning without positive feedback loop M → A
To illustrate the effect on patterning time imposed by the regulatory architecture of the bistable switch, we compared it with a regulatory motif which lacks the cross repressive feedback between A and B and instead relies on strongly cooperative activation of A by M . In this non-feedback case, there is only one steady state and the expression level of A has a sigmoidal dependence on the level of morphogen, so A appears absent for low morphogen values, while A is expressed for high morphogen values.
In order to get a step-like sigmoidal response to the signal for a monostable system, it is necessary to use high cooperativity for the morphogen activation. In the context of the promoter state using the thermodynamical description this can be introduced as h independent binding sites for the activating signal on the promoter each one facilitating the binding of the RNAp when bound by the signal [2], Even though the feedback loop is removed, a certain constant amount of x B is considered to represent the boundary position at the same signal value. Alternatively this is equivalent to redefining the value of ρ A . For the comparison in the text, the values used are h = 25, x B = 3 and K B = 0.12 which give a sharp boundary that reaches its maximum expression at M 1 ( Fig. ST.1). The rest of the parameters are the same as the complete bistable switch to facilitate their comparison (see Fig. 2). We defined this model to generate an outcome that resembles the bistable switch case although quantitative comparison with the feedback case is limited by the specific parameters chosen. Nevertheless, the results suggest that, in the absence of feedback, patterning time is in general faster than the bistable switch for all values of the signal. Importantly, since the stable state never loses stability there is no dramatic increase in patterning time close to the boundary M A (Fig. ST.2). As a consequence, considering the intrinsic fluctuations of the system does not introduce any dramatic change in the patterning time ( Fig. ST.3).
More precisely, the deterministic non-feedback circuit dynamics are determined by the linear ODE, This predicts an exponential decay that for the initial condition n A = 0 is, The comparison in dynamical response of the bistable switch to the non-feedback motif can be extended to changes in parameters α and δ, the relative production and degradation rates respectively. Concretely, keeping the ratio α/δ fixed, the response speed of gene A over B can be controlled. In both cases, reducing the time scale response of B speeds up the whole transient and reduces the patterning time ( Fig. ST.2). In the case of the non-feedback model the reduction of patterning time is linear with respect to changes in α and δ. By contrast, the patterning time of the bistable switch is much more resilient to changes in α and δ (Fig. ST.2). This feature provides some robustness in patterning time to changes in parameters and means that the bistable switch is intrinsically slower than the non-feedback mechanism.

Numeric minimization of the Action
The minimisation of the action (9) was performed by discretising the continuous path ϕ τ in N straight segments between positions {φ 0 , φ 1 , φ 2 , ...φ N } of duration ∆t. With this approximation, the action (10) can be written as, where the index i runs over the two genes and k = 1, ..., N over each linear segment of the path. D i k and f i k are the ith component of the diagonal diffusion tensor (12) and the deterministic field (11) for each segment evaluated at its centre x k = (φ k−1 + φ k )/2. Additionally, the velocity along a segment isφ k = (φ k − φ k−1 )/∆t. For a certain fixed total time τ and fixed initial and final position ϕ 0 and ϕ N , the discretised action (ST.5) for a path can be minimised using a quasi-Newtonian method with the help of the analytical expression for the derivative of the action on the direction of the l-th where D j l,s and f j l,s are the derivative in the l direction of the jth component of the diffusion tensor (12) and the deterministic field (11) of segment s. The minimisation of the action with 2 × (N − 1) degrees of freedom was performed using the LFBGS minimisation strategy. The resulting value of the action S τ is the minimum action for a certain path time τ . The optimal action is found by minimising the resulting S τ with τ by means of a simple custom made line search method. To speed up the convergence to the solution the starting path consists of two straight lines joining the initial and final phenotype with the saddle point. The number of segments for the simulations used in the optimisation is N = 150, a larger number of segments did not exhibit an improvement in the precision of the action while increasing dramatically the computational time.
The actions for both switching directions S AB and S BA is determined by the boundary conditions of the integration fixing φ 0 and φ N at the initial and final expression amounts of the corresponding bistable states A and B. In order to obtain the expression value of the bistable states for different morphogen signals, the bistable switch was reconstructed using continuation strategies through the PyDSTool environment [3].

Simulation of stochastic trajectories
Simulations of the stochastic trajectories were performed using custom made code. Realisations of the Chemical Langevin Equation have been integrated using the Euler-Maruyama algorithm that reproduces the discrete Wiener process [4,5], while the exact realisations of the kinetic reaction equations were performed using the Gillespie algorithm [6]. Both methods require the reproduction of the whole trajectory and the computing time will be proportional to the trajectory time making these calculations very slow when computing stochastic switching events with increasing number of proteins Ω. By contrast, the computation of the action is immediate, independent of the actual transition time between states, making it a useful tool for studying spontaneous cell fate change.
The main drawback in computing the action is the unknown prefactor C (7). This means that the transition time can be computed only to logarithmic precision [7,8]. In cases in which a greater precision is required, the prefactor C can be obtained by fitting an exponential to average transition times obtained with the CLE or Gillespie at adequate values of Ω, such that Ω is large enough that the behaviour can be considered exponential, and low 3/4 enough to obtain reasonable computation times. This range can be guessed from the values of the action S while checking the goodness of the exponential fit. For the case of Fig. 5 the prefactor C was computed by choosing a range where the mean first passage time depended exponentially on Ω, and performing a least squares fitting to an exponential function T = Ce ΩS . The least squares fitting was not linearized to give more weight to points of large Ω, where the exponential behaviour is assured.