Diffusion vs. direct transport in the precision of morphogen readout

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.

Within developing organisms, morphogen profiles provide cells with information about their position relative to other cells. While it has been shown that cells can use this information to determine their position with extremely high precision [1][2][3][4][5], the mechanism by which morphogen profiles are formed is still not well agreed upon. One well-known theory for such a 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 [6][7][8][9][10][11]. Alternatively, a direct transport (DT) model has also been proposed where source cells create protrusions called cytonemes through which morphogen molecules can travel and be delivered directly to the target cells [6,9,[11][12][13]. The presence of these two competing 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 [4,14]. The accumulation times for several morphogens in Drosophila have been measured and found to match the predictions made by the SDC model [15]. In zebrafish, the molecular dynamics of the morphogen Fgf8 have been measured and found to be consistent with Brownian diffusion through extracellular space [16]. 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 has been seen to be highly correlated in both space and time with the formation of cytonemes [17], while Wnt morphogens have been found to be highly localized around cell protrusions such as cytonemes [18,19]. Theoretical studies of both the SDC * Electronic address: amugler@purdue.edu and DT models have examined these measurable effects [12,15,[20][21][22], 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 characteristic length scale of the morphogen profile.
Here we derive fundamental limits to the precision of morphogen concentration sensing for both the SDC and DT models. We focus on the stochastic noise caused by production, transport, and degradation of the morphogen molecules in the steady state regime. Several past studies have focused on the dynamics of morphogen profiles and their accumulation times [12,15,[20][21][22]. Here, we model morphogen profiles in the steady state regime as many of the experimental measurements to which we will later compare our results were taken during stages when the steady state approximation is valid [16,[23][24][25][26].
Intuitively one might expect the DT model to have less noise due to the fact that molecules are directly deposited and cannot be lost to the surrounding environment. Surprisingly, we show below that for sufficiently large lengths of the morphogen profile, the SDC model produces less noise due to it being able achieve a higher effective independent measurement count. This result holds for one-, two-, and three-dimensional geometries. We compare our results with measurements of the above morphogens and find quantitative consistency with our predictions, suggesting that readout precision plays an important role in determining the mechanisms of morphogen profile establishment.
We begin by considering a simple version of the DT model in which there is a single source cell that produces morphogen at rate β (Fig. 1A). Morphogen molecules are then transported to the jth target cell through the cytonemes, a dynamic process which we assume can be well approximated in steady state by a constant rate γ j . The morphogen then degrades within the target cells at rate ν. Letting N be the total number of target cells, the dynamics of the number of morphogen molecules in the source cell m 0 (t) and the number in the jth target cell m j (t) are The first two terms on the righthand side of Eq. 1 represent the production of morphogen and stochastic noise inherent within that process. The remaining two terms represent the transport of morphogen to the target cells and its noise. The same two terms also appear as the first terms on the righthand side of Eq. 2, while the last two terms represent the morphogen degradation and its noise. The noise terms obey η β t η β (t) = βδ t − t , η γ,k t η γ,j (t) = γ jm0 δ jk δ t − t , and η ν,k t η ν,j (t) = νm j δ jk δ t − t respectively [27].
From here, we assume that each cell integrates its morphogen molecule count over a time T [28] such that δm j 2 is the variance in the time average Since Eqs. 1 and 2 are linear with Gaussian white noise, solving for δm j 2 is straightforward: we Fourier transform Eqs. 1 and 2 in space and time, calculate the power spectrum of m j , and approximate δm j 2 by its low-frequency limit [29]. This lowfrequency approximation assumes that T {τ 1 , τ 2 }, where τ 1 = Γ −1 and τ 2 = ν −1 are the characteristic timescales of the morphogen molecules set by their transport and degradation rates respectively. For the jth target cell, this procedure yields a relative error of The two factors in Eq. 3 have an intuitive interpretation: the relative error is seen to decrease with both the mean number of moleculesm j and the number of independent measurements T /τ 2 made in time T , where τ 2 is the autocorrelation time. One important aspect of Eq. 3 is that it is identical to the case in which there is no extrinsic noise inherited from m 0 . This is because the cytoneme transport process is modelled as a simple transfer action. Since no molecules are lost, any noise m j inherits from m 0 is cancelled by the negative cross correlations between the two [29]. The net result is that the precision of morphogen readout in any target cell is equivalent to that of a simple Poisson process, with no extra noise inherited from production and transport. We now consider the SDC model (Fig. 1B). As most morphogen systems are modelled in one-dimensional space within the SDC model [15,20,22], we will begin with this assumption, although we also consider other geometries later on. Consider a system with a single source cell located at the origin which produces morphogen molecules at a rate β. These molecules then diffuse freely through space with diffusion coefficient D to produce a density field of morphogen, ρ (x, t). The morphogen also spontaneously degrades at any point in space with rate ν. The dynamics of ρ are The first two terms on the righthand side of Eq. 4 represent the morphogen diffusion and noise inherent within the diffusion process. The third and fourth terms represent the degradation and its noise, and the last terms represent the production and its noise isolated to the location of the source cell. The noise terms obey βδ t − t respectively [30][31][32]. Hereρ (x) is the steady state mean concentration of morphogen molecules at position x, which for a one-dimensional system with a single source takes the formρ (x) = βλexp −|x| /λ /2D with λ = D/ν. We then imagine a target cell located at x that is perfectly permeable to the morphogen and counts the number m (x, t) = V dy ρ (x + y, t) of morphogen molecules within its volume V as a measure of its local concentration. We use this simpler prescription over explicitly accounting for more complicated structures such as surface receptor binding because it has been shown that it ultimately yields similar concentration sensing results up to a factor of order unity [28]. Assuming |x| ≥ 2a with a being the cell radius enforces the condition that there is no overlap between the source and target cell and produces a mean value ofm (x) = β sinh a/λ exp −|x| /λ /ν.

Once again defining δm (x)
2 as the variance of the time average of m (x, t) requires T τ 3 with τ 3 = ν + D/a 2 −1 to make the low-frequency approximation applicable. With this, we calculate δm (x) 2 to be [29] δm (x) m (x) whereλ = λ/a. Unlike in Eq. 3, Eq. 5 has no negative cross correlation terms to cancel the extrinsic noise. Instead, the noise is comprised of a combination of three positive terms coming from the production, diffusion, and degradation noise terms seen in Eq. 4 [29]. These terms combine to produce the compact form in Eq. 5. Eq. 5 has the same prefactor seen in Eq. 3 but has an additional factor that is necessarily less than 1. This additional factor is due to the fact that morphogen can diffuse away from the cell as well as degrade, which in turn means νT is an underestimate of the number of effectively independent measurements. The additional factor compensates for this underestimation.
By comparing Eqs. 3 and 5, it is possible to determine which model achieves a lower variance and, in turn, higher precision. We take an example system in which there are N target cells in a line such that the position of the jth target cell is x j = 2aj. This creates a line of N + 1 adjacent cells with the first one being the source of the morphogen, as seen in Fig. 1. By setting γ j = γ 0 exp −2j/λ for any given basal transport rate γ 0 , the mean m values of both the SDC and DT models can be made to scale identically. These means can be made precisely equal to each other by appropriate choices of β and ν in each system. We assume the β values in both models are equivalent, as production from the source cell is a mechanism common to both models. We then fix the ν values to achieve equality in the means. With these choices made, we define R to be the ratio of Eq. 5 to Eq. 3, which yields When R is less (greater) than 1, the SDC (DT) model achieves lesser error and is thus the more precise method of morphogen sensing. Fig. 2A shows the boundary in (N ,λ) space that separates these two regimes (solid curve). We also calculate this boundary for a threedimensional system assuming a plane of source cells (dotdashed) and numerically solve for this boundary in a two-dimensional system assuming a line of source cells (dashed), both of which also produce exponentially decaying mean morphogen profiles [29].  [16,24,25,[33][34][35][36]; see [29] for further discussion on experimental data.
As can be seen in Fig. 2A, for sufficiently large population size N and any value ofλ larger than ∼2.35, the SDC model is more precise. This is surprising because unlike the SDC model, the DT model carries no extrinsic noise. Additionally, in the SDC model morphogen molecules diffuse into all of space and may be "lost" in the sense that they could require effectively infinite amounts of time to reach the target cell or, in the case of 3D space, may simply never reach it at all. These effects might suggest that the DT model should be more precise. However, our results contradict this intuition for sufficiently large values of N andλ.
As suggested by Eq. 5, we can understand this result by considering how quickly the morphogen molecules in a target cell are refreshed. In the case of the DT model, morphogen is cleared with rate ν, meaning that in a time T the cell can make ∼ ν DT T independent measurements of its morphogen concentration. In the SDC model, morphogen can be cleared via degradation at rate ν or by simply diffusing away with an effective rate D/a 2 , meaning that a cell can make ∼ ν SDC + D/a 2 T independent measurements. We define R eff to be the ratio of the number of measurements made in the DT model to that made in the SDC model. Using the same parameter relations as in Eq. 6, including the ratio of ν values between the two models, yields (7) As with R, when R eff is less (greater) than 1, the SDC (DT) model achieves a larger number of independent measurements and would thus be more likely to achieve a higher precision. The boundary separating these regimes is plotted in Fig. 2A (dotted) and is seen to be similar to the boundary separating the R regimes. As can also be seen in Fig. 2A, each boundary converges to its own critical value ofλ, which we denote asλ c , as N increases. Since this convergence happens quickly, we focus on the N → ∞ case and ignore systems with a size on the order of only a few cells. This allows us to make the prediction that any exponential morphogen profile in whichλ >λ c should use the SDC model over the DT model as it provides more precise concentration sensing, and vice versa forλ <λ c .
We now test this prediction for various morphogens. In Drosophila, the morphogen Wingless has been shown to be localized near cell protrusions such as cytonemes [18,19], while the Hedgehog gradient has been observed to correlate very highly in both space and time with the formation of cytonemes [17]. These results suggest that these two morphogen profiles could be formed via a method akin to the DT model presented here, and as seen in Fig. 2B, they both have relatively low values of λ [25,33,34] and exist near the regime in which we predict the DT model is more precise than the SDC model (yellow circles). Conversely, Bicoid has been used as a model example of SDC for decades and shown to provide cells with high levels of precision in their positional information [3,4,14]. This is consistent with the fact that the Bicoid gradient has also been measured to have a very largeλ value [24,33], putting it far into the regime where we predict the SDC model achieves a higher precision than the DT model (blue circle). The gradient sizes of Dorsal and Dpp have also been measured [25,33,35] and are less deeply within in the regime where we predict the SDC model is more precise (white circles), although for Dpp there is evidence for a variety of different gradient formation mechanisms [6,9,11].
In zebrafish, the morphogen Fgf8 has been studied at the single molecule level and found to have molecular dynamics closely matching those expected of Brownian movement and SDC model [16]. The gradient size of Fgf8 has also been measured [16,33] and found to be deeply in our predicted SDC regime, as seen in Fig. 2B (rightmost blue square). Similarly, Cyclops, Squint, Lefty1, and Lefty2, all of which are involved in the Nodal/Lefty system, have been shown to spread diffusively and be able to affect cells distant from their source [9,37] (other squares). This would support a gradient formation mechanism similar to the SDC model although measurements of the size of Cyclops and Squint gradients [33,36] put them closer to the regime where we predict that the DT model should be more precise (white squares). This is consistent with the fact that 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 [7,37].
Taken together, the results in Fig. 2B support our prediction that profiles with short length scales should form via DT, whereas profiles with long length scales should form via SDC. Morphogens with experimental evidence for DT fall within a factor of two or three from our predicted DT regime, which is reasonable given the simplicity of our models, while morphogens with experimental evidence for SDC fall squarely within our predicted SDC regime.
We have shown that in the steady-state regime, the SDC and DT 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 inherits no extrinsic noise from the morphogen production in the source cell and molecules are never lost to diffusion, the ability of molecules to diffuse into and away from a target cell in the SDC model allows the cell to make a greater number of effectively independent measurements in the same time frame. By examining how these phenomena affect the cells' sensory precision, we predicted that morphogen profiles with lengths shorter than a critical value 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 largely 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 morphogen gradient sensing to further biological contexts.

Appendix A: Direct Transfer (DT)
Assume there are N + 1 cells, one source cell and N target cells. The source cell produces morphogen molecules at rate β. These molecules can then be sent to the jth target cell at rate γ j , wherein they degrade with rate ν. Let m 0 be the number of molecules in the source cell and m j be the number of molecules in the jth target cell. Approximating these values as continuous allows for the equations where the various η terms are Gaussian white noise terms. For simplification, let Γ = j γ j . With that, linearizing and separating Eqs. A1 and A2 into m j =m j + δm j yields Setting Eqs. A3a and A4a to 0 yields mean values for m 0 and m j to bē From here, Eqs. A3b and A4b can be Fourier transformed to yield Since each η term represents the noise from a simple linear reaction, they must be independent of each other, and their cross spectrums must be δ-correlated in time and frequency space with magnitudes equal to the mean propensity of their respective reactions. This implies These allow for the cross spectrum of m 0 to be Utilizing the zero-frequency approximation here allows for the noise-to-signal ratio for m 0 to be Moving on to m j , Eqs. A8-A11 allow for the cross spectrum between m j and m k to be Once again utilizing the zero-frequency approximation allows for the noise-to-signal ratio for m j to be Eq. A15 gives the full form the noise-to-signal ratio for the jth target cell. Interestingly, it contains no cross correlation terms with the source cell or any other target cell nor any term inherited from the production noise in the source cell. The reason for this can be seen midway through the derivation of Eq. A14. The noise terms that explicitly arise from the production, degredation, and transport processes are highlighted as well as the positive cross correlation term that arises from production noise and negative cross correlation term that arises from transport noise. The negative transport correlation is seen to cancel both the production noise and production correlation, thus leaving the cell with only a degredation noise term and an identical noise term from the transport process.

Appendix B: Synthesis-Diffusion-Clearance (SDC)
We can now compare this model to the SDC model. 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 ρ which must follow where SP is the number of spatial dimensions, SO is the dimensionality of the source, and ∇ 2 is taken over all SP dimensions. Of important note is that δ SP −SO ( x) is a δ function only in the last SP − SO dimensions of the space. So, for example, if there was a 1 dimensional source in 3 dimensional space, then δ 3−1 ( x) would be a δ function in theŷ andẑ directions but not thex direction. This means that β and η β will have units of T −1 L −SO , where T is time and L is space.
We can now assume ρ has reached a steady state and separate it into ρ =ρ + δρ, which in turn allows Eq. B1 to separate into Fourier transforming Eq. B2 in space and dividing it by ν then yields Of similarly important note is that δ SO k is a δ function only in the first SO dimensions of k-space. So in the 1 dimensional source, 3 dimensional space example δ SO k would be a δ function in thex direction of k-space but not theŷ orẑ directions. This allowsρ to be written as where 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 1 dimensional source, 3 dimensional space example, P 3−1 | x| /λ should only take the y and z components of 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, Eq. B3 can be Fourier transformed in space and time to yield where η β k, ω depends only on the first SO dimensions of k-space. Assuming the η terms are all independent of each other allows the cross spectrum of ρ to be δ ρ * k , ω δ ρ k, ω = 1 The cross spectrum of η D can be obtained from its correlation function, which is known from [30] to take the form Fourier transforming Eq. B10 can be easily performed due to the δ functions, integrating the spatial terms by parts, and utilizing Eq. B4 to yield Moving on to η ν , its correlation function is known to take to form Fourier transforming Eq. B12 is again easily performed due to the δ functions and Eq. B4. This yields 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 ofρ. This yields Combining Eqs. B9, B11, B13, and B14 then yields We now define m as where V (a) is a SP -dimensional sphere with radius a. This allows the mean value of m to be written as where Since P N | 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 ( x) as the 0-frequency limit of the cross spectrum in ω-space of m. This allows it to take the form where Wherein once again only the last N dimensions of the input vectors can be taken into account. Combining Eqs. B17 and B19 yields the full 0-frequency noise-to-signal ratio of m to be With Eq. B21, once the forms of P N , M N,N , and Σ N,N are determined for a given SP and SO, the full form of the noise-to-signal ratio can be found.
Now assume that there are N target cells placed in a line such that x j = 2aj. This allows R 1,1 {x j }, N to take the form If the production rates are scaled such that β c = β d then there exists a critical value ofλ = λ/a atλ c ≈ 2.346900061. For anyλ <λ c , R 1,1 must be greater than 1, which in turn implies the DT model must produce lesser noise regardless of population size. The R 1,1 = 1 curve in (N,λ) space with β c = β d is plotted in Fig. 2A (1D case) of the main text.
Eqs. D2 and D3 can also be put into Eq. B21 to obtain the noise-to-signal ratio of the SDC model. However, these expressions do not make clear the distinction between the components that arise from the three distinct noise terms seen in Eq. B9. We now rederive Eq. B21 specifically for the SP = 1 and SO = 0 geometry as an example of how these terms can be held separate. To do so, we go back to Eq. B15 and insert it into the second line of Eq. B19 along with Eqs. B11, B13, and B14 to obtain where drdr dkdk e −ik(x+r) e ik (x+r ) (1 + λ 2 k 2 ) 1 + λ 2 k 2 Utilizing Eq. B4 allows for I D (x) to be solved, yielding Eqs. B6 and D1 can then be used to obtain The same technique can be used to solve for I ν (x) to yeild Combining Eqs. B17, B21, D5, and D10-D12 gives the full form of the SP = 1, SO = 0 noise-to-signal ratio.
Each I (x) term represents the noise inherited from a different source. Importantly, each different I (x) term is nonnegative for any value of x. This shows that unlike Eq. A15 in the DT model, there are no negative correlation terms to cancel any of the inherited noise terms.
Appendix E: 2D space, 0D source For SP = 2 and SO = 0, P 2 , M 2,2 , and Σ 2,2 each take the form 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. Similar problems arise whenever SP = 2 or SP − SO = 2. Again scaling the production rate such that β c = β d , one critical value ofλ can be recovered atλ c ≈ 18.85244498. For anyλ <λ c , R 3,3 must necessarily be greater than 1, which in turn implies the DT model must produce lesser noise. There is also a lower bound on N at N c = 5. For any N ≤ N c , R 3,3 must necessarily be greater than 1, which again implies the DT model must produce lesser noise.
Appendix G: 2D space, 1D source For SP = 2, SO = 1, P 1 and M 2,2 are known from Eqs. D1 and E2. This leaves M 1,2 and Σ 1,2 to take the forms Unfortunately, the remaining integrals are unsolvable analytically. They can, however, be calculated numerically, and by assuming the source line in the SDC model is a line of cells, the production rate of each individual cell in the line can be made equal to that of the source cell in the DT model by setting β c = 2aβ d . The numerically calculated R 1,2 = 1 curve in (N,λ) space with β c = 2aβ d is plotted in Fig. 2A (2D case) of the main text.
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 [16]. This is because the measurements made in [36] we taken during the blastula stage of the zebrafish development while measurements taken in [16] 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 [? ]. 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. Theλ = λ/a value of these ten morphogens can be seen plotted in Fig. 2B of the main text.