Reaction-diffusion with stochastic decay rates

Reaction-diffusion with stochastic decay rates


Understanding anomalous transport and reaction kinetics due to microscopic physical and chemical disorder is a long-standing goal in many fields including geophysics, biology, and engineering. We consider reaction-diffusion characterized by fluctuations in both transport times and decay rates. We introduce and analyze a model framework that explicitly connects microscopic fluctuations with the mescoscopic description. For broad distributions of transport and reaction time scales we compute the particle density and derive the equations governing its evolution, finding power-law decay of the survival probability, and spatially varying decay that leads to subdiffusion and an asymptotically stationary surviving-particle density. These anomalies are clearly attributable to non-Markovian effects that couple transport and chemical properties in both reaction and diffusion terms.

I Introduction

Due to the interaction of diffusion and reaction mechanisms, reaction-diffusion systems in fluctuating environments may develop collective behaviors that are very different from those occurring under well mixed conditions. Smoluchowski’s theory von Smoluchowski (1917) quantifies the interaction of diffusion and reaction for fast bimolecular reactions through an effective rate that is proportional to the molecular diffusion coefficient. This approach is valid under well-mixed conditions. Spatial and temporal fluctuations may lead to the segregation of the reactants Ovchinnikov and Zeldovich (1978) characterized by non-Poissonian encounter processes and broad first-passage time distributions Bénichou et al. (2010); Reuveni et al. (2010a, b), such that reaction kinetics on small and large scales may be very different ben Avraham and Havlin (2005). The sound understanding and quantification of the mechanisms by which heterogeneity on small scales leads to “non-classical” or “anomalous” kinetics on large scales plays a central role in applications as diverse as contaminant degradation and chemical transformations in geological media Steefel et al. (2005); Dentz et al. (2011) and chemical kinetics in crowded intracellular environments Schnell and Turner (2004). A number of approaches have been proposed to model reaction behaviors in heterogeneous environments, including fractional kinetic orders and time-dependent rate coefficients Kopelman (1986); Savageau (1995) as well as delayed-reaction equations Schnell and Turner (2004); Bratsun et al. (2005); Brett and Galla (2013); Tian (2013). Oftentimes, such effective approaches to explain non-exponential survival probabilities are phenomenology-based lumped parameter models Aris and Astarita (1989). Indeed, the variety of mechanisms leading to anomalous diffusion and kinetics precludes general answers to fundamental questions. For instance, are emergent anomalous kinetics better described by non-linear, or by non-Markovian evolution equations? We address this question by solving the random decay model for fluctuations characterized by broad distributions of transport and reaction time scales, obtaining reaction-subdiffusion equations. Insisting on an exact derivation is especially important in this case, since the prima facie reasonable approach of adding reaction terms to known subdiffusion equations Fedotov and Méndez (2002); Méndez, V. et al. (2006); Henry and Wearne (2000, 2002); Henry et al. (2005); Langlands et al. (2007) has been shown to be inconsistent with microscopic dynamics and kinetics Sokolov et al. (2006). We find anomalous kinetics associated with population splitting and identify the cause in non-Markovian, rather than non-linear effects. Furthermore, the transport is highly anomalous, with the particle density approaching a stationary state.

We do not make assumptions regarding the origin of the distribution of transport times, but rather take these properties as given. However, it is important to note there do exist derivations in the literature of transport properties, such as first-passage times, from characteristics of complex media. For instance, first-passage observables have been computed for diffusion on fractals or media with heavy-tailed trap distributions Bénichou et al. (2010), and by applying the mapping between random walks and vibrations to complex elastic networks Reuveni et al. (2010a, b).

The paper is organized as follows. In Sec. II we introduce the random decay model as the continuous time random walk (CTRW) in which the walker is subject to a random decay process during each waiting period. In Sec. III, we derive the generalized master equation for the particle density and derive the solution in Fourier-Laplace space. In Sec. IV, we connect the random rate model to a general stochastic process. We preview the main results in this connection: the mean square displacement approaches a constant. The reaction kinetics follow a power law governed by an evolution equation with a heavy-tailed memory kernel that couples microscopic transport and reaction parameters. In Sec. V.1, we show that ordinary diffusion subject to random decay leads to perfect mixing in the scaling limit, and is equivalent to decay at a single, average rate. In Sec. V.2, we consider random rates combined with heavy-tailed waiting times and derive a generalized fractional Fokker-Planck reaction-diffusion equation with reaction and diffusion kernels that couple reaction and transport. In Sec. V.3, we assume a power-law rate PDF and derive exact asymptotic expressions for the reaction-diffusion equations and solutions. In Sec. VI, we demonstrate localization by deriving the exact expression for the asymptotic, steady-state particle density as a two-sided exponential distribution.

Ii The random decay model

ii.1 Stochastic decay rates

Figure 1: Overview of the model. The particle makes random jumps until it decays. Darker green corresponds to faster jumping. Darker red corresponds to faster decay. For clarity, the decay rate (red) is shown only at sites that the walker occupies. (A) Dark green and dark red: the particle experiences a high decay rate for a short time. (B) Light green and light red: The particle is immobile for a long time with a small decay rate, and so may survive for a long time. Events like (B) cause anomalous kinetics. (C) Light green and dark red: Particle is immobile for a long time with high decay rate, and so has a high probability of decaying in this step.

In this section, we formulate our model of diffusion in a fluctuating physical and chemical environment as a continuous time random walk (CTRW) subject to random decay. That is, a particle of species performs a random walk and at the same time, undergoes an irreversible reaction . The particle waits a random time before making each step, and undergoes decay at a rate during this waiting period, where and vary randomly. Regions where long waiting times and slow decay coincide are responsible for subdiffusion and anomalous kinetics, as illustrated in Fig. 1. This type of quenched disorder in the reaction and diffusion properties may occur in heterogeneous geological media characterized by a spatial distribution of minerals and thus specific reactive surface, and porosity Steefel et al. (2005), which leads to scale effects in the reaction properties Meile and Tuncay (2006); Li et al. (2008).

In this paper, we aim at quantifying the impact of variability in the physical and chemical system properties on the reaction behavior. To this end, we make the simplifying assumption that each random waiting time is independent of all past waiting times, and each decay rate is independent of all past decay rates. In other words, we assume fully annealed disorder, ignoring possible effects of correlations between steps. Models of quenched disorder assume that the medium fluctuates slowly enough that the walker samples a static configuration. However, annealed disorder is inherent in many systems in which the timescales of thermal fluctuations and stochastically varying transport and kinetic parameters are not well separated. For example, the waiting times may arise from long excursions into cul-de-sacs where the duration of the excursion is due to the motion of the particle itself, and so is memoryless. This feature has been exploited in studying reaction-subdiffusion front propagation in spiny dendrites Iomin and Méndez (2013). Furthermore, the stochasticity may stem from a fluctuating internal state of a particle that interacts with the medium. Examples of the latter are conformational fluctuations of proteins Manzo et al. (2015) and enzymes Boehr et al. (2009). It is worth noting that, in distinction from quenched models, annealed models are often solvable. In some cases, these solutions provide insight into features of the corresponding quenched model that do not depend strongly on correlations. This is the case for the model considered here. We derive exact results showing anomalous kinetics due to memory effects. Preliminary numerical results show that exponents characterizing power-law behavior of the mean square displacement  (MSD) and survival probability are the same in annealed and quenched versions of the model Lapeyre and Dentz (), which suggests that the same memory effects are at play in the quenched setting.

The structure of the annealed model allows the following practical simplification, and at the same time, generalization. During the th step, the particle is subject to two random processes, one triggering spatial translation and the other the conversion from species to species . For example, the waiting time before translation may be due to thermally driven escape from a trap with random energy. And the reaction rate may depend on the random local concentration of a catalyst. Escape from the trap will interrupt the reaction. And reaction will effectively interrupt escape by removing the particle from the population whose concentration we are measuring. Since one can view the jumping as interrupting an ongoing reaction and starting a new one, the th step may be simulated as follows. Sample a waiting time from the PDF of and a rate from the PDF of . Then, sample a random decay time from the PDF


If , the reaction is interrupted and the particle jumps. On the other hand, if , the particle indeed decays. But, sampling first a decay rate, and then a decay time is mathematically equivalent to sampling the decay time directly from the PDF of a single-step decay time given by


where is the PDF of . It is worth noting that this procedure is equivalent to sampling a decay time from a PDF that is independent of the step number . In other words, Poissonian decay with a random parameter (the rate) is equivalent to non-Poissonian decay described a single PDF. In fact the argument above also works if 1) is due to averaging an arbitrary PDF over a random parameter, and 2) and are not independent. In the following, we derive the basic results in this generality for two reasons. Firstly, the results may be applied directly to more complicated situations. For example, we will suggest below a stochastic Michaelis-Menten scheme defined by coupling and . Secondly, as we discuss below, results from the study of a broad class of stochastic processes can be applied directly to the general formulation of our model Pal and Reuveni (2017).

ii.2 General formulation

We consider a molecular entity of type (the “particle”) performing a continuous-time random walk while undergoing random conversion to an entity of type . The particle’s position and time are given by


where , , and the initial density, that is the density of is denoted . Here, is the random waiting time before the th step, and is the step displacement. We call the time interval the th renewal period. The particle position at time is then given by


where the renewal process is the number of steps (renewal periods) taken by time . With each renewal period, we associate a single-step random decay time . The decay time may be due to a disordered Poissonian process as in (2) or some more complicated chemical or biological process. It represents the time it would take the particle to decay if it were not subject to transport. Thus, if the particle is still alive at time and , then the particle takes the th step before it decays, and thus survives decay. On the other hand, if , then the particle decays before it has a chance to take a step. The random renewal period during which the particle decays is given by . Thus, the random time at which the particle finally decays after zero or more periods is given by


Note that the clock tracks only the step waiting times, but not the decay time. The analysis is facilitated by this choice, that is, considering an ordinary CTRW for which we mark the special time . Note also, that this framework is different from kinetic Monte-Carlo approaches such as the (spatial) Gillespie method Gillespie (1977); van Kampen (2007), which treats both diffusive and reactive particle events on the same ground. In the present work, particles perform a spatial random walk according to (3), and they may react during the (physical) waiting time with a certain probability as detailed above. This approach is equivalent to the reaction-diffusion equation for the species concentration Sokolov et al. (2006) and to more general non-local reaction and reaction-diffusion equations as developed in the remainder of the paper.

We assume annealed disorder, so that each renewal period is independent. That is, , are independent and identically distributed (iid) copies of . For the moment, we allow that and may be dependent as they may be coupled by chemical and physical properties of the medium. Although we mention such situations in Sec. IX, for the bulk of the paper, we will assume that they are uncoupled. For simplicity we assume that the step displacement is independent of both and and is distributed according to the probability density function (PDF) satisfying


where angle brackets denote averages and is a microscopic length scale characterizing the typical jump length. In this paper we distinguish Laplace transformed quantities by a hat, and Fourier transformed quantities by a tilde. The Laplace-conjugate of is and the Fourier-conjugate of is .

The use of the word coupled above refers only to whether random variables representing microscopic quantities are independent. Below, we shall be concerned with whether the mesoscopic transport and reaction are coupled. There is no direct relation between these two concepts. In fact, we find mescoscopic coupling in the case that all microscopic random variables are independent.

It is worth noting that the random decay time takes the form of the generic first passage time (FPT) under reset Reuveni (2016); Pal and Reuveni (2017). Here the “passage” is completion of a reaction (decay) at time , and the reset time begins the reaction anew. Thus, the random decay model may be viewed as FPT under reset coupled with CTRW by identifying the reset time of the FPT with the waiting time of the CTRW. The identification of the survival time with FPT under reset allows one to immediately apply general results for FPT under reset Pal and Reuveni (2017), including expressions for the mean FPT and fluctuations.

The main quantities of interest are the following. We refer to the probability that the particle has not decayed up to time as the survival probability, denoted by


The evolution of gives information on the chemical kinetics and reaction dynamics averaged over the entire system. For example, under a constant decay rate , it decays exponentially as . The density of surviving particles is given by the particle density of the CTRW conditioned on survival, that is, on . We denote the density of surviving particles by


Recall that we use the word “decay” as a shorthand for any irreversible reaction . Since we assume that the entities are non-interacting, one may interpret as either the probability density for a single entity, or as the local concentration or number density of normalized to one. In the latter case, is the profile observed at time by an imaging technique that detects species , but not species . According to Bayes’ rule and (7), the joint particle density and probability of survival is then given by


For simplicity, we shall refer to as a “density”. Because is normalized to one, we have the marginal


Although and are the physically relevant quantities, is more accessible mathematically. Thus, we will calculate and obtain via (9) by dividing by . The mean square displacement , given by


measures the spatial extent of the surviving particles.

Finally, it is important to note that we focus on rate PDFs with a finite probability that the rate is either zero, or arbitrarily close to zero. This is because the anomalous behaviour is driven by the coincidence of very long waiting times with very long decay times. In Sec. VIII, we discuss how the anomalies are modified in the case that the minimum possible rate is greater than zero.

Iii Generalized master equation and solution

We derive the generalized master equation Klafter and Sokolov (2011) and solution for the random decay model. The particle density defined in (9) satisfies (See Sec. X.1.)

Here is the incoming live-particle flux at position at time . That is, is the probability that the particle is alive and makes a step between times and into the region between and And the factor
is the probability that the particle has survived both translation (i.e. has not escaped from a trap) and decay up to time during a renewal period beginning at time . Thus, the integral counts all particles that arrived at at some time in the past and that have neither jumped away, nor decayed in the time since arriving. The flux satisfies the Chapman-Kolmogorov type integral equation



is the joint probability and probability density to survive both translation and single-step decay until time and then to make a translation (jump) at time . Here is the PDF of the waiting time conditioned on waiting time smaller than decay time, . Eq. (13) expresses particle balance under reaction losses. Note that Eqs. (12a) and (12c) have the same structure as the governing equations of the classical CTRW. But in fact they do not describe a CTRW and the model cannot be cast as one. This is due to the presence of reactions, with the result that is not a waiting time PDF and is not a translation survival probability. Instead the decay of the single-step survival probability includes two loss terms representing translations and decay. Taking the derivative of (12b) we find




is the joint probability to both survive single-step decay and not make a jump until time , and then to decay at time .

The system (12) can be combined into the generalized reaction-diffusion Master equation (GME) (See Sec. X.2.)


where we define the reaction kernel and the diffusion kernel through their Laplace transforms as


The diffusion kernel quantifies the impact of random decay and waiting times on the spatial motion, the reaction kernel on the particle survival. As is evident in the kernels, the reaction and transport processes are intimately coupled.

Integration of (16) over space reveals the dynamics that govern the reaction kinetics,


Eq. (18) is of central importance. It expresses the impact of segregation of the reactants and the different reaction and transport histories on the overall reaction behavior. The Markov property of the reaction process that underlies the exponential model breaks down in the presence of distributed reaction and diffusion rates. It is interesting to note that the Gillespie method Gillespie (1977) is a fully Markovian method, in which interreaction waiting times are exponentially distributed. The non-local nature of 18 indicates that this is no longer valid in spatially heterogeneous systems. This is the subject of ongoing work.

The solution to the GME in Fourier-Laplace space is the generalized Montroll-Weiss equation for the particle density (See Sec. X.3.1.)


Note that (19) involves defined in (13) and defined in (12b), but not defined in (15). This is because both and characterize the event that decay does not occur, while characterizes the event that decay does occur. This is to be expected, since is the density of particles that have not decayed.

Setting in Fourier space is equivalent to integrating over in real space. Thus, putting in (19) and referring to (10), we obtain the expression for the survival probability

The mean survival time is given by , from which we obtain the simple form


As mentioned above, the survival time is formally a FPT under reset. A simple, alternative derivation of (20) from this viewpoint is found in Ref. Pal and Reuveni (2017). From (20) we see that increases with 1) increasing probability of large values of both and , and 2) increasing probability of .

Iv Stochastic rates and anomalous kinetics

Figure 2: (Black solid) Survival probability and (green solid) mean square displacement defined in (11) for heavy-tailed waiting time PDF, and power-law reaction rate PDF with , , and . Left and right ordinate axes differ in physical dimensions, but are numerically equal. (Dotted) Exponential short time behavior of the survival probability, which is characterized by the average rate . (Lower dashed) Asymptotic power-law decay . (Dash-dotted) Short-time power-law behavior . (Upper dashed) Asymptotically constant occurring on the localization time scale .

We now assume that the single-step decay time and the translation waiting time are uncoupled, that is, and are independent. We denote the PDF of waiting times by


where is dimensionless and is the waiting time scale. Furthermore, we adopt the viewpoint of Sec. II.1 that the randomness in decay is due to first-order decay with rates that vary stochastically, but are constant during each renewal period. The variability in rates is characterized by the PDF


where is the time scale of single-step decay, and is a dimensionless PDF.

A main result and key message of this paper is that strong physical disorder expressed via (21) combined with disordered rates expressed via (22) leads to anomalous kinetics as well as anomalies in transport beyond standard subdiffusion. We quantify these anomalies and identify their source in long reaction memory rather than nonlinearity. The anomalous kinetics and transport are clearly evident in Fig. 2, which shows the survival probability and the mean square displacement for a heavy-tailed waiting time PDF that behaves as with , for larger than the characteristic time in (21), and a rate PDF that behaves as with for smaller than the characteristic rate . We observe two remarkable behaviors. Firstly, the survival probability decays as a power-law , where


and secondly, increases proportionally to , as for non-reacting particles, until a characteristic reaction time scale after which it decays towards a constant. These two behaviors indicate a localization of the density of surviving particles.

The power-law decay of the survival probability observed in Fig. 2 can be modeled by a non-linear kinetic rate law as Kopelman (1986)


with as an effective reaction rate. While this equation gives the power-law decay , it implies a conceptual framework that is clearly inconsistent with the correct evolution equation (18). Indeed (18) is linear but non-Markovian, implying history-dependent evolution. This is an important point as the conceptual framework influences the approach taken to more complicated scenarios. A discussion of its importance in interpreting experiments is found in Ref. Sokolov (2012).

V Reaction-diffusion dynamics

Spatial fluctuations in biological and physical systems provide our main motivation for assuming that the random decay has its origin in disordered rates. Thus, in the following analysis we assume that the single-step decay time arises from averaging decay over random rates. However, it may be useful to go in the opposite direction. Given a distribution for , compute the corresponding distribution rates. In this way our results, although explicitly written in terms of random rates, may be applied to random decay times of varying physical origin. The PDF is obtained from that of as follows. Referring to (1), it is easy to see that the single-step decay-survival probability is given by


Since (25) is the Laplace transform of , it may be inverted for any density of for which the inverse Laplace transform exists.1

We begin by writing the GME in terms of random rates. Using the independence of and and referring to (25), we find , , and , where the translation survival probability is given by


is the probability, in the absence of decay, that the particle has not taken a step during a renewal period before time . Thus, the solution (19) to the GME (16) may be written as


Eq. (27) is the basis of the following analysis. We will describe the conditions under which on the one hand, the system becomes well-mixed and exhibits homogeneous kinetics at long times, and on the other hand the system remains poorly-mixed and exhibits persistent physical and chemical anomalies. In the following, we assume that the rate density (22) has weight at, or in the neighborhood of, , leaving the more general case to Sec.VIII.

The main factors determining the evolution of and , and the degree of mixing in particular are 1) whether the mean waiting time between jumps exists, ie . 2) The relative magnitude of the three time scales: the time scale of microscopic transport defined in (21) , the time scale of reactions from (22), and the physical time . If , and


then the system tends to a well-mixed state with homogeneous kinetics as the time scales separate. This is because at long times surviving particles have typically made many steps and will sample many rates before dying. On the other hand, if diverges, then the system never becomes well mixed, no matter how large the scale separation in (28). Instead, we find anomalous kinetics and dynamics due to memory effects. This corresponds to particles that are trapped for long times in regions of low reactivity.

v.1 Well-mixed scenario.  

We first treat the case . We consider the scaling limit in order obtain a mesoscopic picture in which observational length and time scales are much larger than the microscopic scales. In the scaling limit and such that


converges to a positive constant, the GME (16) reduces to


provided . (See Sec. X.4.) Eq. (30) gives a mesoscopic description of evolution of the system. The local density changes little over a short time, but this time represents an infinite number of steps. Of course, physically, displacements and waiting times may be very small, but must be finite. Thus, eq. (30) is an accurate description insofar as the microscopic and observational physical scales are well separated. For a colloidal system, the waiting time is the time between collisions of solvent molecules with a relatively massive particle, so that the ratio of the time required for the particle to move a distance equal to its own radius and the time between collisions may be orders of magnitude or more. On the other hand, if the waiting times are dominated by trapping, then the timescale of the trapping plays the role of the microscopic time scale and the separation between and the mesoscopic scale may not be as large. We present an example of the latter case for in Sec. X.6.

In the present case, (30) describes ordinary diffusion with a constant, homogeneous decay rate. The survival probability satisfies the first-order rate equation


The effective reaction kinetics are determined solely by the characteristic reaction rate. Eq. (30) makes evident that on the mesoscopic level the kinetics are effectively homogeneous in space. Note that this behavior is also observed in general at times shorter than both the reaction and translation time scales and . In this case, the reaction kernel also reduces to . This is obtained from (17) by considering the limit and .

The scaling limit leading to (29) and (30) involves letting approach zero. Since we do not rescale the reactions, this implies , which corresponds to a small Damköhler number. Eq. (31) immediately gives us the mean survival time . The extreme opposite to the scaling limit is and corresponds to large Damköhler number. In this case, the mean survival time is just the mean single-step decay time . This can be seen by noting that for it is highly probable that . Thus, the numerator in (20) is approximately and the denominator approximately . Furthermore, we note that , and use (25) to arrive at


provided the moment exists. There is no mixing at all, and is dominated by particles that never jump, but instead decay in their initial environments. The intermediate behavior between these extremes depends strongly on details of the distributions of both the rate and of the waiting-times, rather than just their asymptotics. We defer a discussion of these more complicated and varied results to Sec. X.5.

Figure 3: Survival probability for a heavy-tailed as in (33) and a reaction-time PDF as in (40) with . (a) , and (from uppermost to lowermost curve) . (Dashed) Asymptotic form (44). Filled circles indicate the localization time . (b) , with (from uppermost to lowermost curve) . (dashed) . For both (a) and (b), solid line curves are numerical inversion of (38), symbols are Monte Carlo simulations of the microscopic model.

v.2 Fractional reaction-diffusion.  

As discussed in the Introduction and indicated in Figs. 1 and 2, broad waiting time distributions lead to anomalously long particle survivals if they coincide with small or vanishing reaction rates. This inhibits mixing and leads to segregation. We mentioned above that only in the case does this segregation persist in the scaling limit of vanishing waiting time scale . We now turn our attention to this scenario. We find that the mesoscopic reaction-diffusion equation possesses kernels that couple the independent microscopic physical and chemical fluctuations, which manifest the non-Markovian reaction kinetics. To illustrate, we consider the heavy-tailed waiting time PDF that behaves as


for times larger than the microscopic time scale . Eq. (33) implies that diverges. Physically, this corresponds to waiting (or trapping) times that occur on all time scales, including the duration of an experiment. The variation in trapping time may be due to thermal activation over a random binding energy, or to long, slow, excursions in inclusions, or many other causes Bouchaud and Georges (1990); Metzler et al. (2014).

The correct scaling limit to employ with (33) is and such that


converges to a positive constant. (See Sec. X.4.) In this limit, the evolution of the particle density is determined by the non-Markovian reaction-diffusion equation (See Sec. X.4.1.)


where the reaction and diffusion kernels are defined by their Laplace transforms


For , (36) and (35) reduce to the well-known fractional Fokker-Planck equation. It is worth noting that the operators in (35) describing subdiffusion with random decay rates may related to fractional calculus via rate-averaged tempered fractional calculus Sabzikar et al. (2015).

The solution to (35) is (See Sec. X.4.)


We have assumed here that for simplicity. The corresponding survival probability obtained by setting assumes the compact form


Setting in (38), we see that the mean survival time of the particle under random diffusion and decay given by (20) takes the form


It is important to note that the scaling limit does not exist if the PDF of the rates decays more slowly than as . In this case the denominator of (39) diverges, so that the mean survival time . Likewise, the kernel in (36) diverges for all , and in (38) is identically zero. Physically, means that the rates are sampled very rapidly and for a heavy-tailed rate PDF there is a high probability of very fast rates. On the other hand, if diverges more rapidly than as , then the numerator of (39) diverges, so that the mean survival time diverges. But, in this case, the scaling limit still exists. For instance, for finite time , (36), (37), and (38) are well defined.

v.3 Broadly distributed mean reaction times.  

In this section, we focus on rate PDFs that decay as a power-law for much smaller than the inverse of the characteristic time


where . This implies a power-law PDF of the mean reaction times for . Substituting (40) into (25) we see that the probability to survive decay in a single step varies asymptotically as (See Sec. X.4.2.)


In general the kernel approaches the inverse of the mean survival time (39) at a time comparable to the reaction time scale . However, for power law rates (40) and , with given by (23), computing (39) gives , and (36) gives . In fact, in this case, both kernels (36) take a particularly simple form, and . Thus, for , (35) becomes the fractional reaction-diffusion equation (See Sec. X.4.3.)


where and . Although the microscopic reactions are first-order, the macroscopic reaction term in (42) is non-Markovian with a memory kernel that couples the microscopic transport and kinetic parameters. This is made clear in the equation governing the evolution of the survival probability


which is obtained by integrating (42) over .

As mentioned earlier, at short times , the survival probability is approximately exponential, . In the case of power-law distributed rates (40) we obtain from (38) the explicit, long time solution (See Sec. X.4.2.)


Eq. (44) shows that, as anticipated in the definition of the fractional-order derivative of (42), the exponent observed in Fig. 2 is given by , which manifests again the intimate coupling of diffusion and reaction mechanisms in the mesoscopic limit. Fig. 3 shows the dependence of on and , and the excellent agreement of the derived analytical expressions with Monte-Carlo simulations of the microscopic model. In Sec X.6 we give a detailed description of the Monte-Carlo algorithms.

Figure 4: Density of surviving particles with and as in Fig. 2, and (black), (yellow), (blue). The curve for agrees with (45). The inset shows the same curves on a semi-log scale. Curves are numerical inversion of (78).
Figure 5: Mean square displacement defined in (11). Symbols are inverse Laplace transform. Lines are stochastic simulations. Heavy-tailed waiting time PDF (33) with , . Rate PDF (54) with and (Crosses) , (Squares) .

Vi Localization of particle density.  

By localization, we mean that the surviving-particle density approaches a stationary density . This density has exponential tails and a well-defined, constant, width . Thus, a measurement that detects the local concentration of the surviving species , but not that of the product species , will be characterized by this width . Since the surviving-particle density ignores the product species , it is independent of the fate of , which depends on the application. For instance may be removed from the system. Or, it may be invisible to the detector but is either immobilized or continues to diffuse. It is interesting to consider the case that species is immobile, but it is detected along with . In this case, the sum of the local concentrations of and approaches the same stationary density obtained by considering species alone Lapeyre and Dentz ().

Localization does not occur in the well-mixed case studied in Sec. V.1. On the contrary, the decay is spatially uniform. This is evident by first noting that satisfies (30) where is the particle density for ordinary diffusion with no decay, ie . Then referring to (9), we see that this implies , which means that the decay is independent of the transport. Finally, substituting this last equality into (11) shows that the MSD evolves exactly as in the non-reactive case, increasing without bound. However, in the case of strong chemical and physical fluctuations, when the system remains poorly-mixed, the particles are localized at long times. The surviving-particle density tends to a stationary state , given by (See Sec. X.4.)


where the localization time and localization length are given by


The corresponding MSD approaches a constant value given by

This localization is clearly verified and illustrated in both the MSD in  Fig. 2 and Fig. 5, and the particle density in Fig. 4. The MSD approaches a constant at long times. As , the density of surviving particles approaches (45) which is represented by the blue curve in Fig. 4. Note that the localization time marks the scale at which the mean square displacement crosses over from the power-law behavior to the constant value, as illustrated in Fig. 2. The deviation of from a power-law for in Fig. 5 is due to corrections to the scaling limit. See Sec. X.6 for details of the numerical methods.

To recap, we have derived the fractional reaction-diffusion equations (35), (42) and fractional kinetic equation (43) in the scaling limit of the random walk. These are exact solutions of the microscopic model with no homogenization or upscaling. The presence of memory kernels coupling the transport and kinetic parameters manifests the poor mixing, even in the scaling limit, in contrast to the perfect mixing in the scaling limit for Brownian diffusion (30). We have derived exact expressions in the scaling limit for the particle density (37) and survival probability (38). We presented the asymptotic solutions for the survival probability (44) and for the localized (stationary) particle density (45) and (46). These derivations and their physical interpretation are the main results of Sec. V and Sec. VI.

Vii Coupled vs. uncoupled reaction.  

To better understand stochastic decay, it is useful to compare the mesoscopic description of the random decay model to that of other models of reaction-subdiffusion. We refer to a model in which the reaction proceeds independently of the transport as “uncoupled”. Otherwise, it is “coupled”. The question of whether a model is coupled or uncoupled is an instructive point of comparison, which we address in the following.

vii.1 Uncoupled reaction

Suppose is the density of surviving particles undergoing subdiffusion and an unspecified decay process. Define via


where is a solution to the fractional Fokker-Planck equation with no decay Metzler and Klafter (2004)


and is the Riemann-Liouville fractional derivative Metzler and Klafter (2004). Substituting into (48) we see that satisfies the equation


By construction, (49) holds formally for any density , with given by (47). But it is evidently only meaningful if results from a particle that diffuses according to (48), and is subject to decay that is independent of the dynamics Sokolov et al. (2006); Sagués et al. (2008); Abad et al. (2010); Fedotov (2010); Abad, E. et al. (2013); Yuste et al. (2013). This becomes clear upon considering the time rate of change of mass at position and time


Using (47) we write (50) as


The first term on the right hand side is the rate due to transport. We are interested in the second term , which is the instantaneous decay rate (times ) at position . The role of the second term as a time and space dependent decay rate is also clear in the last term in (49). The diffusion and decay in (51) are manifestly independent. It is important to note that the reaction term in (49) is Markovian, that is, local in time. Indeed, integrating (49) over space, we find an equation for the survival probability


If we allow to depend on the density itself, then (52) is non-linear Fedotov (2010). Still, the decay is independent of the dynamics and is Markovian. Eq. (49) has been derived for many models of uncoupled dynamics and decay, appearing, for example, as Eq. (20) in Ref. Sokolov et al. (2006), Eq. (29) in Ref. Abad et al. (2010), and Eq. (23) in Ref. Fedotov (2010). Typically, is independent of , so that any -dependence in represents independent, externally imposed, spatially varying decay.

vii.2 Coupled reaction

The random decay model strongly couples chemical kinetics and transport, which results in a very different description and behavior. There is no explicit space-dependent decay in the microscopic model of Sec. II.2 as there is in (49). However the strong coupling results in a non-Markovian reaction term in (35) and in (42), which in turn gives rise to a time- and space-dependent decay rate . The decay rate is that part of the time rate of change of the density that is not due to transport. It is an effective or mesoscale rate that emerges from microscopic kinetics and transport that have no explicit space dependence.

To compute for the random rate model, we begin by dividing (35) by , thereby obtaining an expression for the time rate of change of the mass that is analogous to the expression for independent decay (51). Then is given by the last term in (35) divided by ,


depends on the history of the particle density at through the kernel . Thus, it is the non-Markovian operator that induces a spatial dependence in the effective decay rate.

Figure 6: Decay rate defined in (53). Times from uppermost to lowermost curve (black,gold,light blue,green,yellow,dark blue, orange): . Waiting time PDF (33), , Decay rate PDF (97), , . Scaling limit with generalized diffusivity .

The solution to (53) by numerical inversion of the Laplace transform is shown in Fig. 6. At short times , the decay is uniform and exponential with rate . This corresponds to the dotted line in Fig. 2. At intermediate times, Fig. 6 clearly shows a strongly inhomogeneous decay rate. Because decays as a power at long times, the instantaneous decay rate averaged over space decreases like . As increases, the decay rate near approaches zero, but the asymptotic value as approaches is . This suppression of the decay rate in the central part of the density is responsible for the localization discussed in Sec. VI.

Another case of coupling transport and decay is that in which the walker does not decay while waiting, but rather only before or after making a step Henry et al. (2006); Abad et al. (2013); Yuste et al. (2013). Suppose a fraction of walkers are removed at the beginning of each waiting period. Compare this to the random decay model with rate density


which means that during each waiting period the walker suffers no decay with probability and decays at rate with probability . It can be shown that (42) holds in this case with . For times , the longest trapping times are important, so that the particle decays very early in the waiting period. This is equivalent to removing the walker with probability at the beginning of the step. The fractional reaction-diffusion equation for the latter model given in Ref. Henry et al. (2006) is indeed equal to (42) with .

Viii Lower cut-off in rate PDF.  

Thus far, we have considered rate distributions with rates arbitrarily close to zero. But, suppose we shift , that is, let with so that the probability that is zero. We show in Appendix X.3.2 that the propagator for the shifted reaction rates is


where is the solution for the unshifted density . Note that this leaves unchanged, so that shows the same localization as . However, the asymptotic survival probability (44) now decays exponentially fast with the smallest rate . For instance, for the power law  (40) and the survival probability follows the truncated power law .

Ix Conclusions and Outlook

We derived the mesoscale behavior of a reaction-diffusion system characterized by microscopically fluctuating transport and reaction kinetics, using the framework of a continuous time random walk that samples disordered decay rates. We showed that broadly distributed waiting and reaction times give rise in the scaling limit to a generalized fractional reaction-diffusion equation with non-Markovian reaction and diffusion operators both of which are characterized by intimate coupling of microscopic chemical and physical parameters. This equation describes a system that asymptotically remains poorly mixed leading to power-law kinetics and spatially inhomogeneous reactions. The resulting decay is manifest in a particle density whose profile differs radically from that given by nonreactive subdiffusive CTRW, most notably in a stationary (localized) particle density at long times. This is in stark contrast to the case of ordinary diffusion in the scaling limit, which experiences spatially uniform decay with a rate equal to the average of the disordered rates.

Understanding of the mechanisms by which mesoscale behavior emerges from microscopic disorder plays a key role in diverse physical systems. For example, observed scale effects in reaction laws and decrease in reactivity on large scales in geological media Kump et al. (2000); Li et al. (2008); Meile and Tuncay (2006) can be attributed to spatial heterogeneity in the chemical and physical medium properties. Sometimes these behaviors are modeled by empirical non-linear reaction rate laws Peters et al. (2006). Our results show that physical and chemical system fluctuations are unambiguously attributable to non-Markovian, but linear kinetics. The segregation of reactions, here a mobile and an immobile species, in the presence of fluctuating chemical properties, leads to a broad distribution of effective reaction time scales, which are composed of both transport and reaction times. The reaction process itself is history dependent, as expressed by the non-local kinetic rate law (18). This new understanding of the role of chemical and physical fluctuations provides a systematic way towards quantifying effective large scale reaction behaviors and scale effects in reactivity in terms of the physical and chemical heterogeneity of the host medium in natural and engineered media. Furthermore, the results derived for first-order decay can be generalized to more complex chemical reactions under stochastic reaction and transport rates along the lines of the approach presented in Hansen and Berkowitz (2015).

We have focused on transport in the presence of random translation times and decay rates. However, it is important to point out that the theory presented here is independent of the specific physical context in which it was developed. We derived the main results for a general stochastic framework that combines two processes, CTRW and first passage under restart, by identifying the CTRW waiting time with the restart time. Applications and mathematical properties of CTRW Meerschaert and Straka (2014); Metzler et al. (2014) and first passage time (FPT) under restart Pal and Reuveni (2017); Reuveni (2016) have been studied intensively. But the fruitful union of these two theories remains nearly unexplored. Possible avenues can be found in the many diverse processes that determine statistics of the step displacement  Appuhamillage and Bokil (2011); Ramirez et al. (2013), waiting time  Bénichou et al. (2010); Reuveni et al. (2010a, b), and single-step decay time  Dentz et al. (2011).

However, an important class of chemical processes, namely Michaelis-Menten (MM) reactions Kou et al. (2005), require further generalization of FPT under reset. In recent years, advances in single-molecule spectroscopy have opened the possibility of measuring and controlling Lomholt et al. (2007) catalysis on the level of single, or a few, molecules. This in turn has spurred the development of stochastic approaches to MM reactions. These include considering the effects of fluctuations Grima (2009); Pulkkinen and Metzler (2015), internal states of the enzyme Kolomeisky (2011), and non-Poissonian processes. A stochastic Michaelis-Menten scheme is obtained from the generic FPT under reset by delaying restart of the process by a random time after each interruption. In catalytic reactions, represents the rebinding time. In this stochastic formulation, recent theoretical studies have predicted experimentally accessible Roeffaers et al. (2006), counter-intuitive kinetics by replacing the classical Poissonian processes governing binding, unbinding, and catalysis times with non-Poissonian processes Wu and Cao (2011); Reuveni et al. (2014); Rotbart et al. (2015). The importance of extending this Michaelis-Menten scheme to include heterogeneous catalysis due to a fluctuating environment has been recognized in recent experimental Janssen et al. (2014) and theoretical Reuveni et al. (2014) work. An attractive possibility is to modify the framework presented herein by including the rebinding time . This immediately yields a Michaelis-Menten scheme capable of handling heterogeneous catalysis via diffusion following unbinding events. The challenge of understanding the interplay of transport and Michaelis-Menten-like processes in cellular environments Schnell and Turner (2004); Reuveni et al. (2014) is a particularly promising candidate for such a Michaelis-Menten-CTRW approach, given that macromolecular crowding in cells may lead to both CTRW-like subdiffusion Metzler et al. (2014); Barkai et al. (2012) and modified binding dynamics Morelli et al. (2017); Guigas and Weiss (2008); Zhou et al. (2008).

X Appendix

x.1 Derivation of the integral equations for the propagator

Here we derive (12) from the microscopic model given in Sec. II. The assumption of an unbiased walk with finite step variance [See (6) and (72).], will be employed when passing to the scaling limit, but does not enter here. The particle density (9) may be written

Recall that the indicator function is if the argument is true and otherwise. The factor may be decomposed as follows. From the transition rules (3) and (5) it follows that, at a given time , a particle that is at position has survived until the turning time with probability , and has survived the last time interval from the last turning point to the present time with probability . Referring to (4), this implies that the particle density at time is given by

where the random variable is the number of steps performed up to time . We partition the probability space into disjoint sets, so that the expectation becomes a sum of expectations

The last factor combines the requirements that particle has neither decayed nor jumped during the increment . We now separate explicitly the contributions up to the last turning point at time and during the final resting interval

Because the last factor depends only on and , which are independent of the remaining factors, we split the expectation into two factors obtaining


where we have also used the fact that the step variables share a common distribution. We now define



Eq. (58) denotes the joint probability density for a particle to arrive at position at time on the th step. With these definitions, (56) is rewritten


We analyze by writing in the following form

That this is indeed the expression for can be seen by performing the integrals and eliminating either one of the delta functions for and either one for , and using (3). Note that the only random variables appearing in the last three factors in the expectation are , and , while the first three factors depend only on random variables for . The last three factors are thus independent of the first three and we can again factor the expectation. Furthermore, per the Dirac delta , we have <