# Application of Fractional Calculus to Reaction-Subdiffusion Processes and Morphogen Gradient Formation

###### Abstract

It is a well known fact that subdiffusion equations in terms of fractional derivatives can be obtained from Continuous Time Random Walk (CTRW) models with long-tailed waiting time distributions. Over the last years various authors have shown that extensions of such CTRW models incorporating reactive processes to the mesoscopic transport equations may lead to non-intuitive reaction-subdiffusion equations. In particular, one such equation has been recently derived for a subdiffusive random walker subject to a linear (first-order) death process. We take this equation as a starting point to study the developmental biology key problem of morphogen gradient formation, both for the uniform case where the morphogen degradation rate coefficient (reactivity) is constant and for the non-uniform case (position-dependent reactivity). In the uniform case we obtain exponentially decreasing stationary concentration profiles and we study their robustness with respect to perturbations in the incoming morphogen flux. In the non-uniform case we find a rich phenomenology at the level of the stationary profiles. We conclude that the analytic form of the long-time morphogen concentration profiles is very sensitive to the spatial dependence of the reactivity and the specific value of the anomalous diffusion coefficient.

###### pacs:

02.50.Ey, 82.39.-k, 82.40.-g, 82.33.-z## I Introduction

Fractional diffusion equations are a powerful tool to study anomalous transport processes, i.e., processes in which the mean square displacement of a randomly moving particle displays the long time-behavior , where is the anomalous diffusion exponent and is the so-called anomalous diffusion coefficient. When , one has sublinear growth of (subdiffusion), while for one speaks of superdiffusion. As it is well known, the classical diffusion equation corresponding to the case can be obtained from an average over the trajectories of a Markovian random walk in the limit of large time scales and long displacements. In contrast, stochastic transport processes governed by anomalous diffusion equations reflect memory effects at the microscopic level. In particular, one can show that a suitably defined non-Markovian hopping process, namely a Continuous Time Random Walk (CTRW) with a long-tailed waiting time distribution, yields a subdiffusion equation in terms of the Riemann-Liouville fractional derivative MetzlerKlafterPhysReport . This fractional subdiffusion equation can be taken as a starting point to deal with a number of biologically relevant problems, e.g. the localization of a target protein by a sea of subdiffusively moving ligands in the intracellular environment targetyuste ; target09 . In this case, the complexity of the cell medium results in the ligands encountering a large number of obstacles, barriers, etc. in the course of their trajectories. In the framework of CTRW models, the effect of this crowded environment can be partly captured using waiting time distributions; subsequent averaging of the resulting equations over trajectories then leads to the associated fractional subdiffusion equations.

While anomalous diffusion and in particular subdiffusive processes play a central role in Nature as a manifestation of underlying memory effects at a microscopic level, the situation where the particles simultaneously undergo anomalous transport and reaction (understood as a particle creation, destruction or transformation process) is also very common and important from the point of view of biological applications. In the example of the target protein and the ligands given above, one could allow e.g. for the possibility of the ligands undergoing a degradation process as they sweep the cell medium. Degradation implies a change in chemical structure which results in the ligands losing their ability to interact with the target; for practical purposes this kind of transformation can therefore be regarded as a “death” or “evanescence” process.

In what follows, we shall focus on yet another biological process where degradation/death plays a central role, namely morphogen gradient formation. The location, differentiation and fate of many embryonic cells is governed by the spatial distribution of special signaling molecules called morphogens. Standard models of morphogen gradient formation assume that a specific part of the embryo secrets morphogens at a constant rate. The secreted morphogens then undergo degradation as they disseminate through the tissue and a concentration gradient builds up. Different target genes in the embryonic cells are activated above different morphogen concentration thresholds, implying that the cell response to the local environment will depend on how large the concentration is. Thanks to this differential response, cells are able to interpret the morphogen gradient and translate it into specific “code” for their further development via the expression of the relevant genes.

Traditional models of morphogen gradient formation are based on classical diffusion equations with a linear degradation term. Here, we aim to go one step further and allow for the possibility of anomalous transport, as memory effects are likely to strongly influence the stochastic motion of morphogens in the complex embryonic environment. A remarkable property of morphogen gradients is their robustness against changes or fluctuations in the rate of morphogen production or degradation. An interesting question which we aim to study is the interplay between subdiffusion and robustness with respect to such perturbations.

If one accepts the idea that morphogens perform subdiffusive motion as a result of strong dispersion in their waiting times between consecutive jumps, great caution must be exercised when incorporating the effect of a simultaneous degradation process to the transport equations because of the non-Markovian character of the latter. Several recent works indeed illustrate that heuristic equations where one has separate terms for the reaction and the transport process may lead to unphysical results, e.g. negative particle concentrations (see e.g. Langlands ). Therefore, the derivation of physically correct (but not necessarily intuitive) reaction-subdiffusion equations calls for the use of an extended CTRW formalism where the effect of reaction is incorporated at a mesoscopic level of description. In a recent work Nuestropaper the authors have shown by means of Fourier-Laplace techniques that CTRW models extended in such a way yield equations which (in addition to a standard, purely reactive term) display a mixed reaction-transport term containing both the reaction rate coefficient (reactivity) and a Riemann-Liouville fractional derivative with respect to time.

The remainder of the paper is organized as follows. We first give a brief reminder of classical reaction diffusion equations used for modeling of morphogen gradients and subsequently discuss how to extend such equations to account for anomalous transport via fractional derivatives. We subsequently focus on the specific case of uniform reactivity and assess the robustness of the resulting stationary concentration profiles with respect to perturbations of the incoming morphogen flux. Next, we turn to the non-uniform case and discuss the long-time behaviour of the profiles for several specific situations, namely the case of a piecewise constant reactivity (which not always leads to a stationary profile) and the case of a decaying reactivity respectively given by an exponential and a power law. In some cases, we also provide numerical simulation results based on a CTRW model for evanescent particles and find excellent agreement with the analytic results obtained from the fractional diffusion equation approach. Finally, a summary of results and possible avenues for future research in this area are given in the conclusions section.

## Ii Classical reaction-diffusion equation with linear degradation

The cornerstone of many studies concerning morphogen gradients is the classical one-dimensional reaction-diffusion equation

(1) |

where the evolution of the concentration is described by a Fickian term (characterized by a classical diffusion coefficient ) and a linear degradation term (characterized by the reactivity ). Eq. (1) is then solved subject to the radiation-type boundary condition

(2) |

This boundary condition simply states that a constant flux of morphogens is injected at the origin . The simplest case is given by a constant degradation rate , which yields the exponentially decaying stationary profiles:

(3) |

Despite its simplicity, the exponential dependence of Eq. (3) captures surprisingly well the rapid concentration decay displayed by real profiles. However, the separate determination of and poses significant experimental difficulties and in most cases only the characteristic length can be unambiguously determined. This opens the door to the possibility of considering non-Markovian generalizations of Eq. (1) which are also compatible with the experimental results for the characteristic length.

## Iii Fractional reaction-subdiffusion equation with uniform reactivity

### iii.1 Derivation from a mesoscopic CTRW model with reaction

The starting point to derive the reaction-subdiffusion equation is the fundamental equation of a particle performing a continuous time random walk in the presence of a first order evanescence reaction (degradation). Consider the stochastic motion of an evanescent particle performing nearest-neighbor jumps on an infinite one-dimensional lattice. We shall hereafter denote the distribution of its waiting time between consecutive jumps by .The probability that a particle starting at site arrives at site in the time interval between and can be written as , where the arrival density obeys the following integral difference equation:

(4) |

In this equation, is the probability per unit time that a particle found at a given site at time has performed a jump up to time in the presence of a uniform evanescence reaction (since evanescence and jump are independent processes, the probability of jump is simply multiplied by the survival probability ). Taking the diffusive limit of Eq. (4) one can show via suitable Fourier-Laplace techniques Langlands that the probability to find the particle at position after a time given that it was initially at obeys the following equation:

(5) | |||||

The operator is defined via the equation

(6) |

where is the Laplace transform of the function and denotes the inverse Laplace transform. The operator is closely related to the Riemann-Liouville fractional derivative

(7) |

In fact, and are the same when applied to sufficiently regular functions , as determined by the condition . This condition is actually fulfilled by all functions of relevant to the morphogen problem, hence we shall use in place of in what follows.

If one is dealing with more than one particle, the concentration follows the same kinetics as above, i.e. SSS ; Langlands

(8) |

As one can see, this equation is a non-trivial extension of Eq. (1) for the case of anomalous subdiffusion with constant reactivity . In the normal diffusion limit the Riemann-Liouville operator reduces to unit and one recovers Eq. (1) with a constant . On the other hand, in the absence of reaction () Eq. (8) reduces to the standard fractional diffusion equation, which yields sublinear growth of .

Turning now to the morphogen problem, Eq. (8) is to be solved subject to the boundary condition (2). The solution for the case of a particle source can be obtained from the propagator solution (corresponding the initial condition ) via the relation between the Laplace transforms. The solution in Laplace space is found to be

(9) |

The stationary solution is obtained from the final value theorem for the Laplace transform:

(10) |

### iii.2 Robustness of stationary profiles

Using Eq. (10) it is possible to study the robustness of the concentration profiles with respect to a perturbation in the incoming flux . To this end, we take a reference value of the concentration and assess how large the shift of the associated position

(11) |

becomes when is perturbed; the larger the shift, the smaller the robustness of the profile. The latter can thus be characterized by the inverse of the relative change of with respect to a characteristic length of the problem (e.g. the linear size of a cell), i.e.

(12) |

Inserting Eq. (11) into this definition we find

(13) |

## Iv Fractional reaction-subdiffusion equation with non-uniform reactivity

Seki et al. SekiJPCM07 have shown that a CTRW process described by a generalization of Eq. (4), namely

(14) | |||||

with yields the following reaction-subdiffusion equation:

(15) | |||||

In order to tackle the corresponding morphogen problem, it is first necessary to find the propagator solution of Eq. (15). To this end it is convenient to introduce a new function defined via the transformation

(16) |

in Laplace space. This function is readily found to fulfil the equation

(17) |

In what follows, Eq. (17) will be used to investigate the effect of a non-uniform reactivity for several special cases.

### iv.1 Piecewise constant reactivity

Here, we assume that the reactivity is given by a superposition of Heaviside functions, i.e. . In region () one has , whereas in region () one has . Let us respectively denote by and the solutions of Eq. (17) in the regions and . These functions must fulfil the continuity conditions

(18) |

and

(19) |

In contrast, an integration of Eq. (17) across the origin shows that the solution must be discontinuous there:

(20) |

Using Eqs. (18)-(20) one can find explicit expressions for the Laplace transforms , and . For one gets the stationary biexponential solution

(21) |

with

(22) |

The characteristic constants are

For we shall distinguish two subcases with different physical behaviour. For one asymptotically gets the exponential decay law

(23) |

In contrast, for the behaviour is different. For normal diffusion the profile becomes constant for large , i.e.,

(24) |

However, when the diffusion is anomalous one has

(25) |

### iv.2 Exponentially decaying reactivity

Here, we assume the decay law . While in this case Eq. (17) does not seem exactly solvable for finite , it is possible to find an exact expression of the steady state profile by techniques similar to the ones used above. The final result is

(26) |

where the ’s are modified Bessel functions and . As in the case of piecewise reactivity with , this expression displays a different behaviour for normal and anomalous diffusion. In the normal diffusion case () one gets a monotonically decreasing profile from the concentration value

(27) |

at the origin to the limiting value

(28) |

### iv.3 Power law reactivity

Next, we take with . Since a steady state was already attained for an exponentially decaying reactivity, this will also be the case under the present situation, which describes enhanced particle evanescence.The general solution of Eq. (17) for is given by the modified Bessel functions and with . In order to single out the Bessel function corresponding to the physical solution, we use the fact that the incoming flux must be equal to the amount of particles per unit time that disappear due to degradation, i.e.,

(29) |

This condition leads to different solutions depending on the value of . For one gets

(30) |

where . For large , the above stationary solution can be shown to go to zero as

(31) |

As , this expression tends to a constant limiting value in the normal diffusion case and grows as for .

When the solution is not given by a power law rather than by Bessel functions. One has

(32) |

with . For this solution goes to infinity, a constant value or zero depending on whether is positive, zero or negative.

## V Conclusions and Outlook

In the present work we investigate both analytically and numerically the behaviour of the stationary concentration profiles arising from fractional reaction-subdiffusion equations derived from a CTRW model with a superimposed death process. These fractional equations are a natural extension of classical reaction-diffusion equations traditionally used to study the problem of morphogen gradient formation. We consider the case of linear degradation with both a uniform and non-uniform reactivity. The formulation of the problem in terms of fractional diffusion equations turns out to be a key ingredient in the analysis of the properties of morphogen gradients. This approach allows us to exploit a plethora of powerful analytical techniques available from fractional calculus to tackle the morphogen problem.

In the uniform case one obtains exponentially decaying stationary concentration profiles. Their robustness with respect to changes in the incoming flux increases with increasing . Likewise, one can study the robustness of the profiles with respect to a perturbation in by introducing the quantity

(33) |

This issue will be the subject of future research.

In the non-uniform case the behaviour of the stationary profiles turns out to be very sensitive to the specific spatial dependence prescribed for and to the value of the anomalous diffusion coefficient . Moreover, for the case of a piecewise constant reactivity with and anomalous diffusion, we see that a discontinuous profile arises and no steady state is reached in region . This is a novel effect not seen for normal diffusion. For exponentially decaying reactivity the concentration goes to a constant limiting value far away from the source () when , but it grows without bound for . Finally, when the connectivity decays as a power law, the stationary concentration may go to zero, to a constant or to infinity depending on the values of the characteristic decay exponent and .

In view of the strong inhomogeneities encountered by the diffusing morphogens in the embryonic environment, we believe that the sensitivity of the concentration gradients to the form of may be relevant for the modeling of morphogen gradient formation and interpretation.

Up to the case of piecewise constant reactivity with , in the present work we limit ourselves to study the behaviour of the stationary concentration profiles. However, analytic solutions for transient profiles are available for some of the cases studied, and others can be investigated via numerical techniques for the inversion of Laplace transforms. Besides, a reaction-subdiffusion equation was recently obtained that generalizes the above results to the general case Yadav06 ; Yadav08 . As one could have guessed in view of Eqs. (8) and (15), this equation reads as

(34) |

Our aim is to use the above equation to investigate further problems related to morphogen gradient formation in future. Beyond this field of research, Eq. (34) can be applied to many other problems of interest characterized by different kinds of boundary conditions.

###### Acknowledgements.

This work was partially supported by the Ministerio de Ciencia y Tecnología (Spain) through Grant No. FIS2007-60977, by the Junta de Extremadura (Spain) through Grant No. GRU09038, and by the National Science Foundation under grant No. PHY-0855471.## References

- (1) R. Metzler and J. Klafter, Phys. Rep. 339 1 (2000).
- (2) S. B. Yuste and K. Lindenberg, Phys. Rev. E 76, 051114 (2007).
- (3) R. Borrego, S. B. Yuste, and E. Abad, Phys. Rev. E 80, 061121 (2009).
- (4) M. Ibáñes and J. C. Izpisúa, Molecular Systems Biology 4, 176 (2008).
- (5) O. Wartlick, A. Kicheva and M. González-Gaitán, Cold Spring Harb. Perspect. Biol. 1, a001255 (2009).
- (6) B. I. Henry, T. A. M. Langlands, and S. L. Wearne, Phys. Rev. E 74, 031116 (2006).
- (7) E. Abad, S. B. Yuste, and K. Lindenberg, Phys. Rev. E 81, 031115 (2010).
- (8) I. M. Sokolov, M. G. W. Schmidt, and F. Sagués, Phys. Rev. E 73, 031102 (2006).
- (9) G. Hornung, B. Berkowitz B, and N. Barkai, Phys. Rev. E 72, 041916 (2005).
- (10) K. Seki, A. I. Shushin, M. Wojcik, and M. Tachiya, J. Phys.: Condens. Matter 19, 065117 (2007).
- (11) A. Yadav and W. Horsthemke, Phys. Rev. E 74, 066118 (2006).
- (12) A. Yadav, S. M. Milu, and W. Horsthemke, Phys. Rev. E 78, 026116 (2008).