Diffusion vs. direct transport in the precision of morphogen readout

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record updated: November 6, 2020 (This version)
Version of Record published: November 4, 2020 (Go to version)
Accepted Manuscript published: October 14, 2020 (Go to version)
Accepted: October 13, 2020
Received: May 15, 2020

1. Of interest
Substrate evaporation drives collective construction in termites

Giulio Facchini, Alann Rathery ... Andrea Perna

Research Article Apr 10, 2024
Further reading

Abstract
Introduction
Results
Discussion
Materials and methods
Appendix 1
Appendix 2
Appendix 3
Data availability
References
Article and author information
Metrics

Abstract

Morphogen profiles allow cells to determine their position within a developing organism, but not all morphogen profiles form by the same mechanism. Here, we derive fundamental limits to the precision of morphogen concentration sensing for two canonical mechanisms: the diffusion of morphogen through extracellular space and the direct transport of morphogen from source cell to target cell, for example, via cytonemes. We find that direct transport establishes a morphogen profile without adding noise in the process. Despite this advantage, we find that for sufficiently large values of profile length, the diffusion mechanism is many times more precise due to a higher refresh rate of morphogen molecules. We predict a profile lengthscale below which direct transport is more precise, and above which diffusion is more precise. This prediction is supported by data from a wide variety of morphogens in developing Drosophila and zebrafish.

Introduction

Within developing organisms, morphogen profiles provide cells with information about their position relative to other cells. Cells use this information to determine their position with extremely high precision (Dubuis et al., 2013; Erdmann et al., 2009; Gregor et al., 2007a; Houchmandzadeh et al., 2002; de Lachapelle and Bergmann, 2010). However, not all morphogen profiles are formed via the same mechanism and for some profiles the mechanism is still not well understood. One well-known mechanism is the synthesis-diffusion-clearance (SDC) model in which morphogen molecules are produced by localized source cells and diffuse through extracellular space before degrading or being internalized by target cells (Akiyama and Gibson, 2015; Gierer and Meinhardt, 1972; Lander et al., 2002; Müller et al., 2013; Rogers and Schier, 2011; Wilcockson et al., 2017). Alternatively, a direct transport (DT) model has been proposed where morphogen molecules travel through protrusions called cytonemes directly from the source cells to the target cells (Akiyama and Gibson, 2015; Bressloff and Kim, 2018; Kornberg and Roy, 2014; Müller et al., 2013; Wilcockson et al., 2017). The presence of these two alternative theories raises the question of whether there exists a difference in the performance capabilities between cells utilizing one or the other.

Experiments have shown that morphogen profiles display many characteristics consistent with the SDC model. The concentration of morphogen as a function of distance from the source cells has been observed to follow an exponential distribution for a variety of different morphogens (Driever and Nüsslein-Volhard, 1988; Houchmandzadeh et al., 2002). The accumulation times for several morphogens in Drosophila have been measured and found to match the predictions made by the SDC model (Berezhkovskii et al., 2011). In zebrafish, the molecular dynamics of the morphogen Fgf8 have been measured and found to be consistent with Brownian diffusion through extracellular space (Yu et al., 2009). Despite these consistencies, recent experiments have lent support to the theory that morphogen molecules are transported through cytonemes rather than extracellular space. The establishment of the Hedgehog morphogen gradient in Drosophila is highly correlated in both space and time with the formation of cytonemes (Bischoff et al., 2013), while Wnt morphogens have been found to be highly localized around cell protrusions such as cytonemes (Huang and Kornberg, 2015; Stanganello and Scholpp, 2016). Theoretical studies of both the SDC and DT models have examined these measurable effects (Berezhkovskii et al., 2011; Bressloff and Kim, 2018; Shvartsman and Baker, 2012; Teimouri and Kolomeisky, 2016a), but direct comparisons between the two models have thus far been poorly explored. In particular, it remains unknown whether one model allows for a cell to sense its local morphogen concentration more precisely than the other given biological parameters such as the number of cells or the characteristic lengthscale of the profile.

Here, we derive fundamental limits to the precision of morphogen concentration sensing for both the SDC and DT models. We investigate the hypothesis that sensory precision plays a major role in the selection of a gradient formation mechanism during evolution, and we test this hypothesis by quantitatively comparing our theory to morphogen data. Intuitively one might expect the DT model to have less noise due to the fact that molecules are directly deposited at their target. Indeed, we find below that the noise arises only from molecular production and degradation, with no additional noise from molecular transport. However, we also find below that for sufficiently large morphogen profile lengthscales, the SDC model produces less noise than the DT model due to it being able achieve a higher effective unique molecule count. By elucidating the competing effects of profile amplitude, steepness, and noise, we ultimately conclude that there should exist a profile lengthscale below which the DT model is more precise and above which the SDC mechanism is more precise. We find that this prediction is quantitatively supported by data from a wide variety of morphogens, suggesting that readout precision plays an important role in determining the mechanisms of morphogen profile establishment.

Results

Several past studies have focused on the formation dynamics of morphogen profiles (Berezhkovskii et al., 2011; Bressloff and Kim, 2018; Shvartsman and Baker, 2012; Teimouri and Kolomeisky, 2016a). Here, we model profiles in the steady state regime, as most of the experimental measurements to which we will later compare our results were taken during stages when the steady state approximation is valid (Grimm et al., 2010; Gregor et al., 2007b; Kicheva et al., 2007; Yu et al., 2009; Kanodia et al., 2009). Precision depends not only on stochastic fluctuations in the morphogen concentration, but also on the shape of the mean morphogen profile, as the shape determines concentration differences between adjacent cells that may adopt different fates. Therefore, as in past studies (Gregor et al., 2007a; Tostevin et al., 2007), we define the precision as $P = Δ m_{j} / σ_{j}$ , where $σ_{j}$ is the standard deviation of the number of morphogen molecules arriving at cell j, and $Δ m_{j} = m_{j} - m_{j + 1}$ is the difference between the molecule number in that cell and the adjacent cell. As is typical in studies of both the DT (Teimouri and Kolomeisky, 2016a; Bressloff and Kim, 2018) and SDC (Berezhkovskii et al., 2011; Shvartsman and Baker, 2012; Teimouri and Kolomeisky, 2016b) mechanisms, we focus on a one-dimensional line of target cells. However, we derive analogous results for 2D and 3D systems, and we generally find that the dimensionality does not qualitatively change our results, as we discuss later. In 1D, cells extend in both directions from the source cell, with N cells on each side (Figure 1). We note that in the case of the Bicoid morphogen in the Drosophila embryo, target cells extend only on one side of the source. This will introduce a factor of 2 in the means of both the DT and SDC models and a factor slightly greater than 2 in the variance of the SDC model. This will not affect the agreement of the Bicoid data with our theory in Figure 4B.

Figure 1

Download asset Open asset

Source cell (green) produces morphogen which is delivered to N target cells (blue) on either side via (A) direct transport (DT) or (B) synthesis-diffusion-clearance (SDC).

A common method for determining the statistical properties of stochastic systems is to express the dynamics of their probability distributions in the form of a master equation. The first moment of the master equation then dictates the dynamics of the mean and becomes the rate equation. Higher order moments can be similarly used to calculate other statistical properties such as the variance. As we will show, in the case of the DT mechanism only the rate equation will be needed as the relevant statistical properties are identical to that of a simple birth-death system which is fully characterized by its mean. Conversely, the SDC mechanism has more complicated noise properties, which we will calculate via the first and second moments. Specifically, we will use a Langevin description which is expressed as a rate equation with noise terms Gardiner, 2004.

Direct transport model

We first consider the DT case, where morphogen molecules are transported via cytonemes that connect a single source cell to multiple target cells (Figure 1A). Cytonemes are tubular protrusions that are hundreds of nanometers thick and between several and hundreds of microns long (Kornberg and Roy, 2014; Kornberg, 2014). They are supported by actin filaments, and it is thought that morphogen molecules are actively transported along the filaments via molecular motors (Kornberg and Roy, 2014; Kornberg, 2014; Sanders et al., 2013; Huang and Kornberg, 2015). It was recently shown that a DT model that includes forward and backward transport of molecules within cytonemes reproduces experimentally measured accumulation times (Teimouri and Kolomeisky, 2016a; Bressloff and Kim, 2018), although the noise properties of this model were not considered. Here, we review the steady state properties of this model and derive its noise properties.

Consider a single source cell that produces morphogen at rate $β$ . Morphogen molecules enter each cytoneme at rate $γ$ . The cytoneme that leads to the jth target cell has length $2 j a$ , where a is the cell radius. Once inside a cytoneme, morphogen molecules move forward toward the target cell with velocity $v^{+}$ or backwards toward the source cell with velocity $v^{-}$ , and can switch between these states with rates $ζ^{+}$ (forward-to-backward) or $ζ^{-}$ (backward-to-forward). Once a molecule reaches the forward (backward) end of the cytoneme, it is immediately absorbed into the target (source) cell. Molecules within a target cell spontaneously degrade with rate $ν$ . An alternative model could involve neglecting this degradation step by counting the arrival of morphogen molecules rather than the concentration. This method can at best reduce the variance by a factor of 2 as the noise from degradation is eliminated but the noise from arrival remains. This would cause negligible change to our results presented in Figure 3 and Figure 4 given the order of magnitude difference in precision seen between our two models.

The dynamics of the mean number of morphogen molecules in the source cell $m_{0} (t)$ and jth target cell $m_{j} (t)$ , and the mean density of forward-moving molecules $u_{j}^{+} (x, t)$ and backward-moving molecules $u_{j}^{-} (x, t)$ in the jth cytoneme are (Bressloff and Kim, 2018)

\begin{aligned} \frac{\partial m_{0}}{\partial t} = & β - \sum_{j} [γ m_{0} - v^{-} u_{j}^{-} (0, t)], \\ \frac{\partial u_{j}^{+}}{\partial t} = & - v^{+} \frac{\partial u_{j}^{+}}{\partial x} + ζ^{-} u_{j}^{-} - ζ^{+} u_{j}^{+}, \\ \frac{\partial u_{j}^{-}}{\partial t} = & v^{-} \frac{\partial u_{j}^{-}}{\partial x} - ζ^{-} u_{j}^{-} + ζ^{+} u_{j}^{+}, \\ \frac{\partial m_{j}}{\partial t} = & v^{+} u_{j}^{+} (L_{j}, t) - ν m_{j}, \end{aligned}

where the summation in the first line runs from $j = - N$ to $j = N$ excluding 0 (the source cell). Additionally, the boundary conditions $v^{+} u_{j}^{+} (0, t) = γ m_{0} (t)$ and $v^{-} u_{j}^{-} (L_{j}, t) = 0$ are imposed to reflect the rate at which morphogen molecules enter the cytoneme from either end. This creates a steady-state solution of the form

m_{j}^{D T} = \frac{β Γ_{j}}{2 ν \sum_{k = 1}^{N} Γ_{k}}, w h e r e Γ_{j} = \frac{e^{- 2 | j | κ a} (1 - e^{- ϕ})}{1 - e^{- ϕ - 2 | j | κ a}} .

This solution is identical to that found in Bressloff and Kim, 2018 with the exception of the factor of $1 / 2$ that accounts for target cells existing on both sides of the source cell in our model. Here, $γ Γ_{j}$ is the effective transport rate of morphogen molecules to the jth target cell, and $ϕ = \log (d^{-} / d^{+})$ and $κ = {(d^{+})}^{- 1} - {(d^{-})}^{- 1}$ are defined in terms of the average distance a molecule would move forward $d^{+} = v^{+} / ζ^{+}$ or backward $d^{-} = v^{-} / ζ^{-}$ within a cytoneme before switching direction. The parameter $ϕ$ sets the shape of $Γ_{j}$ , and thus of m_j: when $ϕ ≪ - 1$ the profile is constant, $Γ_{j} = 1$ ; when $ϕ ≫ 1$ it is exponential, $Γ_{j} = e^{- 2 | j | κ a}$ ; and when $| ϕ | ≪ 1$ it is a power law for large j, $Γ_{j} = {(1 + 2 | j | a / d^{+})}^{- 1}$ . The parameter $κ$ sets the lengthscale of the profile, defined as

λ_{DT} = \sum_{j = 1}^{N} \frac{Γ_{j} - Γ_{N}}{Γ_{1} - Γ_{N}} \approx \frac{1}{| κ |} (e^{| ϕ |} - 1) (| ϕ | - \log (e^{| ϕ |} - 1)),

where we approximate the sum as an integral for $N ≫ 1$ . We use this expression to eliminate $κ$ , writing $Γ_{j}$ in Equation 2 entirely in terms of $ϕ$ and $\hat{λ} \equiv λ / a$ .

Despite the complexity of the transport process in Equation 1, we find that it adds no noise to m_j. In fact, here we prove that any system in which molecules can only degrade in the target cells and cannot leave the target cells has the steady-state statistical properties of a simple birth-death process. First, consider the special case of only one target cell. Because each morphogen molecule produced in the source cell acts independently of every other morphogen molecule, we define $p (τ)$ as the probability density that any given molecule will enter the target cell a time $τ$ after it is created in the source cell. Next, we define $Q (δ t)$ as the probability that a morphogen molecule will enter the target cell between t and $t + δ t$ . This event requires the molecule to have been produced between $t - τ$ and $t - (τ + d τ)$ , which occurs with probability $β d τ$ ; to arrive at the target cell a time $τ$ later and to enter the target cell within the window $δ t$ , which occurs with probability $p (τ) δ t$ ; and we must integrate over all possible times $τ$ . Therefore,

Q (δ t) = \int_{0}^{\infty} [β d τ] [p (τ) δ t] = β δ t \int_{0}^{\infty} 𝑑 τ p (τ) = β δ t,

where the last step follows from normalization.

This generalized version of the DT model is visualized in Figure 2. We see that regardless of the form of $p (τ)$ , the probability of a morphogen molecule entering the target cell in any given small time window $δ t$ is simply $β δ t$ . Since the mechanism by which morphogen molecules go from the source cell to the target cell can only affect $p (τ)$ , this result holds regardless of the specifics of the mechanism so long as the condition that the morphogen cannot leave the target cell other than by degradation is maintained. This result also holds when the system is expanded to have multiple target cells, as then $p (τ)$ is replaced with $p_{j} (τ)$ , the probability density that the molecule enters the jth target cell a time $τ$ after being produced. In this case, $\int_{0}^{\infty} 𝑑 τ p_{j} (τ)$ evaluates to $π_{j}$ , the total probability the morphogen molecule is ultimately transported to the jth target cell, and $β δ t$ is simply replaced with $β π_{j} δ t$ . Combined with the constant degradation rate $ν$ of morphogen molecules within the target cell, this is precisely a birth-death process with birth rate $β π_{j}$ and death rate $ν$ . For our system $π_{j} = Γ_{j} / 2 \sum_{k = 1}^{N} Γ_{k}$ in Equation 2. The conditions that morphogen molecules act independently and can only degrade within the target cell are critical as the former ensures the noise sources are linear and the later will be violated by the SDC mechanism.

Figure 2

Download asset Open asset

Diagram outlining a generalized version of the DT model.

Between times $t - (τ + d τ)$ and $t - τ$ a molecule is produced with probability $β d τ$ . This molecule then undergoes a transport process over a time $τ$ with probability density $p (τ)$ . Finally, the molecule is deposited into the target cell at the end of the transport process. Integrating over all possible values of $τ$ then yields $Q (δ t)$ , the probability of the molecule entering the target cell between times t and $t + δ t$ .

We now assume that each cell integrates its morphogen molecule count over a time T (Berg and Purcell, 1977; Gregor et al., 2007a). The variance in the time average $T^{- 1} \int_{0}^{T} 𝑑 t m_{j} (t)$ is simply that of a birth-death process, given by $σ_{j}^{2} = 2 m_{j} / (T / τ)$ (Fancher and Mugler, 2017), so long as $T ≫ τ$ , where $τ = ν^{- 1}$ is the correlation time. We see that, as expected for a time-averaged Poisson process, the variance increases with the mean m_j and decreases with the number $T / τ$ of independent measurements made in the time T. The precision is therefore

P_{DT}^{2} = \frac{m_{j}^{DT} T}{2 τ_{DT}} {(\frac{Δ m_{j}^{DT}}{m_{j}^{DT}})}^{2}, with τ_{DT} = \frac{1}{ν} .

We see that the precision increases with the profile amplitude m_j, the number of independent measurements $T / τ$ , and the profile steepness $Δ m_{j} / m_{j}$ . The transport process influences the precision only via m_j, not $τ$ . For a given N, j, and $\hat{λ}$ , we find that the precision is maximized at a particular $ϕ^{*} > 0$ (Figure 3A). The reason is that an exponential profile ( $ϕ ≫ 1$ ) has constant steepness but small amplitude, whereas a power-law profile ( $ϕ ≪ 1$ ) has low steepness but large amplitude due to its long tail; the optimum is in between.

Figure 3

Download asset Open asset

Comparing theoretical DT precision to SDC precision for a single cell.

(A) DT precision shows a maximum as a function of shape parameter $ϕ$ for any value of the profile lengthscale. (B) Ratio $ρ_{j}$ of DT to SDC precision shows a crossover ( $ρ_{j} = 1$ ) as a function of profile lengthscale $λ / a$ for 1D, 2D, and 3D geometries. Here $j = 50$ is the central cell of $N = 100$ target cells. For each value of $\hat{λ}$ the value of $ϕ$ which maximizes precision in the DT model ( $ϕ^{*}$ ) as seen in A is used.

SDC model

We next consider the SDC case (Figure 1B). Again a single source cell at the origin $x = 0$ produces morphogen at rate $β$ . However, now morphogen molecules diffuse freely along x with coefficient D and degrade spontaneously at any point in space with rate $ν$ . The dynamics of the morphogen concentration $c (x, t)$ are

\frac{\partial c}{\partial t} = D \nabla^{2} c + η_{D} - ν c - η_{ν} + (β + η_{β}) δ (x),

where the noise terms associated with diffusion, degradation, and production obey

\begin{aligned} ⟨ η_{D} (x^{'}, t^{'}) η_{D} (x, t) ⟩ & = 2 D δ (t - t^{'}) {\vec{\nabla}}_{x} \cdot {\vec{\nabla}}_{x^{'}} c (x) δ (x - x^{'}) \\ ⟨ η_{ν} (x^{'}, t^{'}) η_{ν} (x, t) ⟩ & = ν c (x) δ (t - t^{'}) δ (x - x^{'}), \\ ⟨ η_{β} (t^{'}) η_{β} (t) ⟩ & = β δ (t - t^{'}), \end{aligned}

respectively (Gardiner, 2004; Gillespie, 2000; Fancher and Mugler, 2017; Varennes et al., 2017). Here, the time independent $c (x) = β e^{- | x | / λ} / (2 ν λ)$ is the steady state mean concentration, with characteristic lengthscale $λ_{SDC} = \sqrt{D / ν}$ . We imagine a target cell located at x that is permeable to the morphogen and counts the number $m (x, t) = \int_{V} 𝑑 y c (x + y, t)$ of morphogen molecules within its volume V. We use this simpler prescription over explicitly accounting for more realistic mechanisms such as surface receptor binding because it has been shown that the two approaches ultimately yield similar concentration sensing results up to a factor of order unity (Berg and Purcell, 1977). For a cell in steady state at position $x = 2 j a$ , the integral evaluates to

m_{j}^{SDC} = (β / ν) \sinh (1 / \hat{λ}) e^{- 2 | j | / \hat{λ}} .

Importantly, in the model presented here the morphogen molecules can diffuse both into and out of the target cells, thus violating the condition that they can only degrade once in a target cell and disallowing our previous argument depicted in Figure 2. However, because Equation 6 is linear with Gaussian white noise, calculating the time-averaged variance $σ_{j}^{2}$ is straightforward: we Fourier transform Equation 6 in space and time, calculate the power spectrum of $m (x, t)$ , and take its low-frequency limit (Appendix 1). So long as $T ≫ ν^{- 1}$ , we obtain the same functional form as Equation 5,

P_{SDC}^{2} = \frac{m_{j}^{SDC} T}{2 τ_{SDC}} {(\frac{Δ m_{j}^{SDC}}{m_{j}^{SDC}})}^{2},

because diffusion is a Poisson process. This result is distinguished from that of the DT system by the correlation time taking the form

τ_{SDC} = \frac{1}{ν} [1 - \frac{(2 / \hat{λ}) + \sinh (2 / \hat{λ})}{4 \sinh (1 / \hat{λ}) e^{1 / \hat{λ}}}] .

The factor in brackets is always less than one and decreases with $\hat{λ}$ . It reflects the fact that, unlike in the DT model, molecules can leave a target cell not only by degradation, but also by diffusion. Therefore, the rate $τ^{- 1}$ at which molecules are refreshed is larger than that from degradation alone. This effect increases the precision because more independent measurements ( $T / τ$ ) can be made.

To understand this effect more intuitively, consider a simplified SDC model in which diffusion is modeled as discrete hopping between adjacent target cells at rate h. The autocorrelation function is $C_{j} (t) = m_{j} I_{0} (2 h t) e^{- (2 h + ν) t}$ (Appendix 2), where I₀ is the zeroth modified Bessel function of the first kind. The correlation time is $τ = \int_{0}^{\infty} 𝑑 t C_{j} (t) / C_{j} (0) = {[ν (4 h + ν)]}^{- 1 / 2}$ , and we see explicitly that it decreases with both degradation ( $ν$ ) and diffusion (h). In fact, in the limit of fast diffusion ( $h ≫ ν$ ), the expression becomes $τ = {(4 ν h)}^{- 1 / 2}$ . Correspondingly, in the fast-diffusion limit of Equation 9 ( $\hat{λ} ≫ 1$ ), the term in brackets reduces to ${\hat{λ}}^{- 1}$ , and it becomes $τ = {(ν \hat{λ})}^{- 1 / 2} = {[4 ν D / {(2 a)}^{2}]}^{- 1 / 2}$ . These expressions are identical, with $D / {(2 a)}^{2}$ playing the role of the hopping rate h, as expected.

Comparing the models

We now ask which model has the higher precision. We first note that while the precision in both models has the same form when expressed in terms of the means and correlation times (Equations 5 and 8), substituting in the corresponding expressions reveals how the precision depends on the various parameters of each model. Specifically, the precision in the DT model is seen to depend on N, j, the product $β T$ , $\hat{λ}$ , and $ϕ$ . Importantly, it is independent of $ν$ so long as the condition $T ≫ τ_{DT} = ν^{- 1}$ is met. This is due to the mean molecule count scaling as $ν^{- 1}$ (Equation 2) and the number of independent measurements ( $T / τ_{DT} = ν T$ ) scaling as $ν$ , thus causing their product and in turn the precision to be independent of $ν$ . In the SDC model, the precision depends on j, the product $β T$ , and $\hat{λ}$ . For similar reasons as in the DT model, $ν$ does not explicitly appear once Equations 7 and 9 are inserted into Equation 8, but it is present implicitly in the definition of $\hat{λ} = \sqrt{D / ν} / a$ .

To properly compare the two models, we equate several variables. Specifically, we wish to compare the precision of each model within a specific system, which requires N and j to be the same in both models. Additionally, $β$ and $ν$ are restricted by the energy and material costs of producing and degrading morphogen molecules while T is restricted by the need of the system to properly develop in a finite amount of time. These restrictions are assumed to be independent of the method of morphogen profile establishment, which implies that both the DT and SDC systems will adopt the same maximal values of $β$ and T. While the value of $ν$ is restricted in a similar manner, explicitly equating it between the two models does not reduce the number of free parameters, as the precision is independent of $ν$ in the DT model, and $ν$ is insufficient to fully define the value of $\hat{λ}$ in the SDC model. Therefore, we equate the characteristic profile length scale $\hat{λ}$ , as this parameter is consistently defined across both models and is also measured in experiments, allowing us to compare our theory with data as discussed below. Finally, we optimize over $ϕ$ as it is unique to the DT model.

We now observe how the precision of the DT and SDC models compare in a representative system. Figure 3B shows $ρ_{j} = P_{DT}^{2} / P_{SDC}^{2}$ as a function of profile length $\hat{λ}$ for a cell in the center ( $j = N / 2$ ) of a line of $N = 100$ target cells, where for each $\hat{λ}$ we use the $ϕ^{*}$ that maximizes $P_{DT}^{2}$ as seen in Figure 3A. Since $β$ and T have been equated between the two models, $ρ_{j}$ is independent of both. We see that for short profiles the DT model is more precise ( $ρ_{j} > 1$ ) whereas for long profiles the SDC model is more precise ( $ρ_{j} < 1$ ). This effect holds for a single source cell providing morphogen for a 1D line of target cells as well as for a 1D line of source cells with a 2D sheet of target cells and a 2D sheet of source cells with a 3D volume of target cells. For the DT model, the 2D and 3D cases are identical to the 1D case as we assume that cytonemes extend perpendicular to the source cells; for the SDC model see Appendix 1.

The reason that the SDC model is more precise for long profiles is that long profiles correspond to fast diffusion, which increases the refresh rate $τ_{SDC}^{- 1}$ as discussed above. Conversely, the reason that the DT model is more precise for short profiles is that it has a larger amplitude. It also has a smaller steepness, but the larger amplitude wins out. Specifically, whereas the SDC amplitude falls off exponentially, $m_{j} \sim e^{- 2 | j | / \hat{λ}}$ , for sufficiently small $ϕ^{*}$ the DT amplitude falls off as a power law, $m_{j} \sim 1 / | j |$ . The steepness $Δ m_{j} / m_{j}$ of the SDC profile is constant, while the steepness of the DT profile also scales like $1 / j$ . Thus, the product of the ratio of amplitudes and the square of the ratio of steepnesses, on which $ρ_{j}$ depends, scales like $e^{2 | j | / \hat{λ}} / {| j |}^{3}$ . For small $\hat{λ}$ , the exponential dominates over the cubic for the majority of j values. Consequently, the DT model has the higher precision.

Comparison to data

We now test our predictions against data for various morphogens. In Drosophila, the morphogen Wingless (Wg) is localized near cell protrusions such as cytonemes (Huang and Kornberg, 2015; Stanganello and Scholpp, 2016), and the Hedgehog (Hh) gradient correlates highly in both space and time with the formation of cytonemes (Bischoff et al., 2013), suggesting that these two morphogen profiles are formed via a DT mechanism. Conversely, Bicoid has been understood as a model example of SDC for decades (Driever and Nüsslein-Volhard, 1988; Gregor et al., 2007a; Houchmandzadeh et al., 2002). Similarly, Dorsal is spread by diffusion; however, its absorption is localized to a specific region of target cells via a nonuniform degradation mechanism, making it more complex than the simple SDC model (Carrell et al., 2017). Finally, for Dpp there is evidence for a variety of different gradient formation mechanisms (Akiyama and Gibson, 2015; Müller et al., 2013; Wilcockson et al., 2017).

In zebrafish, the morphogen Fgf8 has been studied at the single molecule level and found to have molecular dynamics closely matching the Brownian movement expected in an SDC mechanism (Yu et al., 2009). Similarly, Cyclops, Squint, Lefty1, and Lefty2, all of which are involved in the Nodal/Lefty system, have been shown to spread diffusively and affect cells distant from their source (Müller et al., 2013; Rogers and Müller, 2019). This would support the SDC mechanism, although Cyclops and Squint have been argued to be tightly regulated via a Gierer-Meinhardt type system, thus diminishing their gradient sizes to values much lower than what they would be without this regulation (Gierer and Meinhardt, 1972; Rogers and Müller, 2019).

For all these morphogens, we estimate the profile lengthscales $λ$ from the experimental data (Kicheva et al., 2007; Wartlick et al., 2011; Gregor et al., 2007b; Liberman et al., 2009; Yu et al., 2009; Müller et al., 2012) (Appendix 3). Figure 4A shows these $λ$ values and indicates for each morphogen whether the evidence described above suggests a DT mechanism (red), an SDC mechanism (blue), or multiple mechanisms including DT and SDC (white). We see that in general, the three cases correspond to short, long, and intermediate profile lengths, respectively, which is qualitatively consistent with our predictions.

Figure 4

Download asset Open asset

Comparing theory and experiment.

(A) $λ$ values for morphogens estimated from experiments, colored by whether experiments support a DT (red), SDC (blue), or multiple mechanisms (white). (B) Data from A overlaid with color from theory using values of a and N estimated from experiments. Color indicates percentage of cells for which SDC is predicted to be more precise with dots signifying 0% and 100%, dashed lines signifying 25% and 75%, and solid lines signifying 50%. The $\hat{λ}$ axis is normalized by ${\hat{λ}}_{50}$ , the value of $\hat{λ}$ at which 50% of cells are more precise in SDC.

To make the comparison quantitative, we estimate the values of cell radius a and cell number N from the experimental data (Kicheva et al., 2007; Gregor et al., 2007a; Liberman et al., 2009; Yu et al., 2009; Kimmel et al., 1995) (Appendix 3) in order to calculate $ρ_{j}$ from our theory in each case. The background color in Fig. Figure 4B shows the percentage of cells within each system for which we predict that the SDC model is more precise as a function of $\hat{λ}$ . The data points in Fig. Figure 4B show the values of $\hat{λ}$ from the experiments, also normalized by ${\hat{λ}}_{50}$ , the value of $\hat{λ}$ at which 50% of cells are more precise in SDC, from the theory. For each morphogen species, we assume a 1D system for simplicity as we have checked that considering higher dimensions yields negligible differences to the results presented in Figure 4B. We see that our theory predicts the correct threshold: the morphogens for which the evidence suggests either a DT or an SDC mechanism (red or blue) fall into the regime in which we predict that mechanism to be more precise for most of the cells, and the morphogens with multiple mechanisms (white) fall in between. This result provides quantitative support for the idea that morphogen profiles form according to the mechanism that maximizes the sensory precision of the target cells.

Discussion

We have shown that in the steady-state regime, the DT and SDC models of morphogen profile formation yield different scalings of readout precision with the length of the profile and population size. As a result, there exist regimes in this parameter space in which either mechanism is more precise. While the DT model benefits from larger molecule numbers and no added noise from the transport process, the ability of molecules to diffuse into and away from a target cell in the SDC model allows the cell to measure a greater number of effectively unique molecules in the same time frame. By examining how these phenomena affect the cells’ sensory precision, we predicted that morphogen profiles with shorter lengths should utilize cytonemes or some other form of direct transport mechanism, whereas morphogens with longer profiles should rely on extracellular diffusion, a prediction that is in quantitative agreement with measurements on known morphogens. It will be interesting to observe whether this trend is further strengthened as more experimental evidence is obtained for different morphogens, as well as to expand the theory of multicellular concentration sensing to further biological contexts.

Despite the quantitative agreement between our theory and experiments, it is clear that the models presented here are minimal and thus cannot be directly applied to all systems. This is exemplified by morphogen such as Dorsal, which due to aforementioned diffusive spreading and nonuniform degradation mechanism clearly does not strictly follow either model. Additionally, the SDC model can be violated if the diffusion of morphogen through a biological environment is hindered by the typically crowded nature of such environments, leading to possibly subdiffusive behavior (Ellery et al., 2014; Fanelli and McKane, 2010). For the DT model, we explicitly ignored the dynamics of the cytonemes themselves due to the growth rate of the cytonemes being sufficiently fast so as to traverse the entire system size in significantly less time than is required for the cells to integrate their morphogen counts over (Bischoff et al., 2013; Chen et al., 2017). This assumption is problematic if cytonemes continue to behave dynamically after reaching the source cell. In particular, the process of cytonemes switching between phases of growing and retracting can introduce super-Poissonian noise sources to the morphogen count within the target cells. Super-Poissonian noise can also be introduced by relaxing the assumption that the morphogen molecules behave independently of each other as this was a critical component of our proof that the molecule count in the target cells is statistically identical to a birth-death process. Finally, it is conceivable that a hybrid system which utilizes a combination of diffusion and directed transport could be developed. However, we are unaware of any experimental study into such a possibility and as such have not considered it in this particular work. It will be interesting to explore the implications of each of these complications in future works.

Materials and methods

Methods are described in Appendices 1, 2, and 3 along with more detailed evaluations of 2D and 3D SDC systems and further discussion of data used in Figure 4. Additionally, code used to generate the plots seen in Figure 3 and Figure 4 can be found at Fancher, 2020.

Appendix 1

Time-averaged variance in the SDC model

Here we calculate the time-averaged variance of the morphogen molecule number using the low-frequency limit of the power spectrum. We first introduce the power spectrum, and then we calculate the variance for the 1D, 2D, and 3D geometries.

Power spectrum

We first discuss the correlation function and power spectrum to establish some definitions and notation. Specifically, we show that the variance in the long-time average of a variable is given by the low-frequency limit of its power spectrum. For a one-dimensional function $x (t)$ with mean 0, the correlation function $C (t)$ takes the form

C (t - t^{'}) = ⟨ x (t^{'}) x (t) ⟩ .

Since absolute time is irrelevant in the steady state of any physical system with no time dependent forcing, $t^{'}$ can be set to 0 without loss of generality. This leads to a definition for the power spectrum of $x (t)$ as

\begin{aligned} S (ω) & = \int \frac{d ω^{'}}{2 π} ⟨ {\tilde{x}}^{*} (ω^{'}) \tilde{x} (ω) ⟩ = \frac{1}{2 π} \int d ω^{'} d t d t^{'} ⟨ x (t^{'}) x (t) ⟩ e^{i ω t} e^{- i ω^{'} t^{'}} \\ = \int d t d t^{'} C (t - t^{'}) e^{i ω t} δ (t^{'}) = \int d t C (t) e^{i ω t} . \end{aligned}

Thus, under this definition the power spectrum is seen to be the Fourier transform of the correlation function. Additionally, when $x (t)$ is averaged over a time T, the time averaged correlation function of $x (t)$ takes the form

\begin{aligned} C_{T} (t - t^{'}) & = ⟨ (\frac{1}{T} \int_{t^{'}}^{t^{'} + T} d τ^{'} x (τ^{'})) (\frac{1}{T} \int_{t}^{t + T} d τ x (τ)) ⟩ \\ = \frac{1}{T^{2}} \int_{t}^{t + T} d τ \int_{t^{'}}^{t^{'} + T} d τ^{'} ⟨ x (τ^{'}) x (τ) ⟩ \\ = \frac{1}{T^{2}} \int_{t}^{t + T} d τ \int_{t^{'}}^{t^{'} + T} d τ^{'} C (τ - τ^{'}) . \end{aligned}

Let $y \equiv (τ - τ^{'}) - (t - t^{'})$ and $z \equiv (τ + τ^{'}) - (t + t^{'})$ . This transforms Equation 13 into

\begin{aligned} C_{T} (t - t^{'}) & = \frac{1}{T^{2}} \int_{- T}^{T} d y \int_{| y |}^{2 T - | y |} d z \frac{1}{2} C (y + t - t^{'}) \\ = \frac{1}{T^{2}} \int_{- T}^{T} d y (T - | y |) C (y + t - t^{'}) . \end{aligned}

By inverting the relationship found in Equation 12, $C (y + t - t^{'})$ can be replaced with an inverse Fourier transform of $S (ω)$ to produce

\begin{aligned} C_{T} (t - t^{'}) & = \frac{1}{T^{2}} \int_{- T}^{T} d y \int \frac{d ω}{2 π} (T - | y |) S (ω) e^{- i ω (y + t - t^{'})} \\ = \int \frac{d ω}{2 π} {(\frac{2}{ω T} \sin (\frac{ω T}{2}))}^{2} S (ω) e^{- i ω (t - t^{'})} . \end{aligned}

The factor of ${(ω T)}^{- 2}$ in the integrand of Equation 15 forces only small values of $ω$ to contribute when T is large. Thus, the approximation $S (ω) \approx S (0)$ can be made since only values of $ω$ near 0 are contributing. This causes $C_{T} (0)$ , which we will denote as $σ^{2}$ through this and the main text, to be exactly calculable to

σ^{2} = C_{T} (0) \approx S (0) \int \frac{d ω}{2 π} {(\frac{2}{ω T} \sin (\frac{ω T}{2}))}^{2} = \frac{S (0)}{T} .

Of important note is the fact that this approximation only works if $S (ω)$ varies slowly compared to ${(2 \sin (ω T / 2) / ω T)}^{2}$ near $ω = 0$ . Since $C (t)$ must be time symmetric, $S (ω)$ must also be symmetric and thus an even function of $ω$ . Thus, near $ω = 0$ the lowest order correction term for each function will be the second order term. Normalizing each term by the 0-frequency value of each function then lets us to impose the condition

| \frac{1}{S (0)} \frac{\partial^{2} S (ω)}{\partial ω^{2}} |_{ω = 0} | ≪ | \frac{\partial^{2}}{\partial ω^{2}} {(\frac{2}{ω T} \sin (\frac{ω T}{2}))}^{2} |_{ω = 0} | = \frac{T^{2}}{6} .

So long as this condition is satisfied, the approximation given in Equation 16 is valid.

We now cast Equation 16 into a more intuitive form by considering the correlation time $τ$ , which can be defined as

τ = \int_{0}^{\infty} 𝑑 t \frac{C (t)}{C (0)} .

Continuing the use the fact that $C (t)$ must be time symmetric and thus an even function of t, Equation 12 can be used to produce the result

S (0) = \int 𝑑 t C (t) = 2 \int_{0}^{\infty} 𝑑 t C (t) = 2 τ C (0) .

Inserting this result into Equation 16 produces

σ^{2} \approx \frac{2 τ}{T} C (0),

thus relating the long-time averaged variance, $σ^{2}$ , to the instantaneous variance, $C (0)$ , and the number of correlation times the system averages over, $T / τ$ .

Variance and precision

We now consider a model for the Synthesis-Diffusion-Clearance system. We still assume there is a single source cell which produces morphogen at rate $β$ , but now the morphogen is released into the extracellular environment where it freely diffuses at rate D. The morphogen can also spontaneously degrade at rate $ν$ . Even though in the main text we focus on a zero-dimensional source in a one-dimensional space, here we will look at diffusion in a multitude of different spaces with different dimensions as well as morphogen sources that span a multitude of different dimensions. In each case, the sources will secrete morphogen molecules into a density field c which must follow

\frac{\partial c}{\partial t} = D \nabla^{2} c + η_{D} - ν c - η_{ν} + (β + η_{β}) δ^{S P - S O} (\vec{x}),

where $S P$ is the number of spatial dimensions, $S O$ is the dimensionality of the source, and $\nabla^{2}$ is taken over all $S P$ dimensions. Each $η$ term is a Langevin noise term that represents Gaussian white noise for the diffusion, degradation, and production processes respectively. Of important note is that $δ^{S P - S O} (\vec{x})$ is a $δ$ function only in the last $S P - S O$ dimensions of the space. So, for example, if there was a one dimensional source in three dimensional space, then $δ^{3 - 1} (\vec{x})$ would be a $δ$ function in the $\hat{y}$ and $\hat{z}$ directions but not the $\hat{x}$ direction. This means that $β$ and $η_{β}$ will have units of $T^{- 1} L^{- S O}$ , where T is time and L is space.

We can now assume c has reached a steady state and separate it into $c = \bar{c} + δ c$ , which in turn allows Equation 21 to separate into

0 = D \nabla^{2} \bar{c} - ν \bar{c} + β δ^{S P - S O} (\vec{x})

\frac{\partial δ c}{\partial t} = D \nabla^{2} δ c + η_{D} - ν δ c - η_{ν} + η_{β} δ^{S P - S O} (\vec{x}) .

where

Fourier transforming Equation 22 in space and dividing it by $ν$ then yields

0 = - λ^{2} {| \vec{k} |}^{2} \tilde{\bar{c}} - \tilde{\bar{c}} + \frac{β}{ν} {(2 π)}^{S O} δ^{S O} (\vec{k}) ⟹ \tilde{\bar{c}} = \frac{β}{ν} \frac{{(2 π)}^{S O} δ^{S O} (\vec{k})}{1 + λ^{2} {| \vec{k} |}^{2}},

where

λ = \sqrt{\frac{D}{ν}} .

Of similarly important note is that $δ^{S O} (\vec{k})$ is a $δ$ function only in the first $S O$ dimensions of k-space. So in the one-dimensional source, three-dimensional space example $δ^{S O} (\vec{k})$ would be a $δ$ function in the $\hat{x}$ direction of k-space but not the $\hat{y}$ or $\hat{z}$ directions.

This allows $\bar{c}$ to be written as

\begin{aligned} \bar{c} (\vec{x}) & = \int \frac{d^{S P} k}{{(2 π)}^{S P}} e^{- i \vec{k} \cdot \vec{x}} \tilde{\bar{c}} (\vec{k}) = \frac{β}{ν} \int \frac{d^{S P} k}{{(2 π)}^{S P}} e^{- i \vec{k} \cdot \vec{x}} \frac{{(2 π)}^{S O} δ^{S O} (\vec{k})}{1 + λ^{2} {| \vec{k} |}^{2}} \\ = \frac{β}{ν} \int \frac{d^{S P - S O} k}{{(2 π)}^{S P - S O}} e^{- i \vec{k} \cdot \vec{x}} \frac{1}{1 + λ^{2} {| \vec{k} |}^{2}} = \frac{β λ^{- (S P - S O)}}{ν} P_{S P - S O} (\frac{| \vec{x} |}{λ}), \end{aligned}

where

P_{N} (x) = \int \frac{d^{N} u}{{(2 π)}^{N}} e^{- i \vec{u} \cdot \vec{x}} \frac{1}{1 + {| \vec{u} |}^{2}} .

It is important to note that $P_{N}$ does not integrate over all available dimensions, but only over the last N dimensions of the space. This in turn means that its argument can only depend on the last N dimensions of any input vector. Returning to the one dimensional source, three dimensional space example, $P_{3 - 1} (| \vec{x} | / λ)$ should only take the y and z components of $\vec{x}$ into account. The x component is made irrelevant by the translational symmetry of the system along the x-axis.

Moving on to the noise terms, Equation 23 can be Fourier transformed in space and time to yield

- i ω \tilde{δ c} = - D {| \vec{k} |}^{2} \tilde{δ c} + {\tilde{η}}_{D} - ν \tilde{δ c} - {\tilde{η}}_{ν} + {\tilde{η}}_{β} ⟹ \tilde{δ c} = \frac{{\tilde{η}}_{D} - {\tilde{η}}_{ν} + {\tilde{η}}_{β}}{ν (1 + λ^{2} {| \vec{k} |}^{2} - i \frac{ω}{ν})},

where $η_{β} (\vec{k}, ω)$ depends only on the first $S O$ dimensions of k-space. Assuming the $η$ terms are all independent of each other allows the cross spectrum of c to be

\begin{aligned} ⟨ {\tilde{δ c}}^{*} ({\vec{k}}^{'}, ω^{'}) \tilde{δ c} (\vec{k}, ω) ⟩ = \frac{1}{ν^{2} (1 + λ^{2} {| \vec{k} |}^{2} - i \frac{ω}{ν}) (1 + λ^{2} {| {\vec{k}}^{'} |}^{2} + i \frac{ω^{'}}{ν})} \\ \cdot (⟨ {\tilde{η}}_{D}^{*} ({\vec{k}}^{'}, ω^{'}) {\tilde{η}}_{D} (\vec{k}, ω) ⟩ + ⟨ {\tilde{η}}_{ν}^{*} ({\vec{k}}^{'}, ω^{'}) {\tilde{η}}_{ν} (\vec{k}, ω) ⟩ + ⟨ {\tilde{η}}_{β}^{*} ({\vec{k}}^{'}, ω^{'}) {\tilde{η}}_{β} (\vec{k}, ω) ⟩) . \end{aligned}

The cross spectrum of $η_{D}$ can be obtained from its correlation function. To derive such a correlation function, we first consider a separate Markovian system comprised of a 1-dimensional lattice of discrete compartments that a diffusing species Y can exist in. The dimensionality is chosen purely for simplicity, as the method outlined below can be easily generalized to higher dimensions to produce the same result. Let $y_{i} (t)$ be the number of Y molecules in the ith compartment at time t and d be the rate at which these molecules move to the $i - 1$ or $i + 1$ compartment. Given a sufficiently small time step $δ t$ , the probability of a molecule moving from the ith compartment to the $i \pm 1$ compartment is

P ({y_{i} (t + δ t), y_{i \pm 1} (t + δ t)} = {y_{i} (t) - 1, y_{i \pm 1} (t) + 1}) = y_{i} (t) d δ t .

Higher order interactions in which multiple molecules are transfered within the time step $δ t$ will have probabilites of order ${(δ t)}^{2}$ or higher and can thus be ignored. This allows the mean of $δ y_{i} (t) = y_{i} (t + δ t) - y_{i} (t)$ to take the form

⟨ δ y_{i} (t) ⟩ = (y_{i - 1} (t) + y_{i + 1} (t) - 2 y_{i} (t)) d δ t,

where the first two terms come from molecules moving into theith compartment from the $i - 1$ and $i + 1$ compartments respectively and the third term comes from the two different ways molecules can leave the ith compartment. As $δ t$ is small, each of these transfer processes can be treated as being Poissonianly distributed. This allows the variance of $δ y_{i} (t)$ to simply be the right-hand side of Equation 24 but with each term taken to be its absolute value so there are no subtractions. Additionally, this approximation allows the covariance between $δ y_{i}$ and $δ y_{i \pm 1}$ to be taken as the negative of the sum of the expected number of molecules moving from the ith compartment to the $i \pm 1$ compartment and vice versa. With these, the correlation function between $δ y_{i} (t)$ and $δ y_{j} (t)$ can be written as

⟨ δ y_{j} (t) δ y_{i} (t) ⟩ = (y_{i - 1} (t) + y_{i + 1} (t) + 2 y_{i} (t)) d δ t δ_{i, j} - (y_{i} (t) + y_{i - 1} (t)) d δ t δ_{i - 1, j} - (y_{i} (t) + y_{i + 1} (t)) d δ t δ_{i + 1, j} .

We now take the system to continuous space by letting $y_{i} (t) \to ℓ c (x, t)$ and $δ_{i, j} \to ℓ δ (x - x^{'})$ with any intances of ±1 in the indices also being converted to $\pm ℓ$ . Putting these substitutions into Equation 25 and dividing by ${(ℓ δ t)}^{2}$ yields

\begin{array}{ll} ⟨ \frac{δ c (x^{'}, t)}{δ t} \frac{δ c (x, t)}{δ t} ⟩ = \frac{d}{δ t} ((c (x - ℓ, t) + c (x + ℓ, t) + 2 c (x, t)) δ (x - x^{'}) \\ - (c (x, t) + c (x - ℓ, t)) δ (x - ℓ - x^{'}) - (c (x, t) + c (x + ℓ, t)) δ (x + ℓ - x^{'})) \\ = \frac{d}{δ t} ((c (x + ℓ, t) δ (x - x^{'}) - c (x + ℓ, t) δ (x + ℓ - x^{'})) \\ - (c (x, t) δ (x - ℓ - x^{'}) - c (x, t) δ (x - x^{'})) (c (x, t) δ (x - x^{'}) \\ - c (x, t) δ (x + ℓ - x^{'})) - (c (x - ℓ, t) δ (x - ℓ - x^{'}) - c (x - ℓ, t) δ (x - x^{'}))) . \end{array}

Equation 33 has been rearranged into this form so as to easily apply the operators $\partial_{x}^{\pm}$ defined as

\partial_{x}^{+} f (x) = \frac{f (x + ℓ) - f (x)}{ℓ},

\partial_{x}^{-} f (x) = \frac{f (x) - f (x - ℓ)}{ℓ} .

Using this notation, Equation 33 can be simplfied into

\begin{aligned} ⟨ \frac{δ c (x^{'}, t)}{δ t} \frac{δ c (x, t)}{δ t} ⟩ & = \frac{ℓ d}{δ t} (\partial_{x}^{+} (c (x, t) δ (x - ℓ - x^{'}) - c (x, t) δ (x - x^{'})) \\ + \partial_{x}^{-} (c (x, t) δ (x - x^{'}) - c (x, t) δ (x + ℓ - x^{'}))) \\ = \frac{ℓ^{2} d}{δ t} (\partial_{x}^{+} \partial_{x^{'}}^{+} + \partial_{x}^{-} \partial_{x^{'}}^{-}) (c (x, t) δ (x - x^{'})) . \end{aligned}

Taking the $ℓ \to 0$ limit while holding $D = ℓ^{2} d$ constant allows $\partial_{x}^{\pm}$ and $\partial_{x^{'}}^{\pm}$ to converge to true derivatives, $\partial_{x}$ and $\partial_{x^{'}}$ . Additionally, if the $δ c (x^{'}, t) / δ t$ term on the left-hand side of Equation 35 is replaced with $δ c (x^{'}, t^{'}) / δ t$ for $t^{'} \neq t$ , then the entire right-hand side must go to 0 as the system is Markovian. This can be accomplished by multiplying the right-hand side by a factor of $δ_{t, t^{'}}$ . Taking the $δ t \to 0$ limit then turns the two terms on the left-hand side into true derivatives in time, $\partial_{t}$ and $\partial_{t^{'}}$ , acting on $c (x, t)$ and $c (x^{'}, t^{'})$ respectively while the factor of $δ_{t, t^{'}} / δ t$ on the right-hand side becomes $δ (t - t^{'})$ . Altogether, this transforms Equation 35 into

⟨ \partial_{t^{'}} c (x^{'}, t^{'}) \partial_{t} c (x, t) ⟩ = 2 D δ (t - t^{'}) \partial_{x} \partial_{x^{'}} (c (x, t) δ (x - x^{'})) .

Finally, by approximating the system as being in steady state, $c (x, t)$ can be replaced with $\bar{c} (x)$ and $\partial_{t} c (x, t)$ becomes equivalent to $η_{D} (x, t)$ . Making these substitutions and generalizing Equation 36 to arbitrary dimensions yields

⟨ η_{D} ({\vec{x}}^{'}, t^{'}) η_{D} (\vec{x}, t) ⟩ = 2 D δ (t - t^{'}) \vec{\nabla} \cdot {\vec{\nabla}}^{'} (\bar{c} (\vec{x}) δ^{S P} (\vec{x} - {\vec{x}}^{'})) .

Fourier transforming Equation 37 can be easily performed due to the $δ$ functions, integrating the spatial terms by parts, and utilizing Equation 24 to yield

\begin{aligned} ⟨ {\tilde{η}}_{D}^{*} ({\vec{k}}^{'}, ω^{'}) {\tilde{η}}_{D} (\vec{k}, ω) ⟩ = \int d^{S P} x d^{S P} x^{'} d t d t^{'} e^{i \vec{k} \cdot \vec{x}} e^{- i {\vec{k}}^{'} \cdot {\vec{x}}^{'}} e^{i ω t} e^{- i ω^{'} t^{'}} ⟨ η_{D} ({\vec{x}}^{'}, t^{'}) η_{D} (\vec{x}, t) ⟩ \\ = 2 D \int d^{S P} x d^{S P} x^{'} d t d t^{'} e^{i \vec{k} \cdot \vec{x}} e^{- i {\vec{k}}^{'} \cdot {\vec{x}}^{'}} e^{i ω t} e^{- i ω^{'} t^{'}} δ (t - t^{'}) \vec{\nabla} \cdot {\vec{\nabla}}^{'} (\bar{c} (\vec{x}) δ^{S P} (\vec{x} - {\vec{x}}^{'})) \\ = 2 D (2 π δ (ω - ω^{'})) \int d^{S P} x d^{S P} x^{'} e^{i \vec{k} \cdot \vec{x}} e^{- i {\vec{k}}^{'} \cdot {\vec{x}}^{'}} \vec{\nabla} \cdot {\vec{\nabla}}^{'} (\bar{c} (\vec{x}) δ^{S P} (\vec{x} - {\vec{x}}^{'})) \\ = 2 D (2 π δ (ω - ω^{'})) \int d^{S P} x d^{S P} x^{'} \bar{c} (\vec{x}) δ^{S P} (\vec{x} - {\vec{x}}^{'}) \vec{\nabla} \cdot {\vec{\nabla}}^{'} (e^{i \vec{k} \cdot \vec{x}} e^{- i {\vec{k}}^{'} \cdot {\vec{x}}^{'}}) \\ = 2 D \vec{k} \cdot {\vec{k}}^{'} (2 π δ (ω - ω^{'})) \int d^{S P} x d^{S P} x^{'} \bar{c} (\vec{x}) δ^{S P} (\vec{x} - {\vec{x}}^{'}) e^{i \vec{k} \cdot \vec{x}} e^{- i {\vec{k}}^{'} \cdot {\vec{x}}^{'}} \\ = 2 D \vec{k} \cdot {\vec{k}}^{'} (2 π δ (ω - ω^{'})) \int d^{S P} x \bar{c} (\vec{x}) e^{i \vec{x} (\vec{k} - {\vec{k}}^{'})} \\ = 2 D \vec{k} \cdot {\vec{k}}^{'} \tilde{\bar{c}} (\vec{k} - {\vec{k}}^{'}) (2 π δ (ω - ω^{'})) \\ = \frac{2 λ^{2} \vec{k} \cdot {\vec{k}}^{'}}{1 + λ^{2} {| \vec{k} - {\vec{k}}^{'} |}^{2}} (β {(2 π)}^{S O + 1} δ (ω - ω^{'}) δ^{S O} (\vec{k} - {\vec{k}}^{'})) . \end{aligned}

Moving on to $η_{ν}$ , its correlation function must be $δ$ correlated in time and space since it is a purely local reaction and as such, at steady state, must take the form

⟨ η_{ν} ({\vec{x}}^{'}, t^{'}) η_{ν} (\vec{x}, t) ⟩ = ν \bar{c} (\vec{x}) δ (t - t^{'}) δ^{S P} (\vec{x} - {\vec{x}}^{'}) .

Fourier transforming Equation 39 is again easily performed due to the $δ$ functions and Equation 24. This yields

\begin{array}{ll} ⟨ {\tilde{η}}_{ν}^{*} ({\vec{k}}^{'}, ω^{'}) {\tilde{η}}_{ν} (\vec{k}, ω) ⟩ = \int d^{S P} x d^{S P} x^{'} d t d t^{'} e^{i \vec{k} \cdot \vec{x}} e^{- i {\vec{k}}^{'} \cdot {\vec{x}}^{'}} e^{i ω t} e^{- i ω^{'} t^{'}} ⟨ η_{ν} ({\vec{x}}^{'}, t^{'}) η_{ν} (\vec{x}, t) ⟩ \\ = ν \int d^{S P} x d^{S P} x^{'} d t d t^{'} e^{i \vec{k} \cdot \vec{x}} e^{- i {\vec{k}}^{'} \cdot {\vec{x}}^{'}} e^{i ω t} e^{- i ω^{'} t^{'}} \bar{c} (\vec{x}) δ (t - t^{'}) δ^{S P} (\vec{x} - {\vec{x}}^{'}) \\ = ν (2 π δ (ω - ω^{'})) \int d^{S P} x d^{S P} x^{'} e^{i \vec{k} \cdot \vec{x}} e^{- i {\vec{k}}^{'} \cdot {\vec{x}}^{'}} \bar{c} (\vec{x}) δ^{S P} (\vec{x} - {\vec{x}}^{'}) \\ = ν (2 π δ (ω - ω^{'})) \int d^{S P} x e^{i \vec{x} \cdot (\vec{k} - {\vec{k}}^{'})} \bar{c} (\vec{x}) \\ = ν \tilde{\bar{c}} (\vec{k} - {\vec{k}}^{'}) (2 π δ (ω - ω^{'})) \\ = \frac{1}{1 + λ^{2} {| \vec{k} - {\vec{k}}^{'} |}^{2}} (β {(2 π)}^{S O + 1} δ (ω - ω^{'}) δ^{S O} (\vec{k} - {\vec{k}}^{'})) . \end{array}

Finally, the cross spectrum of $η_{β}$ must be $δ$ correlated in $ω$ -space as well as all source dimensions of k-space since it is merely a uniform production term that does not depend on space or time. This yields

⟨ {\tilde{η}}_{β}^{*} ({\vec{k}}^{'}, ω^{'}) {\tilde{η}}_{β} (\vec{k}, ω) ⟩ = β {(2 π)}^{S O + 1} δ (ω - ω^{'}) δ^{S O} (\vec{k} - {\vec{k}}^{'}) .

Combining Equations 29, 38, 40, and 41 then yields

\begin{aligned} ⟨ {\tilde{δ c}}^{*} ({\vec{k}}^{'}, ω^{'}) \tilde{δ c} (\vec{k}, ω) ⟩ = \frac{β {(2 π)}^{S O + 1} δ (ω - ω^{'}) δ^{S O} (\vec{k} - {\vec{k}}^{'})}{ν^{2} (1 + λ^{2} {| \vec{k} |}^{2} - i \frac{ω}{ν}) (1 + λ^{2} {| {\vec{k}}^{'} |}^{2} + i \frac{ω^{'}}{ν})} \\ \cdot (\frac{2 λ^{2} \vec{k} \cdot {\vec{k}}^{'}}{1 + λ^{2} {| \vec{k} - {\vec{k}}^{'} |}^{2}} + \frac{1}{1 + λ^{2} {| \vec{k} - {\vec{k}}^{'} |}^{2}} + 1) \\ = \frac{β {(2 π)}^{S O + 1} δ (ω - ω^{'}) δ^{S O} (\vec{k} - {\vec{k}}^{'})}{ν^{2} (1 + λ^{2} {| \vec{k} |}^{2} - i \frac{ω}{ν}) (1 + λ^{2} {| {\vec{k}}^{'} |}^{2} + i \frac{ω^{'}}{ν})} \frac{2 + λ^{2} ({| \vec{k} |}^{2} + {| {\vec{k}}^{'} |}^{2})}{1 + λ^{2} {| \vec{k} - {\vec{k}}^{'} |}^{2}} \\ = \frac{\tilde{\bar{c}} (\vec{k} - {\vec{k}}^{'}) (2 π δ (ω - ω^{'})) (2 + λ^{2} ({| \vec{k} |}^{2} + {| {\vec{k}}^{'} |}^{2}))}{ν (1 + λ^{2} {| \vec{k} |}^{2} - i \frac{ω}{ν}) (1 + λ^{2} {| {\vec{k}}^{'} |}^{2} + i \frac{ω^{'}}{ν})} . \end{aligned}

We now define m as

m (\vec{x}, t) = \int_{V (a)} d^{S P} r c (\vec{x} + \vec{r}, t),

where $V (a)$ is a $S P$ -dimensional sphere with radius a. This allows the mean value of m to be written as

\begin{aligned} \bar{m} (\vec{x}) & = \int_{V (a)} d^{S P} r \bar{c} (\vec{x} + \vec{r}) = \frac{β λ^{2 - (S P - S O)}}{D} \int_{V (a)} d^{S P} r P_{S P - S O} (\frac{| \vec{x} + \vec{r} |}{λ}) \\ = \frac{β λ^{S O}}{ν} M_{S P - S O, S P} (\frac{| \vec{x} |}{λ}, \frac{a}{λ}), \end{aligned}

where

M_{N, N^{'}} (x, y) = \int_{V (y)} d^{N^{'}} u P_{N} (| \vec{x} + \vec{u} |) .

Since $P_{N} (| \vec{x} |)$ can only depend on the last N dimensions of its input vectors, the same must be true of $M_{N, N^{'}}$ . From here we define $S (\vec{x})$ as the 0-frequency limit of the cross spectrum in $ω$ -space of m. This allows the time averaged variance, $σ^{2} (x)$ to take the form

\begin{aligned} σ^{2} (x) & = \frac{S (\vec{x})}{T} = \frac{1}{T} lim_{ω \to 0} \int \frac{d ω^{'}}{2 π} ⟨ {\tilde{δ m}}^{*} (\vec{x}, ω^{'}) \tilde{δ m} (\vec{x}, ω) ⟩ \\ = \frac{1}{T} lim_{ω \to 0} \int \frac{d ω^{'}}{2 π} \int_{V (a)} d^{S P} r d^{S P} r^{'} \int \frac{d^{S P} k}{{(2 π)}^{S P}} \frac{d^{S P} k^{'}}{{(2 π)}^{S P}} e^{- i \vec{k} \cdot (\vec{x} + \vec{r})} e^{i {\vec{k}}^{'} \cdot (\vec{x} + {\vec{r}}^{'})} ⟨ {\tilde{δ c}}^{*} ({\vec{k}}^{'}, ω^{'}) \tilde{δ c} (\vec{k}, ω) ⟩ \\ = \frac{1}{{(2 π)}^{2 S P} ν T} \int_{V (a)} d^{S P} r d^{S P} r^{'} \int d^{S P} k d^{S P} k^{'} e^{- i \vec{k} \cdot (\vec{x} + \vec{r})} e^{i {\vec{k}}^{'} \cdot (\vec{x} + {\vec{r}}^{'})} \\ \cdot \frac{\tilde{\bar{c}} (\vec{k} - {\vec{k}}^{'}) (2 + λ^{2} ({| \vec{k} |}^{2} + {| {\vec{k}}^{'} |}^{2}))}{(1 + λ^{2} {| \vec{k} |}^{2}) (1 + λ^{2} {| {\vec{k}}^{'} |}^{2})} \\ = \frac{1}{{(2 π)}^{2 S P} ν T} \int_{V (a)} d^{S P} r d^{S P} r^{'} \int d^{S P} k d^{S P} k^{'} d^{S P} z e^{- i \vec{k} \cdot (\vec{x} + \vec{r})} e^{i {\vec{k}}^{'} \cdot (\vec{x} + {\vec{r}}^{'})} e^{i \vec{z} \cdot (\vec{k} - {\vec{k}}^{'})} \\ \cdot \bar{c} (\vec{z}) \frac{2 + λ^{2} ({| \vec{k} |}^{2} + {| {\vec{k}}^{'} |}^{2})}{(1 + λ^{2} {| \vec{k} |}^{2}) (1 + λ^{2} {| {\vec{k}}^{'} |}^{2})} \\ = \frac{1}{{(2 π)}^{2 S P} ν T} \int_{V (a)} d^{S P} r d^{S P} r^{'} \int d^{S P} k d^{S P} k^{'} d^{S P} z e^{- i \vec{k} \cdot (\vec{x} + \vec{r} - \vec{z})} e^{i {\vec{k}}^{'} \cdot (\vec{x} + {\vec{r}}^{'} - \vec{z})} \\ \cdot \bar{c} (\vec{z}) (\frac{1}{1 + λ^{2} {| \vec{k} |}^{2}} + \frac{1}{1 + λ^{2} {| {\vec{k}}^{'} |}^{2}}) \\ = \frac{1}{{(2 π)}^{2 S P} ν T} \int_{V (a)} d^{S P} r d^{S P} r^{'} \int d^{S P} z \bar{c} (\vec{z}) (\int d^{S P} k e^{- i \vec{k} \cdot (\vec{x} + \vec{r} - \vec{z})} \frac{{(2 π)}^{S P} δ^{S P} (\vec{x} + {\vec{r}}^{'} - \vec{z})}{1 + λ^{2} {| \vec{k} |}^{2}} \\ + \int d^{S P} k^{'} e^{i {\vec{k}}^{'} \cdot (\vec{x} + {\vec{r}}^{'} - \vec{z})} \frac{{(2 π)}^{S P} δ^{S P} (\vec{x} + \vec{r} - \vec{z})}{1 + λ^{2} {| {\vec{k}}^{'} |}^{2}}) \\ = \frac{β λ^{2 - (S P - S O)}}{D ν λ^{S P} T} \int_{V (a)} d^{S P} r d^{S P} r^{'} \int d^{S P} z P_{S P - S O} (\frac{| \vec{z} |}{λ}) \\ \cdot (δ^{S P} (\vec{x} + {\vec{r}}^{'} - \vec{z}) P_{S P} (\frac{| \vec{x} + \vec{r} - \vec{z} |}{λ}) + δ^{S P} (\vec{x} + \vec{r} - \vec{z}) P_{S P} (\frac{| \vec{x} + {\vec{r}}^{'} - \vec{z} |}{λ})) \\ = \frac{β λ^{4 - (2 S P - S O)}}{D^{2} T} \int_{V (a)} d^{S P} r d^{S P} r^{'} P_{S P} (\frac{| \vec{r} - {\vec{r}}^{'} |}{λ}) (P_{S P - S O} (\frac{| \vec{x} + \vec{r} |}{λ}) + P_{S P - S O} (\frac{| \vec{x} + {\vec{r}}^{'} |}{λ})) \\ = \frac{β λ^{4 - (S P - S O)}}{D^{2} T} (\int_{V (a)} d^{S P} r M_{S P, S P} (\frac{| \vec{r} |}{λ}, \frac{a}{λ}) P_{S P - S O} (\frac{| \vec{x} + \vec{r} |}{λ}) \\ + \int_{V (a)} d^{S P} r^{'} M_{S P, S P} (\frac{| {\vec{r}}^{'} |}{λ}, \frac{a}{λ}) P_{S P - S O} (\frac{| \vec{x} + {\vec{r}}^{'} |}{λ})) \\ = \frac{2 β λ^{4 - (S P - S O)}}{D^{2} T} \int_{V (a)} d^{S P} r M_{S P, S P} (\frac{| \vec{r} |}{λ}, \frac{a}{λ}) P_{S P - S O} (\frac{| \vec{x} + \vec{r} |}{λ}) \\ = \frac{2 β λ^{S O}}{ν^{2} T} Σ_{S P - S O, S P} (\frac{| \vec{x} |}{λ}, \frac{a}{λ}) = \frac{2 \bar{m} (\vec{x})}{ν T} \frac{Σ_{S P - S O, S P} (\frac{| \vec{x} |}{λ}, \frac{a}{λ})}{M_{S P - S O, S P} (\frac{| \vec{x} |}{λ}, \frac{a}{λ})}, \end{aligned}

where

Σ_{N, N^{'}} (x, y) = \int_{V (y)} d^{N^{'}} u M_{N^{'}, N^{'}} (u, y) P_{N} (| \vec{x} + \vec{u} |) .

Wherein once again only the last N dimensions of the input vectors can be taken into account. Combining Equations 44 and 46 yields the full precision to be

P^{2} (\vec{x}) = \frac{{\bar{m}}^{2} (\vec{x})}{σ^{2} (\vec{x})} {(\frac{Δ \bar{m} (\vec{x})}{\bar{m} (\vec{x})})}^{2} = \frac{T}{2 τ} \frac{β}{ν} M_{S P - S O, S P} (\frac{| \vec{x} |}{λ}, \frac{a}{λ}) {(1 - \frac{M_{S P - S O, S P} (\frac{| \vec{x} | + 2 a}{λ}, \frac{a}{λ})}{M_{S P - S O, S P} (\frac{| \vec{x} |}{λ}, \frac{a}{λ})})}^{2},

where

τ = \frac{1}{ν} \frac{Σ_{S P - S O, S P} (\frac{| \vec{x} |}{λ}, \frac{a}{λ})}{M_{S P - S O, S P} (\frac{| \vec{x} |}{λ}, \frac{a}{λ})} .

With Equation 48, once the forms of $P_{N}$ , $M_{N, N^{'}}$ , and $Σ_{N, N^{'}}$ are determined for a given $S P$ and $S O$ , the full form of the noise-to-signal ratio can be found. We now calculate these forms for specific choices of $S P$ and $S O$ .

1D space, 0D source

To begin, we start with the simple scenario in which $S P = 1$ and $S O = 0$ . This allows P₁, $M_{1, 1}$ , and $Σ_{1, 1}$ to take the forms

P_{1} (x) = \int \frac{d u}{2 π} e^{- i u x} \frac{1}{1 + u^{2}} = \frac{1}{2} e^{- | x |}

\begin{aligned} M_{1, 1} (x, y) & = \int_{- y}^{y} d u P_{1} (| x + u |) = \frac{1}{2} \int_{- y}^{y} d u e^{- | x + u |} \\ = {\begin{cases} 1 - e^{- y} \cosh (x) & x < y \\ e^{- x} \sinh (y) & x \geq y \end{cases} \end{aligned}

\begin{aligned} Σ_{1, 1} (x, y) & = \int_{- y}^{y} d u M_{1, 1} (u, y) P_{1} (| x + u |) = \frac{1}{2} \int_{- y}^{y} d u (1 - e^{- y} \cosh (u)) e^{- | x + u |} \\ = {\begin{cases} 1 - \frac{1}{4} e^{- y} ((5 + 2 y - e^{- 2 y}) \cosh (x) - 2 x \sinh (x)) & x < y \\ \frac{1}{4} e^{- x} (4 \sinh (y) - e^{- y} (2 y + \sinh (2 y))) & x \geq y \end{cases} \end{aligned}

Equation 51 and Equation 52 can then be put into Equation 48 along with the assumption $| x | > a$ to obtain

P^{2} (x) = \frac{\bar{m} (x) T}{2 τ} {(1 - e^{- \frac{2 a}{λ}})}^{2},

and

τ = \frac{1}{ν} (1 - e^{- \frac{a}{λ}} \frac{\frac{2 a}{λ} + \sinh (\frac{2 a}{λ})}{4 \sinh (\frac{a}{λ})})

as in Equations 9 and 10 of the main text.

Next we apply the condition given by Equation 17 to determine the regime in which these results are valid for the $S P = 1$ and $S O = 0$ case. A similar methodology can be done for each of the other cases we will look at, though this is the only one we do explicitly. To begin, we will reperform the calculation done in Equation 46 but without taking the $ω \to 0$ limit so as to obtain the full form of $S (ω, x)$ .

\begin{aligned} S (ω, x) & = \int \frac{d ω^{'}}{2 π} ⟨ {\tilde{δ m}}^{*} (x, ω^{'}) \tilde{δ m} (x, ω) ⟩ \\ = \int \frac{d ω^{'}}{2 π} \int_{- a}^{a} d r d r^{'} \int \frac{d k}{2 π} \frac{d k^{'}}{2 π} e^{- i k (x + r)} e^{i k^{'} (x + r^{'})} ⟨ {\tilde{δ c}}^{*} (k^{'}, ω^{'}) \tilde{δ c} (k, ω) ⟩ \\ = \frac{1}{{(2 π)}^{2} ν} \int_{- a}^{a} d r d r^{'} \int d k d k^{'} e^{- i k (x + r)} e^{i k^{'} (x + r^{'})} \frac{\tilde{\bar{c}} (k - k^{'}) (2 + λ^{2} (k^{2} + {k^{'}}^{2}))}{(1 + λ^{2} k^{2} - i \frac{ω}{ν}) (1 + λ^{2} {k^{'}}^{2} + i \frac{ω}{ν})} \\ = \frac{1}{{(2 π)}^{2} ν} \int_{- a}^{a} d r d r^{'} \int d k d k^{'} d z e^{- i k (x + r)} e^{i k^{'} (x + r^{'})} e^{i z (k - k^{'})} \bar{c} (z) \frac{2 + λ^{2} (k^{2} + {k^{'}}^{2})}{(1 + λ^{2} k^{2} - i \frac{ω}{ν}) (1 + λ^{2} {k^{'}}^{2} + i \frac{ω}{ν})} \\ = \frac{1}{{(2 π)}^{2} ν} \int_{- a}^{a} d r d r^{'} \int d k d k^{'} d z e^{- i k (x + r - z)} e^{i k^{'} (x + r^{'} - z)} \bar{c} (z) (\frac{1}{1 + λ^{2} k^{2} - i \frac{ω}{ν}} + \frac{1}{1 + λ^{2} {k^{'}}^{2} + i \frac{ω}{ν}}) \\ = \frac{1}{{(2 π)}^{2} ν} \int_{- a}^{a} d r d r^{'} \int d z \bar{c} (z) (\int d k e^{- i k (x + r - z)} \frac{2 π δ (x + r^{'} - z)}{1 + λ^{2} k^{2} - i \frac{ω}{ν}} + \int d k^{'} e^{i k^{'} (x + r^{'} - z)} \frac{2 π δ (x + r - z)}{1 + λ^{2} {k^{'}}^{2} + i \frac{ω}{ν}}) \\ = \frac{1}{2 π ν} \int_{- a}^{a} d r d r^{'} (\int d k e^{- i k (r - r^{'})} \frac{\bar{c} (x + r^{'})}{1 + λ^{2} k^{2} - i \frac{ω}{ν}} + \int d k e^{i k (r^{'} - r)} \frac{\bar{c} (x + r)}{1 + λ^{2} k^{2} + i \frac{ω}{ν}}) \\ = \frac{1}{ν λ} \int_{- a}^{a} d r d r^{'} (\bar{c} (x + r^{'}) Q (\frac{r - r^{'}}{λ}, - \frac{ω}{ν}) + \bar{c} (x + r) Q (\frac{r - r^{'}}{λ}, \frac{ω}{ν})), \end{aligned}

where

Q (x, y) = \int \frac{d u}{2 π} e^{- i u x} \frac{1}{1 + i y + u^{2}} = \frac{1}{2 \sqrt{1 + i y}} e^{- | x | \sqrt{1 + i y}}

and $\sqrt{1 + i y}$ is assumed to be the branch with a positive real component. Plugging this and the explicit form of $\bar{c}$ into Equation 55 then yields

\begin{aligned} S (ω, x) & = \frac{β}{4 ν D} \int_{- a}^{a} d r d r^{'} (\frac{1}{\sqrt{1 - i \frac{ω}{ν}}} e^{- \frac{| x + r^{'} |}{λ}} e^{- \frac{| r - r^{'} |}{λ} \sqrt{1 - i \frac{ω}{ν}}} + \frac{1}{\sqrt{1 + i \frac{ω}{ν}}} e^{- \frac{| x + r |}{λ}} e^{- \frac{| r - r^{'} |}{λ} \sqrt{1 + i \frac{ω}{ν}}}) \\ = \frac{β}{2 ν D} Re [\int_{- a}^{a} d r d r^{'} \frac{1}{\sqrt{1 + i \frac{ω}{ν}}} e^{- \frac{| x + r |}{λ}} e^{- \frac{| r - r^{'} |}{λ} \sqrt{1 + i \frac{ω}{ν}}}] = \frac{2 β}{ν^{2}} Re [Υ (\frac{| x |}{λ}, \frac{a}{λ}, \sqrt{1 + i \frac{ω}{ν}})], \end{aligned}

where

Υ (x, y, w) = \int_{- y}^{y} 𝑑 u 𝑑 u^{'} \frac{1}{4 w} e^{- | x + u |} e^{- w | u - u^{'} |} .

The function $Υ (x, y, w)$ limits to $Σ_{1, 1} (x, y)$ when $w \to 1$ and as such has different forms when x is less or greater than y. As the purpose of this exercise is to determine the regime in which our theoretical approximations are valid and our model obeys $| x | \geq 2 a$ for all cells, here we will only present the $x > y$ solution for simplicity. Using this to perform the integrals in Equation 58 and applying the result toEquation 57 then yields

S (ω, x) = \frac{2 β}{ν^{2}} e^{- \frac{| x |}{λ}} Re [\frac{1}{W^{2}} \sinh (\frac{a}{λ}) - e^{- \frac{a}{λ} W} \frac{W \sinh (\frac{a}{λ} W) \cosh (\frac{a}{λ}) - \cosh (\frac{a}{λ} W) \sinh (\frac{a}{λ})}{W^{2} (W^{2} - 1)}],

where $W = \sqrt{1 + i \frac{ω}{ν}}$ , which in turn implies $ω = - i ν (W^{2} - 1)$ . With this, it is easier to perform all further calculations with respect to W and take the $W \to 1$ limit as that is equivalent to the $ω \to 0$ limit.

We can now combine this with the known form of $S (0, x)$ given in Equation 46 to evaluate Equation 17 to take the form

\begin{aligned} T^{2} & ≫ | {\frac{6}{S} \frac{\partial^{2} S}{\partial ω^{2}} |}_{ω = 0} | = | {\frac{6}{S} \frac{\partial W}{\partial ω} \frac{\partial}{\partial W} (\frac{\partial W}{\partial ω} \frac{\partial S}{\partial W}) |}_{W = 1} | = | {\frac{6}{S} ({(\frac{\partial ω}{\partial W})}^{- 2} \frac{\partial^{2} S}{\partial W^{2}} - {(\frac{\partial ω}{\partial W})}^{- 3} \frac{\partial^{2} ω}{\partial W^{2}} \frac{\partial S}{\partial W}) |}_{W = 1} | \\ = \frac{1}{2 ν^{2}} \frac{96 \sinh (\frac{a}{λ}) - e^{- \frac{a}{λ}} (48 \frac{a}{λ} + 30 {(\frac{a}{λ})}^{2} + 8 {(\frac{a}{λ})}^{3} + 33 \sinh (2 \frac{a}{λ})) + 6 e^{- 3 \frac{a}{λ}} (3 \frac{a}{λ} + {(\frac{a}{λ})}^{2})}{4 \sinh (\frac{a}{λ}) - e^{- \frac{a}{λ}} (2 \frac{a}{λ} + \sinh (2 \frac{a}{λ}))} . \end{aligned}

The right-hand side of Equation 60 is a function that monotonically increases from $9 / 2 ν^{2}$ to $21 / 2 ν^{2}$ as $a / λ$ goes from 0 to $\infty$ . Thus, regardless of the value of $λ$ , $ν$ sets the timescale to which T must be compared.

2D space, 0D source

For $S P = 2$ and $S O = 0$ , P₂, $M_{2, 2}$ , and $Σ_{2, 2}$ each take the form

P_{2} (x) = \int \frac{d^{2} u}{{(2 π)}^{2}} e^{- i \vec{u} \cdot \vec{x}} \frac{1}{1 + {| \vec{u} |}^{2}} = \frac{1}{2 π} K_{0} (x)

\begin{aligned} M_{2, 2} (x, y) & = \int_{V (y)} d^{2} u P_{2} (| \vec{x} + \vec{u} |) = \int_{V (y)} d^{2} u \int \frac{d^{2} u^{'}}{{(2 π)}^{2}} e^{- i {\vec{u}}^{'} \cdot (\vec{x} + \vec{u})} \frac{1}{1 + {| {\vec{u}}^{'} |}^{2}} \\ = y \int_{0}^{\infty} d u^{'} \frac{J_{0} (x u^{'}) J_{1} (y u^{'})}{1 + {u^{'}}^{2}} \end{aligned}

\begin{aligned} Σ_{2, 2} (x, y) & = \int_{V (y)} d^{2} u M_{2, 2} (| \vec{u} |, y) P_{2} (| \vec{x} + \vec{u} |) \\ = y \int_{V (y)} d^{2} u \int_{0}^{\infty} d u^{'} \int \frac{d^{2} u^{''}}{{(2 π)}^{2}} \frac{J_{0} (| \vec{u} | u^{'}) J_{1} (y u^{'})}{1 + {u^{'}}^{2}} \frac{e^{- i {\vec{u}}^{''} \cdot (\vec{x} + \vec{u})}}{1 + {| {\vec{u}}^{''} |}^{2}} \\ = y^{2} \int_{0}^{\infty} d u^{'} d u^{''} \frac{u^{''} J_{0} (x u^{''}) J_{1} (y u^{'}) (u^{'} J_{0} (y u^{''}) J_{1} (y u^{'}) - u^{''} J_{0} (y u^{'}) J_{1} (y u^{''}))}{({u^{'}}^{2} - {u^{''}}^{2}) (1 + {u^{'}}^{2}) (1 + {u^{''}}^{2})}, \end{aligned}

where $J_{n} (x)$ and $K_{n} (x)$ are the Bessel functions of the first kind and modified Bessel functions of the second kind, respectively. Unfortunately, the complicated nature of Bessel functions makes the remaining integrals unsolvable analytically, and therefore we evaluate them numerically. Similar problems arise whenever $S P = 2$ or $S P - S O = 2$ .

3D space, 0D source

For $S P = 3$ and $S O = 0$ , P₃, $M_{3, 3}$ , and $Σ_{3, 3}$ each take the form

P_{3} (x) = \int \frac{d^{3} u}{{(2 π)}^{3}} e^{- i \vec{u} \cdot \vec{x}} \frac{1}{1 + {| \vec{u} |}^{2}} = \frac{1}{4 π x} e^{- x}

\begin{aligned} M_{3, 3} (x, y) & = \int_{V (y)} d^{3} u P_{3} (| \vec{x} + \vec{u} |) = \frac{1}{4 π} \int_{V (y)} d^{3} u \frac{1}{| \vec{x} + \vec{u} |} e^{- | \vec{x} + \vec{u} |} \\ = {\begin{cases} 1 - \frac{1 + y}{x} e^{- y} \sinh (x) & x < y \\ \frac{1}{x} e^{- x} (y \cosh (y) - \sinh (y)) & x \geq y \end{cases} \end{aligned}

\begin{aligned} Σ_{3, 3} (x, y) & = \int_{V (y)} d^{3} u M_{3, 3} (| \vec{u} |, y) P_{3} (| \vec{x} + \vec{u} |) \\ = \frac{1}{4 π} \int_{V (y)} d^{3} u (1 - \frac{1 + y}{| \vec{u} |} e^{- y} \sinh (| \vec{u} |)) \frac{1}{| \vec{x} + \vec{u} |} e^{- | \vec{x} + \vec{u} |} \\ = {\begin{cases} 1 - \frac{1}{4 x} e^{- y} (1 + y) ((5 + 2 y + e^{- 2 y}) \sinh (x) - 2 x \cosh (x)) & x < y \\ \frac{1}{4 x} e^{- x} (4 (y \cosh (y) - \sinh (y)) + e^{- y} (1 + y) (2 y - \sinh (2 y))) & x \geq y \end{cases} \end{aligned}

2D space, 1D source

For $S P = 2$ , $S O = 1$ , P₁ and $M_{2, 2}$ are known from Equations 50 and Equations 62. This leaves $M_{1, 2}$ and $Σ_{1, 2}$ to take the forms

\begin{aligned} M_{1, 2} (x, y) & = \int_{V (y)} d^{2} u P_{1} (| \vec{x} + \vec{u} |) = \frac{1}{2} \int_{0}^{y} d u \int_{0}^{2 π} d θ u e^{- | x_{2} + u_{2} |} \\ = e^{- | x_{2} |} \int_{0}^{2 π} d θ \frac{1 - e^{- y \sin (θ)} (1 + y \sin (θ))}{2 {(\sin (θ))}^{2}} \end{aligned}

\begin{array}{ll} Σ_{1, 2} (x, y) & = \int_{V (y)} d^{2} u M_{2, 2} (| \vec{u} |, y) P_{1} (| \vec{x} + \vec{u} |) \\ = \frac{y}{2} \int_{0}^{y} d u \int_{0}^{2 π} d θ \int_{0}^{\infty} d u^{'} u \frac{J_{0} (u u^{'}) J_{1} (y u^{'})}{1 + {u^{'}}^{2}} e^{- | x_{2} + u \sin (θ) |} \end{array}

Again, we evaluate the remaining integrals numerically.

3D space, 2D source

For $S P = 3$ , $S O = 2$ , P₁ and $M_{3, 3}$ are known from Equations 50 and Equations 65. This leaves $M_{1, 3}$ and $Σ_{1, 3}$ to take the forms

\begin{array}{ll} M_{1, 3} (x, y) & = \int_{V (y)} d^{3} u P_{1} (| \vec{x} + \vec{u} |) = \frac{1}{2} \int_{V (y)} d^{3} u e^{- | x_{3} + u_{3} |} \\ = 2 π {\begin{cases} e^{- y} (1 + y) \cosh (x) + \frac{y^{2} - x^{2}}{2} - 1 & x < y \\ e^{- x} (y \cosh (y) - \sinh (y)) & x \geq y \end{cases} \end{array}

\begin{array}{ll} Σ_{1, 3} (x, y) & = \int_{V (y)} d^{3} u M_{3, 3} (| \vec{u} |, y) P_{1} (| \vec{x} + \vec{u} |) = \frac{1}{2} \int_{V (y)} d^{3} u (1 - \frac{1 + y}{| \vec{u} |} e^{- y} \sinh (| \vec{u} |)) e^{- | x_{3} + u_{3} |} \\ = 2 π {\begin{cases} e^{- y} (1 + y) (\frac{7 + 2 y + e^{- 2 y}}{4} \cosh (x) - \frac{x}{2} \sinh (x) - \cosh (y)) + \frac{y^{2} - x^{2}}{2} - 1 & x < y \\ e^{- x} (\frac{4 y^{2} + 5 y - 1}{8} e^{- y} + \frac{1 + y}{8} e^{- 3 y} + \frac{3 y}{4} \cosh (y) - \frac{5}{4} \sinh (y)) & x \geq y \end{cases} \end{array}

Appendix 2

Hopping model for SDC case

To obtain a more intuitive understanding of why the SDC model results in the scaling properties seen in the various calculations of $M_{SP - SO, SP}$ and $Σ_{SP - SO, SP}$ , we now look at a simpler version of one dimensional diffusion in which we discretize space into compartments of uniform size. Let molecules still be produced in the 0th compartment at rate $β$ and degrade anywhere in space at rate $ν$ . The process of diffusion can be approximated by letting the molecules hop to neighboring compartments with rate h with equal probability of moving left or right. This allows the dynamics of m_j, the number of molecules in the j compartment for $j \in ℤ$ , to be written as

\frac{\partial m_{j}}{\partial t} = β δ_{0 j} + h (m_{j + 1} + m_{j - 1} - 2 m_{j}) - ν m_{j} .

By setting the left-hand side of Equation 45 to 0, the resulting system of equations can be easily solved by assuming ${\bar{m}}_{j} = A exp (- 2 | j | / λ)$ and calculating A and $λ$ . Imposing this assumption on Equation 45 and taking $j > 0$ yields

\begin{aligned} 0 & = h (A e^{- \frac{2 (j + 1)}{λ}} + A e^{- \frac{2 (j - 1)}{λ}} - 2 A e^{- \frac{2 j}{λ}}) - ν A e^{- \frac{2 j}{λ}} = A e^{- \frac{2 j}{λ}} (h e^{- \frac{2}{λ}} + h e^{\frac{2}{λ}} - 2 h - ν) \\ = A e^{- \frac{2 j}{λ}} (4 h \sinh^{2} (\frac{1}{λ}) - ν) ⟹ λ = {asinh}^{- 1} (\sqrt{\frac{ν}{4 h}}) . \end{aligned}

With $λ$ solved for, we solve for the proportionality constant by noting that the total number of molecules in the whole system must follow a simple birth-death process with a mean of $β / ν$ . This in turn implies

\begin{aligned} \frac{β}{ν} & = \sum_{j = - \infty}^{\infty} A e^{- \frac{2 | j |}{λ}} = A (2 (\sum_{j = 0}^{\infty} e^{- \frac{2 j}{λ}}) - 1) = A (\frac{2}{1 - e^{- \frac{2}{λ}}} - 1) = A (\frac{e^{\frac{1}{λ}}}{\sinh (\frac{1}{λ})} - 1) = A \coth (\frac{1}{λ}) \\ ⟹ A = \frac{β}{ν} \tanh (\frac{1}{λ}), \end{aligned}

This in turn gives the average value of m_j to be

{\bar{m}}_{j} = \frac{β}{ν} \tanh (\frac{1}{λ}) e^{- \frac{2 | j |}{λ}} .

Next, we calculate the full distribution of m_j by assuming that at any given moment in time each molecule in the system has probability $P_{j}$ of being in the jth compartment. This can be combined with the aforementioned fact that N, the total number of molecules in the system, must follow a birth-death process and thus to Poissonianly distributed with mean $β / ν$ . For any given value of N, $P (m_{j} | N)$ must be a binomial distribution with success probability $P_{j}$ since each molecule is independent. This allows the marginal distribution $P (m_{j})$ to be calculated to be

\begin{aligned} P (m_{j}) & = \sum_{N = m_{j}}^{\infty} P (N) P (m_{j} | N) = \sum_{N = m_{j}}^{\infty} e^{- \frac{β}{ν}} \frac{{(\frac{β}{ν})}^{N}}{N!} (\binom{N}{m_{j}}) P_{j}^{m_{j}} {(1 - P_{j})}^{N - m_{j}} \\ = e^{- \frac{β}{ν}} \frac{{(\frac{β}{ν} P_{j})}^{m_{j}}}{m_{j}!} \sum_{N = m_{j}}^{\infty} \frac{{(\frac{β}{ν} (1 - P_{j}))}^{N - m_{j}}}{(N - m_{j})!} = e^{- \frac{β}{ν}} \frac{{(\frac{β}{ν} P_{j})}^{m_{j}}}{m_{j}!} e^{\frac{β}{ν} (1 - P_{j})} = e^{- \frac{β}{ν} P_{j}} \frac{{(\frac{β}{ν} P_{j})}^{m_{j}}}{m_{j}!} . \end{aligned}

Thus, m_j is seen to be Poissonianly distributed with mean $β P_{j} / ν$ . Comparing this mean to that derived in Equation 46 then implies

P_{j} = \tanh (\frac{1}{λ}) e^{- \frac{2 | j |}{λ}} .

We now consider the joint distribution of m_j and m_k for $j \neq k$ . Since molecules cannot be in the jth and kth compartment simultaneously, the joint conditional distribution $P (m_{j}, m_{k} | N)$ must be trinomially distributed. This allows for the joint distribution to be calculated in a manner similar to Equation 75 to produce

\begin{aligned} P (m_{j}, m_{k}) & = \sum_{N = m_{j} + m_{k}}^{\infty} P (N) P (m_{j}, m_{k} | N) \\ = \sum_{N = m_{j} + m_{k}}^{\infty} e^{- \frac{β}{ν}} \frac{{(\frac{β}{ν})}^{N}}{N!} (\binom{N}{m_{j}, m_{k}}) P_{j}^{m_{j}} P_{k}^{m_{k}} {(1 - P_{j} - P_{k})}^{N - m_{j} - m_{k}} \\ = e^{- \frac{β}{ν}} \frac{{(\frac{β}{ν} P_{j})}^{m_{j}}}{m_{j}!} \frac{{(\frac{β}{ν} P_{k})}^{m_{k}}}{m_{k}!} \sum_{N = m_{j} + m_{k}}^{\infty} \frac{{(\frac{β}{ν} (1 - P_{j} - P_{k}))}^{N - m_{j} - m_{k}}}{(N - m_{j} - m_{k})!} \\ = e^{- \frac{β}{ν}} \frac{{(\frac{β}{ν} P_{j})}^{m_{j}}}{m_{j}!} \frac{{(\frac{β}{ν} P_{k})}^{m_{k}}}{m_{k}!} e^{\frac{β}{ν} (1 - P_{j} - P_{k})} = (e^{- \frac{β}{ν} P_{j}} \frac{{(\frac{β}{ν} P_{j})}^{m_{j}}}{m_{j}!}) (e^{- \frac{β}{ν} P_{k}} \frac{{(\frac{β}{ν} P_{k})}^{m_{k}}}{m_{k}!}) . \end{aligned}

Thus, the joint probability distribution of m_j and m_k is seen to be separable into the product of the two marginal distribution, meaning that same-time, instantaneous measurements of m_j an m_k must be uncorrelated.

From here we can begin to calculate the full correlation function for m_j and m_k. We start by defining $δ m_{j} (t) = m_{j} (t) - {\bar{m}}_{j}$ and $δ m_{k} (t) = m_{k} (t) - {\bar{m}}_{k}$ . Since ${\bar{m}}_{j}$ is known to set the right-hand side of Equation 45 to 0, the dynamics of $δ m_{j}$ can be written as

\frac{\partial δ m_{j}}{\partial t} = h (δ m_{j + 1} + δ m_{j - 1} - 2 δ m_{j}) - ν δ m_{j},

with the same being true for $δ m_{k}$ . Additionally, we assume the system is at steady state so that all mean expressions are invariant to time translation. Given this, we can without loss of generality take the correlation function between $δ m_{j}$ and $δ m_{k}$ to have the form

C_{j, k} (t) = ⟨ δ m_{k} (t) δ m_{j} (0) ⟩,

where $t > 0$ . Applying the dynamic result given in Equation 78 then yields

\begin{aligned} \frac{\partial C_{j, k}}{\partial t} & = ⟨ \frac{\partial δ m_{k} (t)}{\partial t} δ m_{j} (0) ⟩ = ⟨ (h (δ m_{k + 1} (t) + δ m_{k - 1} (t) - 2 δ m_{k} (t)) - ν δ m_{k} (t)) δ m_{j} (0) ⟩ \\ = h (C_{j, k + 1} + C_{j, k - 1}) - (2 h + ν) C_{j, k} . \end{aligned}

The final form of Equation 80 can be split into the term $- (2 h + ν) C_{j, k}$ which implies $C_{j, k} \propto exp (- (2 h + ν) t)$ and the term $h (C_{j, k + 1} + C_{j, k - 1})$ which is the recursion relation for $I_{ℓ} (2 h t)$ , the modified Bessel function of the first kind, where $ℓ$ is some function of j and k. This means $C_{j, k} (t)$ can be written as

C_{j, k} (t) = A I_{ℓ (j, k)} (2 h t) e^{- (2 h + ν) t},

for some proportionality constant A.

To determine the forms of A and $ℓ (j, k)$ , we can utilize the initial condition that m_j is Poissonianlly distributed and thus has a variance equal to its mean while being completely uncorrelated with m_k when both are measured at the same time. This means $C_{j, k} (0)$ can be written as

C_{j, k} (0) = \frac{β}{ν} P_{j} δ_{j k},

which in turn implies $ℓ (j, j) = 0$ as $I_{n} (0) = δ_{0 n}$ for $n \in ℤ$ . To satisfy the recursion relation term of Equation 80, it must then be the case that $ℓ (j, j + n) = n$ . Setting $k = j + n$ thus yields $ℓ (j, k) = k - j$ . Since k and j are integers, $ℓ (j, k) = j - k$ is equally valid as $I_{n} = I_{- n}$ again for $n \in ℤ$ . Combining these results together yields the final form of $C_{j, k} (t)$ to be

C_{j, k} (t) = \frac{β}{ν} P_{j} I_{k - j} (2 h t) e^{- (2 h + ν) t} .

Next, let $τ$ be the autocorrelation time of m_j. This quantity is typically defined by integrating $C_{j, j} (t) / C_{j, j} (0)$ over all time. Using the known properties of modified Bessel functions, this can be solved to yield

τ = \int_{0}^{\infty} 𝑑 t \frac{C_{j, j} (t)}{C_{j, j} (0)} = \int_{0}^{\infty} 𝑑 t I_{0} (2 h t) e^{- (2 h + ν) t} = \frac{1}{\sqrt{ν (4 h + ν)}} .

If we now define $M = T / τ$ where M is the number of effectively independent measurements that can be made in a time T, we see that for $h ≫ ν$ , $M \approx 2 \sqrt{ν h} T$ . Additionally, from Equation 72 we see that in the $h ≫ ν$ regime $λ \approx 2 \sqrt{h / ν}$ . By equating this $λ$ to the nondimensionalized $λ_{SDC} / a$ from the SDC model we see that $M \approx λ ν T = (λ_{SDC} / a) ν T$ . This is consistent with the fact that for $λ_{SDC} ≫ a$ the right-hand side of Equation 54 becomes approximately $a / λ_{SDC} ν$ , which allows $M = T / τ_{S D C} \approx (λ_{SDC} / a) ν T$ .

In the $h ≪ ν$ regime, we find $M \approx ν T$ . Once again, this consistent with Equation 54 when $λ_{SDC} ≪ a$ as this causes the right-hand side to become approximately $ν^{- 1}$ . Thus, the SDC model is seen to have a correlation time that agrees with Equation 84 in both the large and small h regime.

Appendix 3

Comparison to experimental data

To compare our theory to experimental data, we focus on ten of the morphogens presented in Table 1 of Kicheva et al., 2012 and obtain data from the references therein. For Bicoid, we obtain a value of $λ$ of ∼100 µm from the text of Gregor et al., 2007b with and error of ±10 µm from the finding in Gregor et al., 2007a that cells have a ∼10% error in measuring the Bicoid gradient. We then take the a value of the Drosophila embryo cells that are subjected to the Bicoid gradient to be ∼2.8 µm based on Figure 3A of Gregor et al., 2007a. We use the same figure to estimate the size of the whole embryo to be ∼500 µm or ∼90 cells. This value of a is also used for Dorsal as measurements of both Bicoid and Dorsal occur in the Drosophila embryo at nuclear cycle 14. For the value of $λ$ for Dorsal, we use Figure 3D from Liberman et al., 2009 to obtain a full width at 60% max of 45 ± 10 µm. Since this represents the width of Gaussian fit on both sides of the source whereas our model uses an exponential profile, we assume the appropriate $λ$ value for such an exponential fit would be half this value, 22.5 ± 5 µm. Figure 3A from the same source also shows that the distance from the ventral midline to the dorsal midline is ∼200 µm or ∼35 cells.

For Dpp and Wg, Kicheva et al., 2007 provides explicit measurements of $λ$ for each. These values are 20.2 ± 5.7 µm and 5.8 ± 2.04 µm respectively. For Hh, we use Figure S2C in the supplementary material of Wartlick et al., 2011 to determine $λ$ to be 8 ± 3 µm. Dpp, Wg, and Hh all occur in the wing disc during the third instar of the Drosophila development. As such, we use a common value of a for all three. This value is taken to be 1.3 µm based on the area of the cells being reported as 5.5 ± 0.8 µm² in the supplementary material of Kicheva et al., 2007 and the assumption that the cells are circular. Additionally, the scale bar for Figure 1A in Wartlick et al., 2011 shows the maximal distance from the morphogen producing midline of the wing disc to its edge to be ∼250 µm or ∼100 cells.

The $λ$ value of Fgf8 is reported as being 197 ± 7 µm in Yu et al., 2009. Additionally, based off the scale bars seen in Figure 2C–E of Yu et al., 2009, we estimate the value of a for the cells to be ∼10 µm. For the morphogens involved in the Nodal/Lefty system (cyclops, squint, lefty1, and lefty2), measurements of $λ$ for each are taken from Figure 2C–F of Müller et al., 2012 by observing where the average of the three curves crosses the 37% of max threshold with error bars given by the width of the region in which the vertical error bars of each plot intersect this threshold line. We assume the a value of each morphogen in the Nodal/Lefty system to be equivalent to the a value of cells in the Fgf8 measurements performed in Yu et al., 2009. This is because the measurements made in Müller et al., 2012 were taken during the blastula stage of the zebrafish development while measurements taken in Yu et al., 2009 we taken in the sphere germ ring stage. These stages occur at ∼2.25 and ∼5.67 hpf respectively, but the blastula stage can last until ∼6 hpf based on the timeline of zebrafish development presented in Kimmel et al., 1995. As such, since there is potential overlap in the time frame of these two stages, we assume the cells maintain a relatively fixed size and thus that the value of a for the Nodal/Lefty system can be taken as the same value of a used for Fgf8. Additionally, as seen in Figures 8F and 11B in Kimmel et al., 1995, these two stages also share a rougly equal overall diameter of the embryo of ∼500 µm at the largest point. This creates a circumference of ∼1600 µm or ∼80 cells, which in turn means the morphogen must travel a maximum distance of ∼40 cells away from the source.

Morphogen	Organism	$λ$ (µm)	a (µm)	N
Bicoid	Drosophila	100 ± 10	2.8	90
Fgf8	Zebrafish	197 ± 7	10	40
Lefty2	Zebrafish	150 ± 25	10	40
Lefty1	Zebrafish	115 ± 20	10	40
Dpp	Drosophila	20.2 ± 5.7	1.3	100
Dorsal	Drosophila	22.5 ± 5	2.8	35
Squint	Zebrafish	65 ± 10	10	40
Cyclops	Zebrafish	30 ± 5	10	40
Hh	Drosophila	8 ± 3	1.3	100
Wg	Drosophila	5.8 ± 2.04	1.3	100

Data availability

All data used in this study is simulated via computational methods outlined in the manuscript and appendices.

References

1. Akiyama T
2. Gibson MC
(2015) Morphogen transport: theoretical and experimental controversies
Wiley Interdisciplinary Reviews: Developmental Biology 4:99–112.
https://doi.org/10.1002/wdev.167
- PubMed
- Google Scholar
(2011) Formation of morphogen gradients: local accumulation time
Physical Review E 83:051906.
https://doi.org/10.1103/PhysRevE.83.051906
- PubMed
- Google Scholar
1. Berg HC
2. Purcell EM
(1977) Physics of chemoreception
Biophysical Journal 20:193–219.
https://doi.org/10.1016/S0006-3495(77)85544-6
- PubMed
- Google Scholar
(2013) Cytonemes are required for the establishment of a normal Hedgehog morphogen gradient in Drosophila epithelia
Nature Cell Biology 15:1269–1281.
https://doi.org/10.1038/ncb2856
- PubMed
- Google Scholar
1. Bressloff PC
2. Kim H
(2018) Bidirectional transport model of morphogen gradient formation via cytonemes
Physical Biology 15:026010.
https://doi.org/10.1088/1478-3975/aaa64c
- PubMed
- Google Scholar
(2017) A facilitated diffusion mechanism establishes the Drosophila dorsal gradient
Development 144:4450–4461.
https://doi.org/10.1242/dev.155549
- PubMed
- Google Scholar
1. Chen W
2. Huang H
3. Hatori R
4. Kornberg TB
(2017) Essential basal cytonemes take up hedgehog in the Drosophila wing imaginal disc
Development 144:3134–3144.
https://doi.org/10.1242/dev.149856
- PubMed
- Google Scholar
1. de Lachapelle AM
2. Bergmann S
(2010) Precision and scaling in morphogen gradient read-out
Molecular Systems Biology 6:351.
https://doi.org/10.1038/msb.2010.7
- PubMed
- Google Scholar
1. Driever W
2. Nüsslein-Volhard C
(1988) The bicoid protein determines position in the Drosophila embryo in a concentration-dependent manner
Cell 54:95–104.
https://doi.org/10.1016/0092-8674(88)90183-3
- PubMed
- Google Scholar
1. Dubuis JO
2. Tkacik G
3. Wieschaus EF
4. Gregor T
5. Bialek W
(2013) Positional information, in bits
PNAS 110:16301–16308.
https://doi.org/10.1073/pnas.1315642110
- PubMed
- Google Scholar
(2014) Characterizing transport through a crowded environment with different obstacle sizes
The Journal of Chemical Physics 140:054108.
https://doi.org/10.1063/1.4864000
- PubMed
- Google Scholar
(2009) Role of spatial averaging in the precision of gene expression patterns
Physical Review Letters 103:258101.
https://doi.org/10.1103/PhysRevLett.103.258101
- PubMed
- Google Scholar
Software
1. Fancher S
(2020) morphogen_profile, version 1f1aada
Github.

https://github.com/scfancher/morphogen_profile
1. Fancher S
2. Mugler A
(2017) Fundamental limits to collective concentration sensing in cell populations
Physical Review Letters 118:078101.
https://doi.org/10.1103/PhysRevLett.118.078101
- PubMed
- Google Scholar
1. Fanelli D
2. McKane AJ
(2010) Diffusion in a crowded environment
Physical Review E 82:021113.
https://doi.org/10.1103/PhysRevE.82.021113
- Google Scholar
Book
1. Gardiner CW
(2004) Handbook of Stochastic Methods for Physics, Chemistry, and the Natural Sciences
New York: Springer.
https://doi.org/10.1002/bbpc.19850890629
- Google Scholar
1. Gierer A
2. Meinhardt H
(1972) A theory of biological pattern formation
Kybernetik 12:30–39.
https://doi.org/10.1007/BF00289234
- PubMed
- Google Scholar
1. Gillespie DT
(2000) The chemical Langevin equation
The Journal of Chemical Physics 113:297–306.
https://doi.org/10.1063/1.481811
- Google Scholar
(2007a) Probing the limits to positional information
Cell 130:153–164.
https://doi.org/10.1016/j.cell.2007.05.025
- PubMed
- Google Scholar
(2007b) Stability and nuclear dynamics of the bicoid morphogen gradient
Cell 130:141–152.
https://doi.org/10.1016/j.cell.2007.05.026
- PubMed
- Google Scholar
(2010) Modelling the bicoid gradient
Development 137:2253–2264.
https://doi.org/10.1242/dev.032409
- PubMed
- Google Scholar
(2002) Establishment of developmental precision and proportions in the early Drosophila embryo
Nature 415:798–802.
https://doi.org/10.1038/415798a
- PubMed
- Google Scholar
1. Huang H
2. Kornberg TB
(2015) Myoblast cytonemes mediate wg signaling from the wing imaginal disc and Delta-Notch signaling to the air sac primordium
eLife 4:e06114.
https://doi.org/10.7554/eLife.06114
- PubMed
- Google Scholar
(2009) Dynamics of the dorsal morphogen gradient
PNAS 106:21707–21712.
https://doi.org/10.1073/pnas.0912395106
- PubMed
- Google Scholar
(2007) Kinetics of morphogen gradient formation
Science 315:521–525.
https://doi.org/10.1126/science.1135774
- PubMed
- Google Scholar
(2012) Investigating the principles of morphogen gradient formation: from tissues to cells
Current Opinion in Genetics & Development 22:527–532.
https://doi.org/10.1016/j.gde.2012.08.004
- PubMed
- Google Scholar
(1995) Stages of embryonic development of the zebrafish
Developmental Dynamics 203:253–310.
https://doi.org/10.1002/aja.1002030302
- PubMed
- Google Scholar
1. Kornberg TB
(2014) Cytonemes and the dispersion of morphogens
Wiley Interdisciplinary Reviews: Developmental Biology 3:445–463.
https://doi.org/10.1002/wdev.151
- PubMed
- Google Scholar
1. Kornberg TB
2. Roy S
(2014) Cytonemes as specialized signaling filopodia
Development 141:729–736.
https://doi.org/10.1242/dev.086223
- PubMed
- Google Scholar
1. Lander AD
2. Nie Q
3. Wan FY
(2002) Do morphogen gradients arise by diffusion?
Developmental Cell 2:785–796.
https://doi.org/10.1016/S1534-5807(02)00179-X
- PubMed
- Google Scholar
(2009) Quantitative imaging of the dorsal nuclear gradient reveals limitations to threshold-dependent patterning in Drosophila
PNAS 106:22317–22322.
https://doi.org/10.1073/pnas.0906227106
- PubMed
- Google Scholar
1. Müller P
2. Rogers KW
3. Jordan BM
4. Lee JS
5. Robson D
6. Ramanathan S
7. Schier AF
(2012) Differential diffusivity of nodal and lefty underlies a reaction-diffusion patterning system
Science 336:721–724.
https://doi.org/10.1126/science.1221920
- PubMed
- Google Scholar
1. Müller P
2. Rogers KW
3. Yu SR
4. Brand M
5. Schier AF
(2013) Morphogen transport
Development 140:1621–1638.
https://doi.org/10.1242/dev.083519
- PubMed
- Google Scholar
1. Rogers KW
2. Müller P
(2019) Nodal and BMP dispersal during early zebrafish development
Developmental Biology 447:14–23.
https://doi.org/10.1016/j.ydbio.2018.04.002
- PubMed
- Google Scholar
1. Rogers KW
2. Schier AF
(2011) Morphogen gradients: from generation to interpretation
Annual Review of Cell and Developmental Biology 27:377–407.
https://doi.org/10.1146/annurev-cellbio-092910-154148
- PubMed
- Google Scholar
(2013) Specialized filopodia direct long-range transport of SHH during vertebrate tissue patterning
Nature 497:628–632.
https://doi.org/10.1038/nature12157
- PubMed
- Google Scholar
1. Shvartsman SY
2. Baker RE
(2012) Mathematical models of morphogen gradients and their effects on gene expression
Wiley Interdisciplinary Reviews: Developmental Biology 1:715–730.
https://doi.org/10.1002/wdev.55
- PubMed
- Google Scholar
1. Stanganello E
2. Scholpp S
(2016) Role of cytonemes in Wnt transport
Journal of Cell Science 129:665–672.
https://doi.org/10.1242/jcs.182469
- PubMed
- Google Scholar
1. Teimouri H
2. Kolomeisky AB
(2016a) New model for understanding mechanisms of biological signaling: direct transport via cytonemes
The Journal of Physical Chemistry Letters 7:180–185.
https://doi.org/10.1021/acs.jpclett.5b02703
- PubMed
- Google Scholar
1. Teimouri H
2. Kolomeisky AB
(2016b) Mechanisms of the formation of biological signaling profiles
Journal of Physics A: Mathematical and Theoretical 49:483001.
https://doi.org/10.1088/1751-8113/49/48/483001
- Google Scholar
(2007) Fundamental limits to position determination by concentration gradients
PLOS Computational Biology 3:e78.
https://doi.org/10.1371/journal.pcbi.0030078
- PubMed
- Google Scholar
1. Varennes J
2. Fancher S
3. Han B
4. Mugler A
(2017) Emergent versus Individual-Based multicellular chemotaxis
Physical Review Letters 119:188101.
https://doi.org/10.1103/PhysRevLett.119.188101
- PubMed
- Google Scholar
(2011) Dynamics of Dpp signaling and proliferation control
Science 331:1154–1159.
https://doi.org/10.1126/science.1200037
- PubMed
- Google Scholar
(2017) Control of signaling molecule range during developmental patterning
Cellular and Molecular Life Sciences 74:1937–1956.
https://doi.org/10.1007/s00018-016-2433-5
- PubMed
- Google Scholar
1. Yu SR
2. Burkhardt M
3. Nowak M
4. Ries J
5. Petrásek Z
6. Scholpp S
7. Schwille P
8. Brand M
(2009) Fgf8 morphogen gradient forms by a source-sink mechanism with freely diffusing molecules
Nature 461:533–536.
https://doi.org/10.1038/nature08391
- PubMed
- Google Scholar

Article and author information

Author details

Sean Fancher
1. Department of Physics and Astronomy, Purdue University, West Lafayette, United States
2. Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, United States
Contribution
Conceptualization, Data curation, Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

For correspondence
sfancher@sas.upenn.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-8701-192X
Andrew Mugler
1. Department of Physics and Astronomy, Purdue University, West Lafayette, United States
2. Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, United States
Contribution
Conceptualization, Supervision, Funding acquisition, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing

For correspondence
andrew.mugler@pitt.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-9367-7026

Funding

Simons Foundation (376198)

Sean Fancher
Andrew Mugler

Simons Foundation (568888)

Sean Fancher

National Science Foundation (PHY-1945018)

Andrew Mugler

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 supported by Simons Foundation Grants No. 376198 and 568888 and National Science Foundation Grant No. PHY-1945018. We thank Chris Bairnsfather for useful discussions.

Version history

Received: May 15, 2020
Accepted: October 13, 2020
Accepted Manuscript published: October 14, 2020 (version 1)
Version of Record published: November 4, 2020 (version 2)
Version of Record updated: November 6, 2020 (version 3)

Copyright

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.