Superdiffusion of morphogens by receptor-mediated transport

Binding to cell surface receptors is thought to be essential for the transport of certain morphogens in developing tissues. The finite number of receptors per cell turns the tissue into a subdiffusive medium for the morphogens. We study a simple microscopic model of receptor-mediated transport and find superdiffusive spreading of morphogens. We propose that the superdiffusive spreading in a subdiffusive medium is due to a ratchet effect. A phenomenological model within the framework of the fractional Fokker–Planck equation allows us to analytically study the formation of morphogen gradients. Within this model, we show furthermore that the same features leading to the anomalous transport behavior also result in gradients that are robust against changes in the morphogen secretion rate. Together these findings show that anomalous transport in biological systems can be intimately linked to essential biological features.


Introduction
The development of patterns in biological systems is a long-standing problem. In developing organisms, positional information is often provided by morphogens [1]. These are molecules that affect the fate of a cell in a concentration dependent manner: depending on whether the morphogen concentration in the environment is below or above a certain threshold value, the further development of a cell will be different. During development, morphogen molecules are secreted from a localized source and then spread into the adjacent tissue. As morphogens are taken up by the tissue cells, a spatial gradient in the morphogen concentration forms: in steady state it decreases with the distance from the source. In this way morphogens lead to a patterning of the developing organism.
Even though the concept of morphogens is almost 40 years old by now, the processes underlying morphogen transport are not very well understood. For a long time it was thought that morphogens simply diffuse in the extracellular space between cells [2]. Making use of GFP-fusions to morphogens, this question has recently been re-addressed [3]- [7]. These studies revealed that morphogen transport is more involved than simple diffusion. In particular, there is evidence that morphogen binding to cell surface receptors plays an essential role for transport. In particular, experiments suggest that, for the morphogen decapentaplegic (DPP) in the developing wing of the fruit fly Drosophila melanogaster, transcytosis is the dominant form of transport [4,5]. In this process morphogens bound to surface receptors are internalized into cells, reappear at different position on the cell surface, and are subsequently released into extracellular space.
A theoretical investigation of this transport process has revealed interesting consequences of this form of transport [8,9]. In particular, it was shown that the morphogen gradient does not significantly change if the secretion rate of morphogens from the source is changed as long as it stays beyond a certain value. This feature, known as robustness [10], is linked to there being a finite number of receptors, which limits the maximal current of morphogens through cells. If the morphogens outnumber the receptors, a part of the molecules get trapped in extracellular space and are not carried on.
Robustness is also a feature of other transport processes [10,11]. For example, it was shown that this property emerges if morphogens are immobilized for some time in extracellular space. Recently, Hornung et al [12] suggested that such a process leads to anomalous diffusion. They analyzed the consequences of entrapping with a long-tailed distribution of trapping times within the frame of continuous time random walks (CTRW). This leads to subdiffusive transport and implies robustness. The number of free molecules in the extracellular space between cells n and n + 1 is L n , the number of molecules bound to surface receptors of cell n is S n . Binding to a free receptor occurs at rate k on and unbinding from an occupied receptor at rate k off . Surface receptors are assumed to be very mobile, such that exchange of receptors between the two surfaces is fast and can be adiabatically eliminated. As a consequence morphogen release occurs with equal probability on either side of a cell.
In this work, we investigate a simplified model for morphogen transport mediated by cell surface receptors. We find long-tailed waiting time distributions for morphogens in the extracellular space. In striking contrast to the mechanism presented in [12], however, we find that receptor-mediated transport implies superdiffusive spreading of an initially localized distribution. To reveal the mechanism generating superdiffusive transport, we propose a phenomenological model of receptor-mediated transport in terms of a CTRW. By analyzing the corresponding fractional Fokker-Planck equation (FFPE), we calculate the stationary state and obtain a simple expression for the robustness of the stationary profile with respect to changes in the morphogen secretion rate.

Superdiffusion in the absence of degradation
We start by introducing a model describing morphogen transport via cells mediated by carriers. For simplicity, we will consider a one-dimensional (1D) system where all cells are arranged in a line. They are represented as sites on a lattice, see figure 1. The number of morphogens between cells n and n + 1 is denoted by L n , the number of receptors on cell n occupied by morphogens is denoted by S n . Each cell is assumed to contain the same number R of receptors. The dynamics of these quantities is then given bẏ where n = . . . , −2, −1, 0, 1, 2, . . .. Furthermore, k on is the binding rate of extracellular morphogens to unoccupied surface receptors of an adjacent cell and k off is the rate at which morphogens bound to receptors detach. Morphogen degradation occurs with a rate b deg and is restricted to receptor-bound molecules. This is motivated by the well-established lysosomal  (1) and (2). We find n 2 ∼ t 1.1 , i.e. superdiffusive spreading. Inset: concentration L tot n at the initial site n = 0 and at site n = 5. We find L tot 0 ∼ t −0.59 which corresponds to a long-tailed waiting time distribution, equation (3). In a CTRW this would lead to subdiffusive behavior. Parameters are k on R/k off = 10 −4 and there are 10 5 particles in the system. pathway for degradation of morphogens, e.g. in embryos of the fruit-fly Drosophila [13]. For simplicity, we neglect extracellular degradation in the above equations. If small, taking it into account does not qualitatively change our results. The dynamic equations (1) and (2) are obtained from the model of transcytosis introduced in [8] in the limit of infinite fast exchange between receptors on the cell surface and in the cell To study the transport generated by the dynamic equations, we first investigate the spreading of ligands in the absence of degradation, b deg = 0. We assume that initially all ligands are located on site n = 0 between two cells and that the lattice is infinite. If the initial number of morphogens is small then the mean square displacement of the total morphogen distribution L tot n = L n + (S n + S n+1 )/2 increases linearly with time, i.e. the distribution spreads diffusively. If, however, the number of morphogens is sufficiently large, then transport is superdiffusive, n 2 ∼ t β with β = 1.1 for k on R/k off = 10 −4 and 10 5 particles in the system, see figure 2. Only after the number of morphogens has decreased everywhere to sufficiently small values, transport crosses over to being diffusive. The larger the number of morphogens, the longer the transport is superdiffusive.
In figure 3, we present the values of β as a function of k on R/k off and for different numbers of particles in the system. For a fixed value of k on R/k off , the value of β increases with the number of particles in the system, while for fixed particle number, the value increases with the value of k on R/k off . While for 1000 particles in the system, transport is diffusive, β = 1, for 10 4 and 10 5 particles we clearly find superdiffusive spreading.
The origin of anomalous transport apparently lies in the finite number of receptors present on the surface of each cell. In the presence of a large number of morphogens, the primary effect Exponents characterizing the spreading of an initially localized particle distribution as a function of system parameters. We have n 2 ∼ t β and L tot 0 ∼ t −α . The inset displays the ratio β/α. Number of particles in the system are, respectively, 1000 (black diamonds), 10 4 (blue boxes), and 10 5 (red circles).
of having a finite number of carriers, however, is to hamper transport, which should lead to a slowing down of morphogen spreading. Indeed, if all receptors adjacent to an extracellular site are occupied then transport of further unbound morphogens is blocked. That this is the case can be seen in the decay of the morphogen number on extracellular sites which follows a power law L 0 ∼ t −α with 0 < α < 1, see figure 2 (inset) and figure 3.
Although the inhibition of transport through the non-availability of empty receptor sites is a genuine non-linear effect, some insight can be obtained by considering a linear phenomenological description of the dynamics generated by equations (1) and (2). To this end, we interpret the algebraic decay of the number of morphogens mentioned above to result from a long-tailed waiting-time distribution w(t) with asymptotic power law behavior Then, the spreading of morphogens by surface receptors can be described by a CTRW similarly to [12] where binding of morphogens to extracellular heparan sulfate proteoglycans was assumed to result in a long-tailed waiting-time distribution. The time-evolution of the morphogen distribution can then be described by the FFPE [14] ∂ t λ(x, t) where λ(x, t) is the density of morphogens at the space-time point (x, t) and L FP ≡ −K α ∂ 2 x is the Fokker-Planck operator with the generalized diffusion coefficient K α . For the dynamics specified by (1) and (2), K α = a 2 /τ α , where a is the distance between the two extracellular sites, i.e. the step size of the random walker, and τ a characteristic time, which is related to the binding and unbinding constants k on and k off as well as to the number of receptors and morphogens. Due to the latter, τ is different for every site and changes with time. For simplicity, we will consider in the following only the case with τ = const. More details of the connection between the CTRW and the FFPE are presented in the appendix.

The process in equation (4) differs from the usual diffusion equation by the action in time of the Riemann-Liouville operator
see, e.g. [15], which reflects non-Markovian properties of random walks through algebraic distributions of waiting times between steps of the walker [14,16]. The waiting-time distribution (3) with 0 < α < 1 corresponds to subdiffusive transport with mean-square displacement x(t) 2 = K α t α , see [14,17]. This is at odds with the observed superdiffusive behavior obtained from numerical integration of the dynamic equations (1) and (2). The key to resolving this dilemma lies in realizing that blocking due to occupied receptors not only prevents the forward, but also the backward transport: consider a highly occupied extracellular site, say n = 0, with empty receptors on the neighboring cells. Once all receptors are occupied, further binding of free morphogens at n = 0 is blocked. Only if bound morphogens are released into the extracellular space can further free morphogens bind. As there are many more morphogens at n = 0, the newly binding morphogens come with an overwhelming probability from this site rather than form n = ±1. The occupation of receptors thus leads to a kind of ratchet effect.
In the phenomenological FFPE, this effect can be taken into account by adding a convective current in the direction opposite to concentration gradients. For a distribution initially localized at x = 0, we write where sign(x) = 1 if x > 0 and −1 for x < 0 and v α = vτ 1−α . Here, we have used the Galilean variant form of the fractional diffusion-advection equation to assure that the distributionlike the solutions of equations (1) and (2)-stays unimodal for all times, see [14]. As for the generalized diffusion constant K α and in order to illustrate the consequences of the ratchet effect, we assume the convection velocity v to be constant. In principle, it depends on the morphogen distribution and thus on time and space. Note that due to symmetry, x = 0 for all times. Fractional diffusion-advection equations have been studied in various contexts and can lead to superdiffusive transport [14], [18]- [23]. This also holds for equation (6). Multiplication of equation (6) with x 2 and integration with respect to x yields d dt where Z (t) = ∞ 0 xλ(x, t) dx. A dynamic equation for Z is obtained by multiplication of equation (6) with x and integration with respect to x from 0 to ∞. It reads Obviously, λ(x = 0, t) decays with time and can for long times be neglected compared to 1/2. From the definition of the fractional derivative, see equation (5), one gets and 7 see, e.g. [15]. In addition, from equation (8) one obtains Z (t) ∼ v α t α /2 (α + 1), such that equation (7) leads to Therefore, if α > 1/2, transport is superdiffusive with a transport exponent 2α > 1. Note that from the numerical solutions presented in figure 2 we obtain α = 0.59, while the transport exponent is β = 1.1 ≈ 2α = 1.18 supporting our phenomenological approach to superdiffusion by receptor-mediated transport. The inset in figure 3 shows that the linear phenomenological description breaks down, however, as α approaches 1. Here, the blocking mechanism is not of subdiffusive nature anymore and the non-linear effects in the dynamics are no longer captured by our phenomenological linear model. Fundamentally, however, the mechanism underlying superdiffusion is also in this case a combination of inhibition of transport and the ratchet effect.

Robustness
So far, we have neglected the degradation of morphogens. For the formation of morphogen gradients in developing tissues, though, it is essential. We now describe how this process can be incorporated into the phenomenological framework of equation (6). In the microscopic dynamic equations (1) and (2) morphogens decay at rate b deg while they are bound to receptors. Correspondingly, in the phenomenological framework of the CTRW introduced above, degradation only occurs during particle transport. We therefore phenomenologically describe degradation by a sink term of the form The Riemann-Liouville operator in front of the usual degradation term restricts degradation to occur only simultaneously to transport. For the special case α = 1/2, this form of the degradation term has been derived in [24] in the context of the comb model [23,25]. We now investigate the spreading of morphogens that are secreted from a source at x = 0 and that spread into the adjacent tissue according to the fractional diffusion-advection equation (6) with sink term (12). Without loss of generality only the half-space x > 0 will be considered. The effect of the source can then be captured by the boundary condition at x = 0. The rate of morphogen production can be expressed by a current j 0 into the system. However, due to the blocking of free morphogens by occupied receptors, only some of the molecules introduced into the system will be transported. We therefore choose as boundary condition Note that the exact form of the boundary condition will not be important in the following. In the present framework, it is reasonable to assume that the current of morphogens through transport by carriers is only a fraction of the morphogens secreted per unit time into the system. Stationary solutions to the problem are obtained by setting ∂ t λ = 0 which leads to = 0 if f (t) ≡ 0 this equation can be easily solved yielding 8 The characteristic decay length a/B α of the morphogen profile is determined by the phenomenological parameters α, τ and v: An important feature of biological systems is their robustness against external perturbations. Morphogen gradients generated by transcytosis have been found to be robust against changes in the morphogen secretion rate [8]. In order to investigate robustness in the framework of the fractional diffusion-advection equation (6), we calculate the robustness R as defined in [8], R( j 0 ) = a( j 0 ∂ j 0 x(λ; j 0 )) −1 , where λ is the stationary profile (14) and x(λ; j 0 ) is obtained by inverting equation (14). Note that since equation (13) is linear, the value of R is the same for all values of x. It measures the shift of the stationary profile along x as the current j 0 is changed by a factor 2. If R > 1 then the morphogen profile shifts by less than a cell diameter and the pattern determined by the gradient is robust. For the distribution (14) we find Note that the appearance of α in the denominator of R is a consequence of the chosen boundary condition at the source x = 0. For different choices, the dependence is less important if not absent. The value of B α and hence of R can be increased by increasing b deg τ . For τ → ∞, we find B α = ab deg /v.

Discussion
In summary, we have shown that receptor-mediated transport of morphogens is superdiffusive if the number of morphogens is much larger than the number of receptors. In that case, free morphogens get trapped in the extracellular space, turning the tissue into a subdiffusive medium with a waiting time distribution t −1−α . But since the blocking of transport by occupied receptors induces an advective current through a ratchet effect, superdiffusive spreading of morphogens can occur as we have shown by studying a phenomenological CTRW model. Our numerical results support the phenomenological approach which gives 2α for the transport exponent. Finally, based on the corresponding fractional advection-diffusion equation, the robustness of morphogen gradients to changes in the morphogen secretion rate was calculated. Note that while we haven chosen here to present our results for the simple model given by (1) and (2), the same phenomena can also be observed in the full transcytosis model of Bollenbach et al [8].
They really depend only on the finite number of carriers and are therefore a generic feature of receptor-mediated transport processes. It might seem puzzling that the above waiting time distribution, which implies a diverging mean first passage time, goes along with superdiffusive transport. Intuitively, the superdiffusive growth of the mean square displacement can be traced back to the fact that some particles are trapped for a long time at the initial location, while others get advected through the ratchet effect. So both features of the transport, subdiffusive waiting time and advection, act together to produce the superdiffusive behavior.
The phenomenological parameters of the CTRW are, of course, linked to the parameters of the microscopic model, equations (1) and (2). Numerics shows that the exponents α and β increase with larger values of k on R/k off . They also increase with increasing particle numbers. For low enough values of k on R/k off and small enough particle numbers, transport tends to become subdiffusive. The value of τ is expected to increase with κ = k on (k off R) −1 . For large values of κ binding to receptors is strongly hampered. Consequently, the characteristic time 9 needed to be transferred by a cell is larger. This is also in agreement with the analysis of robustness in [8], where it was shown that robustness increases with κ. The detailed relation between the phenomenological and the microscopic parameters still requires further studies.