Patterning precision under nonlinear morphogen decay and molecular noise
Abstract
Morphogen gradients can instruct cells about their position in a patterned tissue. Nonlinear morphogen decay has been suggested to increase gradient precision by reducing the sensitivity to variability in the morphogen source. Here, we use cellbased simulations to quantitatively compare the positional error of gradients for linear and nonlinear morphogen decay. While we confirm that nonlinear decay reduces the positional error close to the source, the reduction is very small for physiological noise levels. Far from the source, the positional error is much larger for nonlinear decay in tissues that pose a flux barrier to the morphogen at the boundary. In light of this new data, a physiological role of morphogen decay dynamics in patterning precision appears unlikely.
Editor's evaluation
The authors use analytic calculations and numerical simulations to convincingly show that the purported benefits of nonlinear decay in morphogen gradients may be marginal in some cases and completely reversed in others (far from the concentration source). This is a valuable contribution to the field, as it questions common assumptions about the biological function of nonlinear morphogen decays during development.
https://doi.org/10.7554/eLife.84757.sa0Introduction
According to Wolpert’s famous French flag model (Wolpert, 1969), morphogen gradients encode readout positions ${x}_{\theta}$ via concentration thresholds ${C}_{\theta}=C({x}_{\theta})$, and differentiating cells base their fate decisions on whether the local morphogen concentration lies above or below such thresholds (Figure 1A). Thus, these readout positions mark the boundary locations between domains of different cell fates. Variations in the morphogen profile result in variations in the readout positions. The accuracy of the spatial information carried by morphogen gradients can be quantified with the positional error, which is defined as the standard deviation of the readout positions over different gradient realizations (Vetter and Iber, 2022):
How the observed precision of tissue patterns arising from this principle is achieved, in spite of natural molecular noise in morphogen production, transport, decay, internalization, turnover and other sources of variability, is a key question in developmental biology (Lander, 2011; Vetter and Iber, 2022; Iber and Vetter, 2022).
Morphogen dynamics are often described by reactiondiffusion equations of the form (Lander et al., 2009).
with morphogen concentration $C$, diffusion coefficient $D$, and decay rate $d$. $C}_{ref$ is a constant reference concentration that we introduce to make all units independent of $n$. The exponent $n$ models linear ($n=1$) or nonlinear ($n>1$) decay of the morphogen. Linear decay leads to exponential gradient profiles (Figure 1B) of the form (Lander et al., 2002).
with an amplitude C_{0} at the source at $x=0$, and a characteristic decay length $\lambda $, set by the diffusion coefficient $D$ and the degradation rate $d$. Nonlinear decay, on the other hand, results in shifted powerlaw gradients (Figure 1C; Eldar et al., 2003)
where ${\lambda}_{m}$ is a gradient length scale that depends on $D/d$, $n$ and ${C}_{0}/{C}_{\mathrm{ref}}$ (see Appendix 1). Nonlinear decay would for instance arise in case of cell lineage transport, when ligands interact with receptor clusters, or if ligand binding results in receptor upregulation, as is the case for several morphogens, most prominently for Hedgehog (Hh) (Eldar et al., 2003; Nahmad and Stathopoulos, 2009; Wartlick et al., 2009; Balaskas et al., 2012). Most reported morphogen gradient profiles have been fitted assuming linear decay ($n=1$) (Gregor et al., 2007; Gregor et al., 2008; Kicheva et al., 2007; Wartlick et al., 2011; Wartlick et al., 2014; Cohen et al., 2015; Zagorski et al., 2017; Mateus et al., 2020). For the FGF8 gradient in the developing mouse brain, $n\approx 4$ has been reported (Chan et al., 2017).
Embryos are subject to molecular noise, which can cause fluctuations in morphogen production and transport rates, and consequently, in the gradient amplitudes and morphogen fluxes from the source to the patterned cells. This results in shifts $\mathrm{\Delta}x$ between different gradient realisations (Figure 1B and C). In the case of linear decay, the shift is only related to the relative morphogen influx or amplitude and not to the absolute morphogen levels (Appendix 1, Equation 11). However, with nonlinear decay, the shift depends on the absolute levels, with higher influxes resulting in smaller shifts (Appendix 1, Equation 12). Previous research (Eldar et al., 2003) suggested that the circumstance that this shift vanishes for powerlaw gradients at sufficiently large morphogen influx values leads to more robust patterning, because the readout position becomes independent of the influx in this limit when the morphogen decay is nonlinear (Appendix 1). In other words, if two tissues are patterned by two different noisefree powerlaw gradients, both with high (but different) morphogen influxes from the source, the resulting gradients will be nearly identical, resulting in a reproducible pattern. For exponential gradients, the shift will not disappear, and the gradients will thus differ. However, the gradients that result from nonlinear decay also possess significantly shallower tails, relative to the higher concentration (Figure 1C). Their usefulness for patterning has therefore been questioned (Lander et al., 2009), and it has remained unclear whether nonlinearity in the morphogen decay would in fact help achieving higher positional accuracy. Indeed, to first order, the positional error of variable gradients is inversely proportional to the magnitude of their slope (Gregor et al., 2007; Vetter and Iber, 2022) according to
where ${\sigma}_{C}$ is the standard deviation of local morphogen concentration and $x$ denotes the patterning axis. This suggests that for patterning precision, the benefit of a smaller positional shift of gradients with $n>1$ might be offset or even overcompensated by their flatter shape further from the source.
In the previous analysis (Eldar et al., 2003), molecular noise was considered only in the form of a foldchange in the morphogen amplitude or influx between different gradient realisations, resulting in shifted deterministic gradients, as shown in Figure 1B and C. To account for the intrinsic stochasticity of biological systems, we now extend this deterministic view by incorporating kinetic variability into the model, as depicted in Figure 1D and E. This is achieved by introducing randomness into all kinetic parameters of the reactiondiffusion equation. Our model is cellbased, meaning that each cell in the tissue is assigned its own specific variable kinetic parameters, emulating intercellular variability (Methods). With this quantitative statistical tool, we demonstrate numerically that the positional error of noisy morphogen gradients does not significantly improve with nonlinear decay. In the contrary, if the morphogen cannot leave the patterned tissue opposite of the source, the powerlaw gradients become shallow in a substantial part of the domain, leading to reduced positional accuracy with nonlinear decay.
Results
Noisy gradient model
We simulated steadystate diffusion to study the impact of nonlinear decay on the precision of noisy morphogen gradients. Our model uses a onedimensional cellular domain composed of a source of length ${L}_{\mathrm{s}}$ and a patterning region of length ${L}_{\mathrm{p}}$ (Figure 1D and E). To represent morphogensecreting source cells explicitly, the diffusion equation (Equation 2) was extended by a morphogen production term, resulting in
Here, $H(x)$ is the Heaviside function, ensuring that production at rate $p$ only occurs in the source ($x<0$). Zeroflux boundary conditions were used at both outer ends of the tissue, mimicking a situation in which morphogen molecules are restricted to the patterned tissue by an impermeable boundary:
We generated variable morphogen gradients by numerically solving Equation 6 with kinetic parameters p_{i}, d_{i} and ${D}_{i}$, and cell areas ${A}_{i}$ independently drawn from lognormal distributions for each cell $i=1,\mathrm{\dots},N$ in the domain (Vetter and Iber, 2022; Adelmann et al., 2022; Figure 2, for details see Methods). The individual gradient realisations ${C}_{j}(x)$ can be thought of as representing different embryos, denoted by the index $j$. They exhibit inter and intratissue variability due to the stochastic nature of the three kinetic parameters that vary from cell to cell. Cells have to convert the spatial morphogen distribution they are exposed to into a single concentration value, which determines their fate in the tissue according to the French flag model. There are several ways cells may achieve this, such as averaging the morphogen signal over their entire cell surface, beyond their cell surface via a cilium, or reading out the signal at a single point. In a recent study, we found that the readout mechanism has little impact on the gradient precision perceived by the cells (Adelmann et al., 2022). We therefore only analysed the case where cells average the morphogen signal over their cell surface, or over their diameter in the 1D model here, respectively. The concentration in each cell is then compared to the threshold concentration ${C}_{\theta}$ and the location ${x}_{\theta ,j}$ of the first cell whose sensed concentration subceeds this threshold is recorded. This process is repeated for all gradients $j$, allowing to compute the positional error according to its definition (Equation 1), ${\sigma}_{x}={\mathrm{stddev}}_{j}\{{x}_{\theta ,j}\}$ to quantify the precision of the positional information conveyed by the morphogen gradients.
Model parameters
To define the stochastic nature of the morphogen kinetics involved in the formation of the gradients, we express the mean value and standard deviation of a parameter $q$ by ${\mu}_{q}$ and ${\sigma}_{q}$, respectively (Figure 2). Based on measurements of the Hedgehog morphogen gradient in the Drosophila wing disc and mouse neural tube (Kicheva et al., 2007; Cohen et al., 2015; Vetter and Iber, 2022), we used a mean diffusivity of ${\mu}_{D}=0.033$ µm^{2}/s and a mean gradient length ${\mu}_{\lambda}=20$ µm. We furthermore set the average degradation rate to ${\mu}_{d}={\mu}_{D}/{\mu}_{\lambda}^{2}$, and the average production rate to ${\mu}_{p}={\mu}_{d}{C}_{\mathrm{ref}}$, where ${C}_{\mathrm{ref}}=1$ arb. units to normalise the concentrations. Other specific values of gradient parameters would not change the results reported here, which are for the steady state, but would only alter the timescale it takes for the steady state to be reached. The noisetosignal ratio in each quantity $q$ is given by the corresponding coefficient of variation, ${\mathrm{CV}}_{q}={\sigma}_{q}/{\mu}_{q}$. Reported physiological noise levels in morphogen production, decay, and transport differ between morphogens and tissues, but are around ${\mathrm{CV}}_{p,d,D}\approx 0.3$ (Vetter and Iber, 2022), which we use to define the distribution widths of the kinetic parameters.
In addition to the morphogen kinetics, our simulations also include celltocell variability in the cell areas. The widths and crosssectional areas of cells vary in all layers along the apicalbasal axis (Gómez et al., 2021). Most quantifications have been carried out on the apical surface. One of the highest reported values for the apical area CV is found in the vertebrate neural tube (${\mathrm{CV}}_{A}\approx 0.9$) (Escudero et al., 2011; Guerrero et al., 2019; BocanegraMoreno et al., 2022), but most values are considerably lower (Kokic et al., 2019). We therefore used ${\mathrm{CV}}_{A}=0.5$ in all simulations unless specified otherwise.
The diffusion coefficient, $D$, and the degradation rate, $d$, set the steadystate patterning length scale, $\lambda =\sqrt{D/d}$. Thus, our results are independent of the absolute values chosen for $D$ and $d$, and only depend on their ratio. Positional quantities, such as the positional error, are reported relative to the average cell diameter, which in turn was chosen to be a fixed multiple of the average gradient decay length. We fixed the average cell diameter at a fourth of the exponential gradient length, ${\mu}_{\delta}/{\mu}_{\lambda}=1/4$, as found in the developing mouse neural tube (Kicheva et al., 2014; Cohen et al., 2015).
Impact of nonlinear decay on gradient precision
We previously showed that in case of linear decay, there is a negligible impact of cell area variability as long as ${\mathrm{C}\mathrm{V}}_{A}<1$ (Adelmann et al., 2022). We now find that this holds similarly for nonlinear decay (Figure 3A), justifying the use of a fixed ${\mathrm{CV}}_{A}=0.5$ in the remainder of our analysis.
Much as for linear decay (Adelmann et al., 2022), the positional error scales with the square root of the mean cell diameter also for nonlinear decay (Figure 3B). Small cell diameters, as observed in all known tissues that employ gradientbased patterning (Adelmann et al., 2022), are therefore important for high spatial precision also in case of nonlinear decay.
The positional error increases from less than one cell diameter close to the source to about two cell diameters at a distance of 75 cell diameters away from the source (Figure 3C). Close to the distant domain boundary opposite of the source, where a noflux condition was imposed, the positional error rapidly increases for nonlinear decay, while remaining relatively low for linear decay. If only the production in the source is varied (${\mathrm{CV}}_{p}=0.3$, ${\mathrm{CV}}_{d,D}=0$), the positional error remains constant as the readout distance from the source increases, but increases again sharply close to the distant end in case of nonlinear decay (Figure 3D). But even for strong nonlinearity ($n=4$), the positional error remains in the subcellular range when only production noise is considered, as long as the readout position is further than about $\lambda $ away from the distal end.
Independent of whether all parameters are varied or only the production rate, the positional error drops in close vicinity to the source with stronger nonlinearity in the decay (insets of Figure 3C and D). However, with less than 20% of a single cell diameter from $n=1$ to $n=4$, the effect is likely too small to be physiologically relevant. Further away from the source, linear decay yields a smaller positional error than nonlinear decay (Figure 3D–E). No matter how long the patterning domain is, nonlinearity always increases the positional error as the distal tissue boundary is approached (Figure 3F).
What then causes the increased positional errors with nonlinear decay near the distal domain boundary? A zeroflux boundary condition there implies shallower gradients than on infinite domains: ${C}^{\prime}(x)\to 0$ as $x\to {L}_{\mathrm{p}}$. This effect occurs irrespective of $n$, but the spatial range over which the gradient flattens (and thus deviates from the pure exponential and shifted powerlaw forms for infinite domains, Equations 3; 4) increases with $n$. By virtue of Equation 5, nonlinear decay thus leads to greater positional errors at readout positions in the vicinity of the distal boundary compared to linear decay.
In summary, our computer simulations of noisy morphogen gradients suggest that it is insufficient to quantify gradient robustness and patterning precision by considering variability in the morphogen production alone. Moreover, if the morphogen cannot exit the patterning domain opposite of the source, shifted powerlaw gradients that result from nonlinear morphogen decay flatten over a significantly larger range than exponential gradients, leading to increased positional errors. The gain in positional accuracy close to the source for nonlinear decay is negligible and therefore barely physiologically relevant. Overall, exponential gradients lead to more robust patterning.
Impact of boundary condition at the source
Given the impact of the distal domain boundary, we wondered whether the representation of the morphogen source by either a spatial production domain (Figure 4A), by a flux boundary condition $D{C}^{\prime}(0)={j}_{0}$ (Figure 4B) as used by Eldar et al., 2003, or by a fixed gradient amplitude $C(0)={C}_{0}$ (Figure 4C) would affect the positional error predicted by the model. While there are small quantitative differences, the gradient shapes (Figure 4A–C) and positional errors (Figure 4D–F) are overall very similar.
As we increase the variability in the production rate via ${\mathrm{CV}}_{p}$ (Figure 4D), in the influx from the source via ${\mathrm{CV}}_{{j}_{0}}$ (Figure 4E), or in the gradient amplitude at the source boundary via ${\mathrm{CV}}_{{C}_{0}}$ (Figure 4F), we find the smallest increase in the positional error for the production rate and the largest increase for the gradient amplitude. Neumann or Dirichlet boundary conditions thus overestimate the positional error when the variability in the source is high. Instead of using such boundary conditions, a spatial source domain should explicitly be modeled, where applicable. With the physiological values ${\mathrm{CV}}_{p}\approx 0.3$ and ${\mathrm{C}\mathrm{V}}_{{C}_{0}}\u2a850.3$ (Vetter and Iber, 2022), however, variability in the morphogen production plays merely a subordinate to moderate role in the overall gradient variability. Molecular noise in morphogen degradation and diffusivity dominates the patterning precision.
Impact of the morphogen source strength
As the gradient amplitude determines the sensitivity of the readout position to amplitude changes for nonlinear decay (Appendix 1, Equation 12) but not for linear decay (Appendix 1, Equation 11), the average morphogen production rate is expected to affect the patterning accuracy in the case of nonlinear decay, but not for linear decay. We put this theoretical prediction to the test by varying the mean relative production rate, the mean influx from the source, and the mean morphogen amplitude in the three different simulated morphogen production models. Changes in these parameters have no effect on the positional error if morphogen degradation is linear, which is consistent with the theory (Figure 5A–F, blue lines). With nonlinear decay, on the other hand, we indeed observe the positional error to be highly dependent on morphogen abundance. Precision arguments previously brought forward for deterministic morphogen gradients (Eldar et al., 2003) do not appear to directly quantitatively translate to the positional error in settings where celltocell variability is included, and where morphogen production remains at physiological levels.
For high morphogen supply levels, nonlinear decay leads to a smaller positional error close to the source (Figure 5A–C). The effect is, however, substantially less pronounced in the model that includes a spatial morphogen source domain (Figure 5A) than in those that do not (Figure 5B and C), highlighting once again the limitations of the latter. With an explicit source domain, nonlinear decay yields only marginally more spatial accuracy, when production is high ($p/d{C}_{\mathrm{r}\mathrm{e}\mathrm{f}}\u2a860.4$). Lower production levels increase the positional error close to the source substantially in all three models, reaching several cell diameters, for $n>1$. The gradients effectively flatten out at low production, reducing their usefulness for spatial tissue patterning. The stronger the nonlinearity in the degradation, the more pronounced this loss of patterning precision.
Further away from the source, the benefit of nonlinear decay is lost entirely, and exponential gradients remain more precise than shifted powerlaw gradients also at high morphogen supply levels (Figure 5D–F).
In summary, simplified models without explicit representation of morphogensecreting cells overestimate the beneficial impact of nonlinear decay on patterning precision. In all models considered here, the benefit of nonlinear morphogen decay is restricted to a close vicinity of the morphogen source, where patterning precision is high anyway (Vetter and Iber, 2022) and may thus not be as critical for robust development, and to a regime of very strong morphogen production. Further into the tissue, and at moderate morphogen abundance, linear decay yields more accurate patterning.
Discussion
Nonlinear morphogen decay was proposed as a potential precisionenhancing mechanism for tissue patterning in the seminal theoretical work by Eldar et al., 2003 in a deterministic setting, where morphogen gradients are devoid of noise. Here we have explored this idea with a stochastic model, taking noisy gradients into account, as they arise from celltocell variability in morphogen kinetics. The surprising outcome of our quantitative analysis is that, while a small advantageous effect indeed exists near the morphogen source, this gain is outweighed by a substantial loss of precision in the spatial information that signalling gradients provide to cells in the interior and distal parts of a patterned tissue when morphogen decay is nonlinear. In tissues that pose a diffusion barrier to the signalling molecule at their boundary, shifted powerlaw gradients that emerge with selfenhanced degradation, flatten out over a substantial portion of the spatial domain, whereas exponential gradients remain more graded (Figure 1). This leads to greater spatial precision with linear decay (Figure 3), and is contrary to the original expectation (Eldar et al., 2003).
This longrange boundary effect is not the only reason why linear morphogen decay is favourable for precise pattern formation. The positional error, which is the decisive quantity that measures the spatial accuracy with which cells can determine their location in the pattern, and ultimately their fate in differentiation, is highly sensitive to morphogen supply levels when morphogen decay is nonlinear, but largely insensitive when decay is linear (Figure 5). This implies that patterning is more robust to variations in the size and strength of the morphogensecreting source, if decay is linear. These results challenge the established view that powerlaw gradients buffer fluctuations in morphogen production (Eldar et al., 2003). We find that the positional error behaves in the opposite way, buffering production fluctuations only with linear, but not with nonlinear decay. From an evolutionary perspective, the linear case may be favoured, as patterning precision is unaffected by changes in the size and kinetics of the morphogensecreting source only if $n=1$.
Our study demonstrates that a stochastic approach is required to quantify patterning precision of real noisy gradients. Moreover, we find the positional error to be overestimated in simplified models that replace the morphogensecreting cells by a Neumann or Dirichlet boundary condition (Figure 4). Based on this, we recommend to include an explicit representation of the source in future theoretical or numerical work on the subject, as we did with Equation 6.
Distinguishing exponential gradients from shifted power laws can be very difficult in practice, as they can appear similar over the short distances over which they can be reliably measured with classical imaging techniques. The FGF8 gradient in the developing mouse brain is the only case we are aware of where $n>1$ has been reported robustly (Chan et al., 2017), and whether this is linked to patterning precision in any way remains unclear. Available gradient data in other systems, such as Sonic Hedgehog and Bone Morphogenetic Protein in the neural tube (Zagorski et al., 2017), is too variable to confidently reject the hypothesis that $n=1$. Most further reports of morphogen gradient shapes (Gregor et al., 2007; Gregor et al., 2008; Kicheva et al., 2007; Wartlick et al., 2011; Wartlick et al., 2014; Cohen et al., 2015; Zagorski et al., 2017; Mateus et al., 2020) are consistent with exponentials within measurement accuracy. New measurement techniques are needed to determine whether nonlinear decay is at work in the formation of known morphogen gradients during development. In light of our findings, a physiological role of nonlinear ligand decay in patterning precision appears implausible. If anything, our data suggest an overall advantage of linear decay, also considering the evolutionary aspect of tissue size and protein synthesis rate differences between species.
The morphogen concentration declines significantly over several orders of magnitude, independent of whether there is linear or nonlinear decay. At low morphogen concentrations, thermal fluctuations and stochastic binding kinetics of ligands and receptors will affect gradient and readout precision (Berg and Purcell, 1977; Lauffenburger and Linderman, 1996; Lander et al., 2009). Cells can, in principle, achieve high readout precision despite such fluctuations via spatial and temporal averaging (Lauffenburger and Linderman, 1996). To assess such effects, quantitative measurements of morphogen numbers and cellular responses would be required far away from the source. This requires the further development of more sensitive measurement technology (Lelek et al., 2021). Once the absolute concentration levels of the morphogen gradients can be determined, it can be assessed whether the approximation of the gradients by a continuous functions is valid along the whole tissue or, whether discrete models have to be considered.
In future work, the simulated gradients can be used as inputs to complex downstream networks, and the effect of noise in the readout can be studied. However, these downstream networks would not alter the relative precision of gradients generated by linear and nonlinear decay. In conclusion, nonlinear decay may slightly enhance precision close to the source, but it rapidly deteriorates far from the source.
Methods
Generation of variable morphogen gradients
To generate noisy morphogen gradients numerically, we constructed the onedimensional cellular domains in an iterative process, cell by cell. For each cell $i$, an area ${A}_{i}$ was drawn from a lognormal distribution with specified mean value ${\mu}_{A}$ and coefficient of variation ${\mathrm{CV}}_{A}$ (Adelmann et al., 2022). The drawn area was then converted to a cell diameter ${\delta}_{i}=2\sqrt{{A}_{i}}/\pi $. Using the transformation properties of lognormal distributions, the cell areas was drawn according to
allowing to accurately fix the mean cell diameter ${\mu}_{\delta}$. This procedure was repeated for cells $i=1,2,3\mathrm{\dots}$ until the sum of the diameters equaled the source length ${L}_{\mathrm{s}}$ or the patterning domain length ${L}_{\mathrm{p}}$. The spatial axis was then discretized into cellular subintervals accordingly (Figure 2). We used a patterning domain length of 200 cells (${L}_{\mathrm{p}}=200{\mu}_{\delta}$) and a source domain length of 5 cells (${L}_{\mathrm{s}}=5{\mu}_{\delta}$), unless otherwise stated.
Once the patterning axis was constructed, the kinetic parameters ${p}_{i},{d}_{i},{D}_{i}$, were drawn from lognormal distributions for each cell $i$ independently. For simulations without explicit source domain, a random morphogen influx j_{0} or an amplitude C_{0} was also drawn from a lognormal distribution. We then numerically solved Equation 6 using Matlab’s builtin fourthorder boundary value problem solver bvp4c (version R2020b). At cell boundaries, we imposed continuity of both the morphogen concentration and flux. Repeating this procedure 10^{3} times using independent random parameters and cell areas yielded statistically independent realisations of noisy morphogen gradients. To estimate the standard errors of the positional errors as shown in the plots, we used bootstrapping.
Choice of parameter distribution
In this article, we assume lognormally distributed cell areas and kinetic parameters, analogous to our previous works (Vetter and Iber, 2022; Adelmann et al., 2022). For the cell areas, this choice is rooted in the reported distributions of apical areas in the Drosophila larval and prepupal wing discs, and in the mouse neural tube (SánchezGutiérrez et al., 2016; Guerrero et al., 2019). The results reported here are, however, largely independent of the probability distribution, as long as it satisfies certain physiological criteria:
The random parameters must be strictly positive. This rules out probability distributions which allow for negative values, including for example a normal distribution.
The probability of drawing a nearzero parameter must vanish quickly. This is because tiny diffusion coefficients, fluxes, or amplitudes do not allow for successful patterning over biologically relevant distances or timescales. A normal distribution truncated at zero, for example, is ruled out because minuscule diffusion coefficients would occur frequently.
In recent related work (Adelmann et al., 2022), we demonstrated that other distributions which fulfill the above criteria yield similar results.
If the morphogen source is not modeled explicitly (omitting the production term in Equation 6), the gradient amplitude or morphogen influx levels at the source boundary can serve as a proxy for the production of the morphogen. For these simulations, the amplitudes and fluxes were also drawn from lognormal distributions. The width of these distributions is controlled via their coefficients of variation, ${\mathrm{CV}}_{{\mathrm{C}}_{0}}$ or ${\mathrm{CV}}_{{j}_{0}}$ as specified in the respective figures.
Appendix 1
Qualitative difference between linear and nonlinear morphogen decay
In this section we present some theoretical considerations about the consequences of nonlinear decay for noisefree morphogen gradients. We consider deterministic steadystate gradients obtained by analytically solving
Solving Equation 7 on an infinite onedimensional domain for linear morphogen decay ($n=1$) with a concentration that drops to zero at infinite distance from the source ($C(x)\to 0$ as $x\to \mathrm{\infty}$), results in exponential gradient profiles (Figure 1A) of the form
with an amplitude C_{0} at the source at $x=0$. The amplitude can be set by Dirichlet boundary conditions or by flux boundary conditions at the source, ${D\partial C/\partial x}_{x=0}={j}_{0}.$ Imposing flux boundary conditions leads to an amplitude ${C}_{0}={j}_{0}\lambda /D$. Thus, with linear decay, influx and amplitude are proportional.
Nonlinear decay (Equation 7, $n>1$), results in shifted powerlaw gradients (Figure 1A) that can be expressed as
where
is a length scale determining the shift in the power law, and ${C}_{0}=C(0)$ is the amplitude analogous to Equation 8. As the linear decay is approached ($n\to 1$), $m$ diverges ($m\to \mathrm{\infty}$), the powerlaw length scale approaches the exponential length scale (${\lambda}_{m}\to \lambda $), and the powerlaw gradients (Equation 9) become exponential (Equation 8). For a flux boundary condition at the source, the morphogen amplitude is
Amplitude and influx at the source boundary are thus not proportional for nonlinear morphogen decay, unlike in the linear case. Moreover, powerlaw gradients do not possess a constant gradient decay length $\lambda $ that quantifies a distance over which a foldchange in morphogen concentration occurs. Nevertheless, if one were to locally fit an exponential to the powerlaw gradient (Eldar et al., 2003),
one would observe the “gradient decay length” $\lambda (x)$ to increase with the distance from the source according to $\lambda (x)=x/(m\mathrm{ln}x)$ (Figure 1C).
Morphogen gradients define readout positions ${x}_{\theta}$ via concentration thresholds ${C}_{\theta}=C({x}_{\theta})$ (Figure 1A, D and E). For linear decay, the readout position follows from Equation 8 as
and for nonlinear decay from Equation 9 as
Due to inevitable molecular noise in morphogen production, transport, and decay, morphogen gradients differ between embryos, and hence readout positions ${x}_{\theta ,i}$ vary between different gradient realisations $i$ for both linear and nonlinear morphogen decay (Vetter and Iber, 2022). In the past, the impact of changes in morphogen production on readout precision has been studied for gradients that remain otherwise unchanged between embryos (Eldar et al., 2003). We now revisit this scenario. In response to a change in the morphogen amplitude from C_{0} to ${C}_{0}^{*}$, the readout position shifts along the patterning axis (Figure 1B and C). For linear decay, this shift $\mathrm{\Delta}x$ is independent of the absolute gradient amplitude C_{0} and depends only on the relative amplitude change, ${C}_{0}^{*}/{C}_{0}$, and the characteristic gradient length $\lambda $:
For nonlinear decay, the shift is given by
According to Equation 12, the shift $\mathrm{\Delta}x$ is proportional to ${\lambda}_{m}$ which in turn is proportional to ${C}_{0}^{1/m}$, implying that the shift increases with decreasing amplitude (Appendix 1, Figure 1A,B). This dependency of nonlinear decay on the gradient amplitude qualitatively distinguishes linear from nonlinear decay. Alternatively, the shift may be expressed as a function of a change in morphogen influx from the source from j_{0} to ${j}_{0}^{*}$. For linear decay, it simply reads
because flux and amplitude are proportional, making the shift again independent of absolute morphogen levels. For nonlinear decay, however, amplitude and influx are related as
The resulting readout shift is therefore
with a powerlaw length scale ${\lambda}_{m}$ that can be expressed in terms of the influx j_{0} as
Therefore, since $\mathrm{\Delta}x$ is proportional to ${\lambda}_{m}$, which is in turn proportional to ${j}_{0}^{1/(m+1)}$, the shift also increases with decreasing influx (Appendix 1 Figure 1A,C), albeit slower than with the amplitude.
Data availability
The current manuscript is a computational study, so no data has been generated for this manuscript. The source code is released under the 3clause BSD license. It is available as a public git repository at https://git.bsse.ethz.ch/iber/Publications/2022_adelmann_vetter_nonlinear_decay (copy archived at Adelmann, 2023).
References

Software2022_adelmann_vetter_nonlinear_decay, version swh:1:rev:9dd42cc39087ce2edcfccd4af67a8ea40f471194Software Heritage.

Physics of chemoreceptionBiophysical Journal 20:193–219.https://doi.org/10.1016/S00063495(77)855446

Ptch1 and Gli regulate Shh signalling dynamics via multiple mechanismsNature Communications 6:6709.https://doi.org/10.1038/ncomms7709

Epithelial organisation revealed by a network of cellular contactsNature Communications 2:526.https://doi.org/10.1038/ncomms1536

Relationship between epithelial organization and morphogen interpretationCurrent Opinion in Genetics & Development 75:101916.https://doi.org/10.1016/j.gde.2022.101916

Do morphogen gradients arise by diffusion?Developmental Cell 2:785–796.https://doi.org/10.1016/s15345807(02)00179x

The measure of success: constraints, objectives, and tradeoffs in morphogenmediated patterningCold Spring Harbor Perspectives in Biology 1:a002022.https://doi.org/10.1101/cshperspect.a002022

BookReceptors: Models for Binding, Trafficking, and SignalingOxford University Press.

SingleMolecule localization microscopyNature Reviews. Methods Primers 1:39.https://doi.org/10.1038/s4358602100038x

Precision of morphogen gradients in neural tube developmentNature Communications 13:1145.https://doi.org/10.1038/s41467022288343

Morphogen gradient formationCold Spring Harbor Perspectives in Biology 1:a001255.https://doi.org/10.1101/cshperspect.a001255

Positional information and the spatial pattern of cellular differentiationJournal of Theoretical Biology 25:1–47.https://doi.org/10.1016/s00225193(69)800160
Decision letter

Gordon J BermanReviewing Editor; Emory University, United States

Aleksandra M WalczakSenior Editor; CNRS, France

Gordon J BermanReviewer; Emory University, United States
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Patterning precision under nonlinear morphogen decay and molecular noise" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Gordon J Berman as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Aleksandra Walczak as the Senior Editor.
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
Although the reviewers all agreed that the study was rigorously performed and reaches an interesting and valuable conclusion, the paper is difficult read for people unfamiliar with the current literature. The reviewers ask that the authors clarify the text of the article (see reviews below), allowing a reader uninformed in the field of morphogen gradient calculations (e.g., a biophysicist working in another field) to have an easier time understanding the key findings here.
Reviewer #1 (Recommendations for the authors):
These morphogens often serve as transcription factors to help set up other, also noisy, downstream morphogen gradients. Thus, it would be interesting to have the authors' comment on whether having a linear + a nonlinear or a nonlinear + a nonlinear set of decay functionals might change the readout accuracy further away from the source. As these calculations might be timeconsuming, I'm not specifically asking for them to be done, but it would be good to get the authors' take on the matter for the Discussion section.
My main concern with the submission is that the article is written in a very short format that is very difficult for someone outside the field (myself included) to read. All of the information was there and I was able to figure things out, but only with significant effort. I would recommend expanding the text somewhat and providing a bit more detail as to where previous mathematical approaches ended and where this one starts.
Reviewer #2 (Recommendations for the authors):
1. This is a bit of author preference, but I find that I prefer a freestanding methods section as opposed to having them interleaved into the results, as they are now. I think the freestanding section allows for a bit more elaboration on details without detracting from the narrative.
2. I generally found this paper a nice exploration of the question of how nonlinear decay affects patterning. The paper is written in such a way that it may be challenging to read for those less familiar with the equations and conventions of mathematical modeling. If the authors wish to make the paper more accessible, I would suggest (in addition to the separation of methods) spending a bit more time defining the variables and parameters used throughout the paper, the justification for the different types of x/y axes used in the different graphs, and perhaps adding a cartoon or two illustrating a hypothetical tissue being patterned by a morphogen with some of the parameters/variables depicted. In addition, giving specific biological examples that correspond to some of the model considerations (e.g. boundary conditions) may be interesting for some readers.
3. It might be interesting to add two points to the Discussion section. First, the simulations here explore the effect of parameter variation, which might be a good representation of, e.g. how genetic variability between individuals affects parameters and therefore patterning. However, even with consistent mean parameter values, molecular stochasticity can also impact these processes – is this worth exploring in future work? Second, given that most (all?) morphogens seem to operate in more complex regulatory networks, with features like opposing gradients or mutual inhibition of morphogen targets, might the interaction between nonlinear decay and these other features be worth exploring in the future?
Reviewer #3 (Recommendations for the authors):
I find that this paper is very hard to read. The results of this paper are interesting, but, as written, this paper will appeal to a very narrow audience of researchers who have worked or work on modeling the production/diffusion/degradation model of morphogen gradients. This stems from two reasons. First, there is a lack of conceptual explanation of what the relevant questions really are. The authors refer to a previous paper but do not explain what the arguments in that paper really were. As far as I understand, the previous paper claimed that the nonlinear degradation will buffer against fluctuations in the levels of the source. That is interesting in principle, but there are other sources of noise and all of it is hard to get. I would urge the authors to rewrite their paper putting themselves in the shoes of someone who is not deeply familiar with all the details of the literature. The second reason why I find this paper confusing is that it is in part difficult to understand what the authors did. I had to check a previous paper to find some explanation of how molecular noise was introduced in the simulations. However, even after quickly checking that paper I remain confused about what the real assumptions were and how much I can trust the choice of parameters, etc. Again I would strongly recommend that the authors write a paper that is easier to read and contains all the relevant information.
https://doi.org/10.7554/eLife.84757.sa1Author response
Reviewer #1 (Recommendations for the authors):
These morphogens often serve as transcription factors to help set up other, also noisy, downstream morphogen gradients. Thus, it would be interesting to have the authors' comment on whether having a linear + a nonlinear or a nonlinear + a nonlinear set of decay functionals might change the readout accuracy further away from the source. As these calculations might be timeconsuming, I'm not specifically asking for them to be done, but it would be good to get the authors' take on the matter for the Discussion section.
Unfortunately, it is unclear to us which patterning scenario the referee is referring to specifically. In a scenario where a morphogen gradient with nonlinear decay establishes a second gradient with nonlinear decay further away from the source, for example Hh instructing Dpp expression in the Drosophila wing disc [Tonimato et al. 2000, Molecular Cell], the first gradient would already lead to an imprecise readout, and thus an imprecise setup of the downstream gradient. Another possibility would be opposing gradients, as observed in the neural tube. However, we previously demonstrated that opposing gradients are not required for high patterning precision in the mouse neural tube [Vetter et al. 2022, Nat. Commun.]. Furthermore, a nonlinear decay gradient could pattern the tissue close to the source and a second linear decay gradient could pattern the tissue further in the domain, but we are not aware of a biological example.
As we only see a minuscule improvement of nonlinear decay in comparison to linear decay close to the source and a deterioration of precision further away from the source, we would not expect any of the mentioned combinations to be of biological relevance — a prospect that may be too limited to warrant a discussion in the paper.
My main concern with the submission is that the article is written in a very short format that is very difficult for someone outside the field (myself included) to read. All of the information was there and I was able to figure things out, but only with significant effort. I would recommend expanding the text somewhat and providing a bit more detail as to where previous mathematical approaches ended and where this one starts.
We have now expanded the manuscript to include more information about the numerical simulations and the setup. A new panel in Figure 1, for example, illustrates the French flag model, making it easier for nonexperts to understand how morphogen gradients pattern a tissue. The newly added Figure 2 explains the entire modeling setup and how noise is introduced into the simulations. The previous Figure 2 has been turned into an appendix Figure 1, as it is concerned with a technical aspect that may be hard to understand by a general audience. We have also added graphical representations of the different boundary conditions to Figures 4 and 5 to make them more accessible. Moreover, we have included a freestanding methods section, further explaining how we generated the noisy gradients and discussing the choice of parameter distribution. Additionally, we have moved the majority of the mathematical derivations to a newly created appendix, which we believe makes the manuscript more accessible to a broader readership and allows the interested reader to look up the details. Finally, we explain in greater detail how the paper by Eldar et al. (2003) addresses noise and how our simulations differ from their view.
Overall, the manuscript was expanded to make it more accessible to a broader audience. Thank you for helping us improve the manuscript.
Reviewer #2 (Recommendations for the authors):
1. This is a bit of author preference, but I find that I prefer a freestanding methods section as opposed to having them interleaved into the results, as they are now. I think the freestanding section allows for a bit more elaboration on details without detracting from the narrative.
We have now included a methods section to provide a more detailed explanation of the technical aspects. In the methods section, we explain how the variable gradients were generated and how the choice of parameter distribution influences the results. We separated most of the mathematical derivations from the main text and moved it to a new appendix to make the main text more accessible to a wider audience.
Thank you for pointing this out. We agree that these steps helped make the manuscript more accessible.
2. I generally found this paper a nice exploration of the question of how nonlinear decay affects patterning. The paper is written in such a way that it may be challenging to read for those less familiar with the equations and conventions of mathematical modeling. If the authors wish to make the paper more accessible, I would suggest (in addition to the separation of methods) spending a bit more time defining the variables and parameters used throughout the paper, the justification for the different types of x/y axes used in the different graphs, and perhaps adding a cartoon or two illustrating a hypothetical tissue being patterned by a morphogen with some of the parameters/variables depicted. In addition, giving specific biological examples that correspond to some of the model considerations (e.g. boundary conditions) may be interesting for some readers.
We improved the manuscript in several ways. Firstly, we added a panel (Figure 1A) illustrating tissue patterning with the French flag model. Secondly, we added Figure 2, which describes the simulation setup and displays all relevant variables graphically. This makes it easier to understand how the simulations generate noisy gradients. The previous, more technical Figure 2 has been made an appendix Figure 1. Additionally, we moved the mathematical derivations to the appendix. We also added a discussion on how the choice of kinetic parameter distribution affects the results and how the variable gradients were generated to a dedicated methods section. To make the meaning of the boundary conditions clearer, we added small panels above Figures 4 and 5 that graphically illustrate them. Overall, we believe that these additions make the manuscript more informative and accessible to a broad readership also outside of the field of mathematical modeling. We thank the referee for pointing this out.
The main reason for choosing the boundary conditions is their mathematical convenience, which might also explain their widespread use elsewhere. Unfortunately, the boundary conditions in real tissues can be rather complicated as there is diffusion out of the domain, diffusion in 3D, and the boundary conditions might well differ between different systems. Here, we focused on the representation of the source in three different ways (explicit production domain, Dirichlet or vonNeumann boundary conditions as we now better illustrate (Figures 4 and 5)) while using zeroflux boundary conditions otherwise.
3. It might be interesting to add two points to the Discussion section. First, the simulations here explore the effect of parameter variation, which might be a good representation of, e.g. how genetic variability between individuals affects parameters and therefore patterning. However, even with consistent mean parameter values, molecular stochasticity can also impact these processes – is this worth exploring in future work? Second, given that most (all?) morphogens seem to operate in more complex regulatory networks, with features like opposing gradients or mutual inhibition of morphogen targets, might the interaction between nonlinear decay and these other features be worth exploring in the future?
We thank the referee for pointing out these aspects, which are indeed relevant and which require further research.
We assume that by “molecular stochasticity” the reviewer refers to the effect of a low number of individual morphogen molecules on patterning precision. In this work we assume that molecule numbers remain sufficiently large to allow the morphogen concentration to be represented by a continuous field variable C(x). Together with this simplification goes the assumption that cells are always able to read out a gradient concentration, even if it is very low, and that the concentration never drops to zero. Within this framework, molecular stochasticity, if interpreted in this way, is an aspect that lies outside the scope of our model.
At the moment there are only a few reliable measurements of morphogen gradients, and all that we are aware of are limited to relative (normalized) units. The absolute morphogen count remains elusive. We hope that the experimental community will soon be able to reliably measure the morphogen levels also in absolute terms (numbers of molecules). With such measurements, it would be feasible (and indeed important) to determine the influence of molecular stochasticity on gradient precision in future work. We discuss this now in the revised discussion.
In case that the referee is instead referring to stochasticity in the molecular nature of the morphogen kinetics, then this is an aspect that we hope to have made much clearer in the revised manuscript with the added illustrations and additional explanations. In our simulations, we represent such molecular stochasticity explicitly, and they are in fact at the heart of our model. Each cell in the simulated patterning domain possesses its own kinetic parameters p, d, and D, drawn independently from random distributions. It is this very aspect of our model that gives rise to noisy gradients, and to different gradients between embryos or tissues. We have revised the first section in Results, and added further details to the Methods section, to make this central aspect clearer. With the new Figure 2, we hope to also make this more visually apparent.
Regarding the referee’s second point: In this work, we focus solely on the precision of morphogen gradients and do not address their downstream readout in complex regulatory networks. Recently, it was demonstrated that opposing gradients are not required for high patterning precision in the mouse neural tube [Vetter et al. 2022, Nat. Commun.]. As we numerically show that close to the source, nonlinear decay does not lead to a significant increase in precision, and that precision deteriorates further away from the source, we have no reason to assume that interactions of nonlinear decay with the proposed features would significantly change the results. However, if network interactions are explicitly modeled, it would be straightforward to study the impact of nonlinear decay on gradient precision using the existing framework.
We have added a new paragraph to the discussion summarizing the two important points raised by the reviewer.
Reviewer #3 (Recommendations for the authors):
I find that this paper is very hard to read. The results of this paper are interesting, but, as written, this paper will appeal to a very narrow audience of researchers who have worked or work on modeling the production/diffusion/degradation model of morphogen gradients. This stems from two reasons. First, there is a lack of conceptual explanation of what the relevant questions really are. The authors refer to a previous paper but do not explain what the arguments in that paper really were. As far as I understand, the previous paper claimed that the nonlinear degradation will buffer against fluctuations in the levels of the source. That is interesting in principle, but there are other sources of noise and all of it is hard to get. I would urge the authors to rewrite their paper putting themselves in the shoes of someone who is not deeply familiar with all the details of the literature. The second reason why I find this paper confusing is that it is in part difficult to understand what the authors did. I had to check a previous paper to find some explanation of how molecular noise was introduced in the simulations. However, even after quickly checking that paper I remain confused about what the real assumptions were and how much I can trust the choice of parameters, etc. Again I would strongly recommend that the authors write a paper that is easier to read and contains all the relevant information.
In the revised manuscript, we now describe in greater detail how the previous paper treated variability in the morphogen gradients as a foldchange in the morphogen influx, and how this approach is limited, as the gradients only differ by the amplitudes/fluxes, but are not variable or noisy within the tissue. Our approach allows for a broader and more realistic range of sources of noise, which indeed exists, as the referee correctly notes. Our model introduces celltocell variability in all parameters (the production rate, degradation rate, diffusion coefficient, and cell areas) and also stochasticity in the morphogen influx or gradient amplitude, where applicable. We assume that these parameters follow lognormal distributions, and now provide a detailed explanation for this choice in the methods section.
We have also made several improvements to the manuscript to make it more accessible and selfcontained. Firstly, we added a panel explaining conceptually how a tissue is patterned according to the French flag model (Figure 1A). In the newly added Figure 2 we illustrate the computational model with all its variables, and how noise is introduced in the simulations. Secondly, we have now added a methods section that details how the noisy morphogen gradients were generated and how the parameter distributions were chosen. Thirdly, we split the mathematical derivations off from the Results section and moved them to the newly created appendix. Fourthly, to further facilitate interpretation of the results, we now normalize the positional error by the cell diameter, which we hope is an intuitive relative scale for a broad audience. Lastly, for an easier understanding of the different boundary conditions, we added small graphical depictions of them to Figures 4 and 5. All of these changes make the paper more selfcontained and accessible, hopefully allowing a much broader readership to follow the logic without consulting the references.
We thank the reviewer for helping us improve the manuscript.
https://doi.org/10.7554/eLife.84757.sa2Article and author information
Author details
Funding
Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (CRSII5_170930)
 Roman Vetter
 Dagmar Iber
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work was partially funded by SNF Sinergia grant CRSII5_170930.
Senior Editor
 Aleksandra M Walczak, CNRS, France
Reviewing Editor
 Gordon J Berman, Emory University, United States
Reviewer
 Gordon J Berman, Emory University, United States
Version history
 Preprint posted: November 4, 2022 (view preprint)
 Received: November 7, 2022
 Accepted: April 10, 2023
 Version of Record published: April 27, 2023 (version 1)
Copyright
© 2023, Adelmann et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 486
 Page views

 72
 Downloads

 4
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
A deep analysis of multiple genomic datasets reveals which genetic pathways associated with atherosclerosis and coronary artery disease are shared between mice and humans.

 Computational and Systems Biology
Mouse models have been used extensively to study human coronary artery disease (CAD) or atherosclerosis and to test therapeutic targets. However, whether mouse and human share similar genetic factors and pathogenic mechanisms of atherosclerosis has not been thoroughly investigated in a datadriven manner. We conducted a crossspecies comparison study to better understand atherosclerosis pathogenesis between species by leveraging multiomics data. Specifically, we compared genetically driven and thus CADcausal gene networks and pathways, by using human GWAS of CAD from the CARDIoGRAMplusC4D consortium and mouse GWAS of atherosclerosis from the Hybrid Mouse Diversity Panel (HMDP) followed by integration with functional multiomics human (STARNET and GTEx) and mouse (HMDP) databases. We found that mouse and human shared >75% of CAD causal pathways. Based on network topology, we then predicted key regulatory genes for both the shared pathways and speciesspecific pathways, which were further validated through the use of single cell data and the latest CAD GWAS. In sum, our results should serve as a muchneeded guidance for which human CADcausal pathways can or cannot be further evaluated for novel CAD therapies using mouse models.