Implications of diffusion and time-varying morphogen gradients for the dynamic positioning and precision of bistable gene expression boundaries

The earliest models for how morphogen gradients guide embryonic patterning failed to account for experimental observations of temporal refinement in gene expression domains. Following theoretical and experimental work in this area, dynamic positional information has emerged as a conceptual framework to discuss how cells process spatiotemporal inputs into downstream patterns. Here, we show that diffusion determines the mathematical means by which bistable gene expression boundaries shift over time, and therefore how cells interpret positional information conferred from morphogen concentration. First, we introduce a metric for assessing reproducibility in boundary placement or precision in systems where gene products do not diffuse, but where morphogen concentrations are permitted to change in time. We show that the dynamics of the gradient affect the sensitivity of the final pattern to variation in initial conditions, with slower gradients reducing the sensitivity. Second, we allow gene products to diffuse and consider gene expression boundaries as propagating wavefronts with velocity modulated by local morphogen concentration. We harness this perspective to approximate a PDE model as an ODE that captures the position of the boundary in time, and demonstrate the approach with a preexisting model for Hunchback patterning in fruit fly embryos. We then propose a design that employs antiparallel morphogen gradients to achieve accurate boundary placement that is robust to scaling. Throughout our work we draw attention to tradeoffs among initial conditions, boundary positioning, and the relative timescales of network and gradient evolution. We conclude by suggesting that mathematical theory should serve to clarify not just our quantitative, but also our intuitive understanding of patterning processes.

at coordinates x and times t. We assume there is some set of coordinates where the system is instantaneously bistable, with boundaries denoted by x 1 (t), x 2 (t); that is, x < x 1 (t) implies (1) is instantaneously monostable (high U ) at time t, x 1 (t) ≤ x ≤ x 2 (t) implies (1) is instantaneously bistable at time t, and x 2 (t) < x implies (1) is instantaneously monostable (low U ). Assume several realizations of the process are given, each ultimately representing a horizontal slice of a 2D "embryo". For convenience we define X high (t) = {x : where the subscripts are derived with respect to the expected steadystate concentration of U .
We will consider the case where α(t, x) varies smoothly and monotonically in time, that is, t 2 > t 1 =⇒ α(t 2 , x) ≤ α(t 1 , x) ∀x, and asympototically approaches a finite value: α(t → ∞, x) → α * (x). For most of this discussion we will assume α(t, x) decreases, although it is conceptually straightforward to generalize our discussion to the case where α(t, x) increases. For our setup, decrease of α(t, x) in time means t 2 > t 1 =⇒ x 1 (t 2 ) < x 1 (t 1 ) and x 2 (t 2 ) < x 2 (t 1 ). We will further assume that X bi (∞) ⊂ X high (0) and X bi (0) ⊂ X low (∞), although this requirement may be relaxed.
If the temporal evolution of α(t, x) is much faster than the dynamics of u, then the behavior of the system is well approximated by (1) with α(x) = α * (x). If the dynamics of u(t, x) are much faster than the dynamics of α(t, x) and the system satisfies certain technical requirements including slow evolution of α(t, x) [1], then we expect that u(t, x) will rapidly converge to the instantaneous system attractors and will in effect traverse the stable paths of those attractors. Thus, x ∈ X high (0) =⇒ u(∞, x) high, and x ∈ X low (∞) =⇒ u(∞, x) low regardless of initial condition. In other words, at infinite time there should be no remaining irregularity at the boundary between high and low U .
If the dynamics are on similar timescales, however, then variation in initial conditions may exert influence on system outcome: Two trajectories u(t, x) with x ∈ X bi (∞) ⊂ X high (0) may converge to different steady states, such that the final boundary is imprecise. Guaranteeing a precise boundary, therefore, becomes a matter of restricting variation in initial conditions for x ∈ X bi (∞) to the basin of attraction of the stable path corresponding to high U .
To examine this question intuitively, we now assume a particular form of α(x, t) that allows us to rewrite (1) as a time-independent (autonomous) system. In particular, suppose that the morphogen at coordinate x decreases in concentration exponentially with rate r: where for simplicity of notation we have dropped the dependence of α on x. The time evolution of α(t) may be described by for initial condition α(0). We can then rewrite (1) as where α(0) ∈ R + =⇒ α(t) ∈ R + ∀t, b ∈ R, and we have again dropped the dependence of u and b on x for simplicity of notation. The system (2) has saddle-node bifurcations at each of . Since α(t → ∞) → 0, the steady states in this system are set by b: • one stable high state if b < − 2 . In a toy model of an embryo, we would expect r to remain constant across the axis while α(0) and b form spatial gradients determining the initial morphogen concentration and mono-/bistability respectively. Since monostable systems have outcomes that do not depend on initial conditions, to ensure boundary regularity we must ensure that bistable cells end up in the high state. Specifically, we seek to determine for which range of initial morphogen concentrations α(0) and chemical concentrations u(0) trajectories terminate in the high state. This is the definition of the basin of attraction for the high state, represented in Fig 1 by the region above the thick red line in each plot. Decreasing b or r increases the size of the basin of attraction for the high state, which matches the intuition that a longer time spent with a higher input ought to bias the system toward a higher state. As r increases, the system behaves effectively as though α(t) = 0, since trajectories evolve much more slowly than the change in the input level.

Convergence rate to steady state
The convergence rate to a stable steady state u * = (u * 1 , u * 2 ) at coordinate x is determined by the eigenvalues λ i of the Jacobian J(u * ) of Eq (3) in the main text for α(x) = be − x λ , as where denotes the real part of the argument. Then where the second minimum is taken across the set of stable steady states existing at coordinate x. For the model with a static gradient considered in Eq (3) in the main text, the eigenvalues are analytically determined to be  Figure 1: Slowing the rate of decrease in morphogen concentration causes more trajectories to reach the high steady state. Phase spaces (gray vector field) overlaid with sample trajectories (gold) for (2) with varying b, r. Darker trajectories start with higher α(0). Red circles denote stable stady states and the red x denotes the unstable steady state. The thick red line is the separatrix, so named because it separates initial conditions that end in the high state (above the line) from those that end in the low state (below the line). As b or r is decreased (to the top or left), the basin of attraction of the high state increases in size; i.e., trajectories starting with a wider range of initial conditions will end up in the high state. For fixed b, r, a boundary can form when the initial gradient α(0) spans the separatrix.