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
Article 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.
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

 566
 Page views

 90
 Downloads

 6
 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
 Neuroscience
Closedloop neuronal stimulation has a strong therapeutic potential for neurological disorders such as Parkinson’s disease. However, at the moment, standard stimulation protocols rely on continuous openloop stimulation and the design of adaptive controllers is an active field of research. Delayed feedback control (DFC), a popular method used to control chaotic systems, has been proposed as a closedloop technique for desynchronisation of neuronal populations but, so far, was only tested in computational studies. We implement DFC for the first time in neuronal populations and access its efficacy in disrupting unwanted neuronal oscillations. To analyse in detail the performance of this activity control algorithm, we used specialised in vitro platforms with high spatiotemporal monitoring/stimulating capabilities. We show that the conventional DFC in fact worsens the neuronal population oscillatory behaviour, which was never reported before. Conversely, we present an improved control algorithm, adaptive DFC (aDFC), which monitors the ongoing oscillation periodicity and selftunes accordingly. aDFC effectively disrupts collective neuronal oscillations restoring a more physiological state. Overall, these results support aDFC as a better candidate for therapeutic closedloop brain stimulation.

 Cancer Biology
 Computational and Systems Biology
Currently, the identification of patientspecific therapies in cancer is mainly informed by personalized genomic analysis. In the setting of acute myeloid leukemia (AML), patientdrug treatment matching fails in a subset of patients harboring atypical internal tandem duplications (ITDs) in the tyrosine kinase domain of the FLT3 gene. To address this unmet medical need, here we develop a systemsbased strategy that integrates multiparametric analysis of crucial signaling pathways, and patientspecific genomic and transcriptomic data with a prior knowledge signaling network using a Booleanbased formalism. By this approach, we derive personalized predictive models describing the signaling landscape of AML FLT3ITD positive cell lines and patients. These models enable us to derive mechanistic insight into drug resistance mechanisms and suggest novel opportunities for combinatorial treatments. Interestingly, our analysis reveals that the JNK kinase pathway plays a crucial role in the tyrosine kinase inhibitor response of FLT3ITD cells through cell cycle regulation. Finally, our work shows that patientspecific logic models have the potential to inform precision medicine approaches.