Anomalous reaction-diffusion equations for linear reactions
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug Anomalous reaction-diffusion equations for linear reactions
Sean D. Lawley ∗ University of Utah, Department of Mathematics, Salt Lake City, UT 84112 USA (Dated: August 27, 2020)Deriving evolution equations accounting for both anomalous diffusion and reactions is notoriouslydifficult, even in the simplest cases. In contrast to normal diffusion, reaction kinetics cannot beincorporated into evolution equations modeling subdiffusion by merely adding reaction terms tothe equations describing spatial movement. A series of previous works derived fractional reaction-diffusion equations for the spatiotemporal evolution of particles undergoing subdiffusion in one spacedimension with linear reactions between a finite number of discrete states. In this paper, we firstgive a short and elementary proof of these previous results. We then show how this argumentgives the evolution equations for more general cases, including subdiffusion following any fractionalFokker-Planck equation in an arbitrary d -dimensional spatial domain with time-dependent reactionsbetween infinitely many discrete states. In contrast to previous works which employed a variety oftechnical mathematical methods, our analysis reveals that the evolution equations follow from (i)the probabilistic independence of the stochastic spatial and discrete processes describing a singleparticle and (ii) the linearity of the integro-differential operators describing spatial movement. Wealso apply our results to systems combining reactions with superdiffusion. I. INTRODUCTION
The signature of a normal diffusive process is that themean-squared displacement grows linearly in time. Thatis, if X ( t ) denotes the one-dimensional position of thediffusive particle at time t ≥
0, then E (cid:2)(cid:0) X ( t ) − X (0) (cid:1) (cid:3) ∝ t, (1)where E denotes expected value. However, the mean-squared displacement in complex systems often deviatesfrom the linear behavior in (1) and instead grows as apower law, E (cid:2)(cid:0) X ( t ) − X (0) (cid:1) (cid:3) ∝ t α , α > , (2)in a phenomenon called anomalous diffusion if α = 1.Subdiffusion is defined by (2) with α < α > ∂∂t c ( x, t ) = D − αt K α ∂ ∂x c ( x, t ) , x ∈ R , t > , (3) ∗ [email protected] for the subdiffusive chemical concentration c ( x, t ) at posi-tion x at time t . In (3), the parameter K α > (time) − α )and D − αt is the Riemann-Liouville fractional derivative[9], defined by D − αt f ( t ) = 1Γ( α ) dd t Z t f ( s )( t − s ) − α d s, (4)where Γ( α ) is the Gamma function. Note that D − αt issometimes denoted by ∂ − α ∂t − α . As a technical aside, theoperator appearing in the derivation of (3) is actuallythe Gr¨unwald-Letnikov derivative, but this operator isequivalent to (4) for sufficiently smooth functions [10].Generalizing (3), fractional Fokker-Planck equationsmodel the spatiotemporal evolution of subdiffusivemolecules under the influence of an external force [11].A fractional Fokker-Planck equation takes the form ∂∂t c ( x, t ) = D − αt L x c ( x, t ) , x ∈ V ⊆ R d , t > , (5)where V ⊆ R d is a d -dimensional spatial domain and L x is the forward Fokker-Planck operator, L x f ( x ) := − d X i =1 ∂∂x i (cid:2) µ i ( x ) f ( x ) (cid:3) + 12 d X i =1 d X k =1 ∂ ∂x i ∂x k h(cid:0) σ ( x ) σ ( x ) ⊤ (cid:1) i,k f ( x ) i , (6)where µ ( x ) : V R d is the external force (drift) vectorand σ ( x ) : V R d × m describes the space-dependenceand anisotropy in the diffusivity. Of course, if α = 1,then D − αt is the identity operator and (5) reduces tothe familiar equation of integer order, ∂∂t c ( x, t ) = L x c ( x, t ) . (7)A fundamental and now longstanding question is howto model reaction kinetics for subdiffusive molecules (seethe review [12] and [13–22]). In the classical case of nor-mal diffusion, reaction terms can simply be added to theevolution equations describing spatial movement. Moreprecisely, consider the vector of n chemical concentra-tions, c ( x, t ) = ( c i ( x, t )) ni =1 ∈ R n , (8)where c i ( x, t ) denotes the concentration of species i atposition x ∈ R d at time t ≥
0. In the absence of spatialmovement, suppose the concentrations obey the mean-field reaction equations ∂∂t c = f ( c ) , (9)where f : R n R n . In the case of normal diffusion whereeach chemical species moves by (7), one incorporates thereaction kinetics (9) into spatiotemporal evolution equa-tions by the simple addition of f ( c ) to the righthand side, ∂∂t c = L x c + f ( c ) . (10)However, this simple procedure fails for subdiffusion. In-deed, it was shown that the following attempt to combinesubdiffusion with degradation at rate λ > ∂∂t c = D − αt K α ∂ ∂x c − λc, x ∈ R , t > , leads to an unphysical negative concentration, c ( x, t ) < R with an irreversible reaction be-tween n = 2 chemical species. Equivalent equationswere then derived in [23] and [25] using different for-malisms. These results were generalized in [26] to al-low reversible reactions between any finite number n ofchemical species. In particular, in the case that (i) thereactions in (9) are linear, f ( c ) = R c , where R ∈ R n × n is a constant reaction rate matrix, and(ii) each chemical species moves by the one-dimensionalfractional diffusion equation in (3), it was found that [26] ∂∂t c = K α e Rt D − αt h e − Rt ∂ ∂x c i + R c , x ∈ R , (11)where e Rt is the matrix exponential. In contrast tothe simple form in (10) with decoupled reaction andmovement terms, notice that the reactions modify themovement term in (11). Interestingly, for the caseof L´evy flights with an irreversible reaction between n = 2 species, it was shown in [25] that the reaction-superdiffusion equations have the usual decoupling of re-action and movement terms. The results in [23–26] werederived using continuous-time random walks and Fourier-Laplace transform theory.In this paper, we first give a short and elementaryproof of (11). We then show how this argument givesthe evolution equations for more general cases, includ-ing subdiffusion following any fractional Fokker-Planckequation in an arbitrary d -dimensional spatial domainwith time-dependent reactions between infinitely manydiscrete states. This analysis reveals that the evolutionequations follow from (i) the probabilistic independenceof the stochastic spatial and discrete processes describ-ing a single particle and (ii) the linearity of the integro-differential operators describing spatial movement. Inaddition, under mild assumptions on initial and bound-ary conditions, the evolution equations imply that thespatial and discrete processes are independent. That is,under some mild conditions, the evolution equations holdif and only if the spatial and discrete processes are inde-pendent.The rest of the paper is organized as follows. In sec-tion II, we give a simple argument that yields (11). Insection III, we generalize this argument to yield the evo-lution equations describing more complicated spatial anddiscrete processes. In section IV, we apply this more gen-eral result to some examples. We conclude by discussingour results and highlighting future directions. II. SIMPLIFIED SETTING
We first consider a setup that is equivalent to themain problem considered in [26]. Assume { J ( t ) } t ≥ is acontinuous-time Markov jump process on the finite statespace { , . . . , n } . Suppose the matrix R ∈ R n × n containsthe transition rates, meaning the distribution of J ( t ) sat-isfies the linear ordinary differential equation,dd t r = R r , (12)where r ( t ) is the vector of probabilities, r ( t ) = ( r i ( t )) ni =1 := (cid:0) P ( J ( t ) = i ) (cid:1) ni =1 ∈ R n . (13)Of course, the solution to (12) is the matrix exponential, r ( t ) = e Rt r (0) , t ≥ . (14)In the language of Markov chain theory, R is the forwardoperator and the transpose R ⊤ is the backward operator(i.e. R ⊤ is the infinitesimal generator [27]).Assume that { X ( t ) } t ≥ is a one-dimensional subdiffu-sive process taking values in R . Let q ( x, t ) denote theprobability density that X ( t ) = x and assume that itsatisfies the fractional diffusion equation, ∂∂t q = D t L x q, x ∈ R , t > , (15)where D t = D − αt , α ∈ (0 , , (16)is the fractional derivative of Riemann-Liouville typegiven in (4), and L x = K α ∂ ∂x , (17)is the one-dimensional Laplacian with generalized diffu-sivity K α >
0. We use the subscripts in (16) and (17) toemphasize that D t acts only on the time variable t and L x acts only on the space variable x .Langlands et al. [26] developed a mesoscopiccontinuous-time random walk argument to derive the fol-lowing system of fractional reaction-diffusion equations, ∂∂t p = e Rt D t h e − Rt L x p i + R p , x ∈ R , t > , (18)for the joint density p ( x, t ) = ( p i ( x, t )) ni =1 , where p i ( x, t ) d x = P ( X ( t ) = x, J ( t ) = i ) . The derivation in [26] implicitly assumed that X and J are independent processes.We now prove that (18) follows immediately from (12),(15), the independence of X and J , and the linearity of D t and L x . Note first that independence ensures thatthe joint probability distribution is the product of theindividual distributions, p i ( x, t ) d x = P ( X ( t ) = x, J ( t ) = i )= P ( X ( t ) = x ) P ( J ( t ) = i )= q ( x, t ) r i ( t ) d x. (19)Therefore, differentiating p ( x, t ) = q ( x, t ) r ( t ) with re-spect to time and using (12) and (15) yields ∂∂t p ( x, t ) = r ( t ) D t h L x q ( x, t ) i + q ( x, t ) R r ( t )= r ( t ) D t h L x q ( x, t ) i + R p ( x, t ) . (20)Using (14), the first term in the righthand side of (20)becomes r ( t ) D t h L x q ( x, t ) i = e Rt r (0) D t h L x q ( x, t ) i = e Rt D t h r (0) L x q ( x, t ) i = e Rt D t h e − Rt r ( t ) L x q ( x, t ) i = e Rt D t h e − Rt L x p ( x, t ) i . (21)Combining (20) and (21) yields (18). III. MORE GENERAL SETTING
It is easy to see that the calculation in (19)-(21) and theresulting evolution equation in (18) holds in much greatergenerality. First, the spatial domain need not be R , andwe will instead take it to be any d -dimensional open set V ⊆ R d with d ≥
1. Second, the operator L x need notbe the Laplacian and the operator D t need not be theRiemann-Liouville fractional derivative. Instead, we willtake L x to be any linear operator acting on functions ofspace x ∈ V ⊆ R d and D t to be any linear operator actingon functions of time t ∈ [0 , ∞ ). That is, if ϕ ( t ) , ψ ( t ) arereal-valued functions of time t ∈ [0 , ∞ ) in the domainof D t and f ( x ) , g ( x ) are real-valued functions of space x ∈ V ⊆ R d in the domain of L x , then we assume L x ( ϕf + ψg ) = ϕ L x f + ψ L x g, D t ( ϕf + ψg ) = f D t ϕ + g D t ψ. (22)Third, the jump process J ( t ) need not have constantjump rates or a finite state space. We summarize thisin the following theorem. Equation (27) in Theorem 1and its proof is the main result of this paper. Theorem 1.
Assume { J ( t ) } t ≥ is a stochastic pro-cess on the possibly infinite, countable state space, { , , . . . , n } , where n ∈ N ∪ {∞} . Suppose the distribution, r ( t ) := ( r i ( t )) ni =1 := ( P ( J ( t ) = i )) ni =1 ∈ R n , satisfies dd t r ( t ) = R ( t ) r ( t ) , t > , (23) for some function R ( t ) : [0 , ∞ ) R n × n , and r ( t ) = Ψ( t ) r (0) , t ≥ , (24) where Ψ( t ) : ( −∞ , ∞ ) R n × n satisfies Ψ( t )Ψ( − t ) = id , t ∈ ( −∞ , ∞ ) , (25) where id is the identity operator.Assume { X ( t ) } t ≥ is a stochastic process taking valuesin the closure of the open set V ⊆ R d whose probabilitydensity, q ( x, t ) d x = P ( X ( t ) = x ) , satisfies ∂∂t q = D t L x q, x ∈ V ⊆ R d , t > , (26) where L x and D t are linear operators satisfying (22) . If X and J are independent, then the joint probabilitydensity p ( x, t ) = ( p i ( x, t )) ni =1 , p i ( x, t ) d x = P ( X ( t ) = x, J ( t ) = i ) , satisfies ∂∂t p = Ψ( t ) D t h Ψ( − t ) L x p i + R ( t ) p , x ∈ V, t > . (27) Proof of Theorem 1.
Since X and J are independent, thejoint probability density is simply the product, p ( x, t ) = q ( x, t ) r ( t ) , and the proof then follows exactly as in (20)-(21) with e ± Rt replaced by Ψ( ± t ).Theorem 1 states that if X and J are independent,then their joint density p ( x, t ) satisfies the evolutionequations in (27). To investigate the converse of The-orem 1, assume that the joint density p ( x, t ) of X and J satisfies the evolution equations in (27). Now, notice thatthe product q ( x, t ) r ( t ) also satisfies (27) if q ( x, t ) satisfies(26) and r ( t ) satisfies (23). Therefore, if (i) p ( x, t ) and q ( x, t ) r ( t ) satisfy the same initial conditions and bound-ary conditions (or growth conditions if the domain V is unbounded) and if (ii) the solution to equation (27)with these initial/boundary conditions is unique, then p ( x, t ) = q ( x, t ) r ( t ). Therefore, X and J must be in-dependent. In conclusion, the joint density of X and J satisfies (27) if and only if X and J are independent, aslong as dependencies between X and J are not imposedat t = 0 or on the spatial boundary. IV. EXAMPLES
In this section, we illustrate Theorem 1 by applying itto some examples of interest.
A. Some previous results
To get the result (11) of Langlands et al. in [26], thenwe apply Theorem 1 with V = R , L x = K α ∂ ∂x , D t = D − αt , Ψ( t ) = e Rt , where R = R ( t ) is constant in time and n < ∞ . B. Fractional Fokker-Planck equations
To find the evolution equations for fractional Fokker-Planck equations with linear reactions, we apply Theo-rem 1 with L x given by the Fokker-Planck operator in(6) and D t = D − αt . C. Other memory kernels
Theorem 1 shows that the form of the evolution equa-tions in (27) holds for more general operators than theRiemann-Liouville fractional derivative D − αt . For ex-ample, we can take the time operator D t to be theintegro-differential operator [28–30], D t ϕ ( t ) = dd t Z t M ( t − t ′ ) ϕ ( t ′ ) d t ′ , (28)where M ( t ) : [0 , ∞ ) R is the so-called memory kernel. D. Superdiffusion
Using a continuous-time random walk argument andproperties of Fourier-Laplace transforms, Schmidt et al.[25] found that the reaction and movement terms are de-coupled in the reaction-superdiffusion equations for L´evyflights with a single irreversible reaction. In this case, theequation describing the movement (without reaction) ofa single molecule is ∂∂t q = K µ ∆ µ/ x q, where ∆ µ/ x is the Riesz symmetric fractional derivativeacting on x . Therefore, the decoupling of reaction andmovement terms in the reaction-superdiffusion equationsfollows from Theorem 1 upon taking the spatial operatorto be L x = K µ ∆ µ/ x and the time operator D t to be theidentity. E. Time-dependent rates
Theorem 1 allows the reaction rate matrix to vary intime. In particular, suppose that the reaction rate matrixis some given function of time { R ( t ) } t ≥ (which does notdepend on spatial position). Starting from the result ofLanglands et al. in (11) for constant reaction rates, onemight conjecture that the evolution equations for suchtime-dependent reaction rates are ∂∂t p = e R t R ( s ) d s D t h e − R t R ( s ) d s L x p i + R ( t ) p , (29)where the integration R t R ( s ) d s is performed componentwise. Indeed, (29) has been used to model some physi-cal systems involving a single irreversible reaction with atime-dependent rate [31, 32]. However, Theorem 1 showsthat the conjecture in (29) can fail, since the solution op-erator Ψ( ± t ) in (24) for the equation (23) is not alwaysgiven by the matrix exponential e ± R t R ( s ) d s .In fact, the two-state irreversible reaction,1 λ ( t ) → , (30)is a rare case of time-dependent reaction rates for which(29) holds, since this is one of the few instances of time-dependent reaction rates in which Ψ( ± t ) = e ± R t R ( s ) d s (see Appendix A.2.4 in [33]). To illustrate, suppose J ( t ) ∈ { , } models (30) for some reaction rate λ ( t ),and thus assume that the distribution r ( t ) ∈ R satisfiesthe nonautonomous linear system of ordinary differen-tial equations in (23) with time-dependent reaction ratematrix, R ( t ) = (cid:18) − λ ( t ) 0 λ ( t ) 0 (cid:19) ∈ R × . In this case, one can check that the solution operator in(24) is indeed the matrix exponential,Ψ( t ) = e R t R ( s ) d s = e − R t λ ( s ) d s − e − R t λ ( s ) d s ! , for t ≥ , and Ψ( − t ) = e − R t R ( s ) d s for t > ± t ) = e ± R t R ( s ) d s , and thus Theorem 1shows that (27) holds rather than (29). For example,suppose (30) is now reversible,1 λ ( t ) ⇋ λ ( t ) , and thus r ( t ) ∈ R satisfies (23) with R ( t ) = (cid:18) − λ ( t ) λ ( t ) λ ( t ) − λ ( t ) (cid:19) ∈ R × . In this case, it is straightforward to check that the solu-tion operator isΨ( t ) = (cid:18) Ψ ( t ) 1 − Ψ ( t )1 − Ψ ( t ) Ψ ( t ) (cid:19) , t ≥ , (31)where for i ∈ { , } and t ≥ ii ( t ) = e − R t ( λ ( s )+ λ ( s )) d s × (cid:16) Z t λ − i ( s ) e R s ( λ ( σ )+ λ ( σ )) d σ d s (cid:17) . Also, the condition (25) implies that the operator evalu-ated at a negative time argument is the matrix inverseΨ( − t ) = (Ψ( t )) − , for t > . Note that matrix Ψ( t ) is invertible for each t ≥ e R t R ( s ) d s = (cid:18) − χ ( t ) χ ( t ) χ ( t ) 1 − χ ( t ) (cid:19) , t ≥ , (32) where χ ij ( t ) = R t λ j ( s ) d s (cid:16) − e − R t ( λ ( s )+ λ ( s )) d s (cid:17)R t λ ( s ) d s + R t λ ( s ) d s . It is straightforward to check that (31) and (32) are gen-erally not equal, except in special cases (such as constantrates, λ j ( t ) ≡ λ j >
0, or equal rates, λ ( t ) = λ ( t ) for all t ≥ J ( t ) ∈ { , , } has two irreversible reactions,1 λ ( t ) → λ ( t ) → , and its distribution r ( t ) ∈ R satisfies (23) with R ( t ) = − λ ( t ) 0 0 λ ( t ) − λ ( t ) 00 λ ( t ) 0 ∈ R × . The corresponding solution operator for t ≥ t ) = e − R t λ ( s ) d s ( t ) e − R t λ ( s ) d s − e − R t λ ( s ) d s − Ψ ( t ) 1 − e − R t λ ( s ) d s , whereΨ ( t ) = e − R t λ ( s ) d s Z t λ ( s ) e R s ( λ ( σ ) − λ ( σ )) d σ d s. For this example, one can check thatΨ( t ) = e R t R ( s ) d s , if t > , except for special cases, and thus (29) is invalid.Summarizing, except for a single irreversible reaction,the evolution equation (29) is typically false for time dep-pendent rates and is corrected by (27) in Theorem 1. V. DISCUSSION
We have given a short and elementary proof of the evo-lution equations for a general class of systems which cancombine anomalous motion with linear reaction kinetics.Our results generalize some previous results in [23–26].The derivations of these previous results employed a va-riety of mathematical techniques, including continuous-time random walk theory, Fourier and Laplace trans-forms, Tauberian theorems, and asymptotic expansions.In light of these previous derivations, one might concludethat the form of the evolution equations depends on thesefiner details. However, we have shown that the evolutionequations follow directly from (i) the independence of thestochastic spatial and discrete processes describing a sin-gle particle and (ii) the linearity of the integro-differentialoperators describing particle motion.Of course, in the present work and in the previous work[23–26], the evolution equations are not strictly necessaryin the sense that the solution to the equations is merelythe product of the distributions of the spatial and discreteprocesses. Nevertheless, these results are expected to beuseful for developing models where the independence as-sumption breaks down. Indeed, evolution equations of avery similar form to (27) have been derived in [34, 35] forpure subdiffusion with certain space-dependent reactionrates. Furthermore, we agree with Refs. [23, 26] thatthese results could provide a platform for investigatingnonlinear reactions, such as those stemming from mass-action kinetics.For example, a natural starting place is an irreversiblebimolecular reaction of the form [36] A + A → ∅ , which describes particles that can annihilate each other.However, while the general form of the evolution equa-tions in (27) may be instructive for this nonlinear exam- ple, it is clear that the approach of the present work can-not be applied directly. Indeed, the present work reliedon the independence of the spatial position and discretestate of a single particle. However, it is clear for thisexample that a single particle is more likely to be in thediscrete “annihilated” state if it is in a region of spacecontaining a high concentration of particles. Similarly, ifwe consider a unimolecular reaction of the form A λ ( x ) → ∅ , where the first order rate λ ( x ) > x of the particle, it is clear that the particle ismore likely to be in the “annihilated” state if it is in aregion of space where λ ( x ) is large. ACKNOWLEDGMENTS
The author gratefully acknowledges support fromthe National Science Foundation (DMS-1944574, DMS-1814832, and DMS-1148230). [1] H. Scher and E. W. Montroll, Physical Review B , 2455(1975).[2] B. Berkowitz, J. Klafter, R. Metzler, and H. Scher, Wa-ter Resources Research , 9 (2002).[3] F. Amblard, A. C. Maggs, B. Yurke, A. N. Pargellis, andS. Leibler, Physical review letters , 4470 (1996).[4] I. Golding and E. C. Cox, Physical review letters ,098102 (2006).[5] F. H¨ofling and T. Franosch, Reports on Progress inPhysics , 046602 (2013).[6] J. Klafter and I. M. Sokolov, Physics world , 29 (2005).[7] A. Caspi, R. Granek, and M. Elbaum, Physical ReviewLetters , 5655 (2000).[8] R. Metzler and J. Klafter, Physics reports , 1 (2000).[9] S. G. Samko, A. A. Kilbas, O. I. Marichev, et al. , Fractional integrals and derivatives , Vol. 1 (Gordon andBreach Science Publishers, Yverdon Yverdon-les-Bains,Switzerland, 1993).[10] I. Podlubny,
Fractional differential equations: an intro-duction to fractional derivatives, fractional differentialequations, to methods of their solution and some of theirapplications (Elsevier, 1998).[11] R. Metzler, E. Barkai, and J. Klafter, Physical reviewletters , 3563 (1999).[12] A. Nepomnyashchy, Mathematical Modelling of NaturalPhenomena , 26 (2016).[13] G. Hornung, B. Berkowitz, and N. Barkai, Physical Re-view E , 041916 (2005).[14] V. Gafiychuk and B. Datsko, Physical Review E ,066210 (2008).[15] J. P. Boon, J. F. Lutsko, and C. Lutsko, Physical ReviewE , 021126 (2012).[16] C. N. Angstmann, I. C. Donnelly, and B. I. Henry, Math-ematical Modelling of Natural Phenomena , 17 (2013).[17] T. Koszto lowicz and K. Lewandowska, Mathematical Modelling of Natural Phenomena , 44 (2013).[18] S. K. Hansen and B. Berkowitz, Physical Review E ,032113 (2015).[19] P. Straka and S. Fedotov, Journal of theoretical biology , 71 (2015).[20] M. A. Dos Santos, Journal of Statistical Mechanics: The-ory and Experiment , 033214 (2019).[21] W.-B. Zhang and M. Yi, Physica A: Statistical Mechanicsand its Applications , 121347 (2019).[22] G. Li, H. Zhang, and B. Zhang, Physica A: StatisticalMechanics and its Applications , 121917 (2019).[23] B. Henry, T. Langlands, and S. Wearne, Physical ReviewE , 031116 (2006).[24] I. M. Sokolov, M. Schmidt, and F. Sagu´es, Physical Re-view E , 031102 (2006).[25] M. Schmidt, F. Sagu´es, and I. Sokolov, Journal ofPhysics: Condensed Matter , 065118 (2007).[26] T. Langlands, B. I. Henry, and S. L. Wearne, PhysicalReview E , 021111 (2008).[27] J. Norris, Markov Chains , Statistical & ProbabilisticMathematics (Cambridge University Press, 1998).[28] I. M. Sokolov and J. Klafter, Physical review letters ,140602 (2006).[29] M. Magdziarz, Journal of Statistical Physics , 763(2009).[30] S. Carnaffan and R. Kawai, SIAM Journal on ScientificComputing , B886 (2017).[31] E. Abad, S. Yuste, and K. Lindenberg, Physical ReviewE , 061120 (2012).[32] E. Abad, S. Yuste, and K. Lindenberg, MathematicalModelling of Natural Phenomena , 100 (2013).[33] O. Aalen, O. Borgan, and H. Gjessing, Survival andevent history analysis: a process point of view (SpringerScience & Business Media, 2008).[34] E. Abad, S. Yuste, and K. Lindenberg, Physical Review E , 031115 (2010).[35] S. Yuste, E. Abad, and K. Lindenberg, Journal of Statis-tical Mechanics: Theory and Experiment , P11014(2014). [36] S. B. Yuste and K. Lindenberg, Physical Review Letters87