Stochastic Simulation of Fractional Dynamics

ABSTRACT

Theoretical formalisms based on fractional differential equations are increasingly being used to model complex kinetic processes, as they can accurately capture memory effects and temporal heterogeneity. This study proposes a stochastic simulation methodology for processes governed by such equations. The numerical implementation uses a discrete-time Monte Carlo algorithm to model the evolution of an irreversible unimolecular reaction by integrating Mittag-Leffler statistics directly into the dynamics of microscopic events. The study concludes that the proposed method is robust and consistently recovers the classical limit, thereby establishing a rigorous link between fractional chemical kinetics and computational simulation. Consequently, this approach constitutes a powerful predictive framework for chemical systems in complex media.

Keywords:
fractional kinetics; time-driven Monte Carlo methods; Mittag-Leffler distribution

1 INTRODUCTION

Traditionally, the theoretical description of chemical kinetics has relied on deterministic rate equations formulated in terms of integer-order ordinary differential equations, which are derived from the law of mass action. These equations implicitly assume homogeneous conditions and large particle numbers. While this framework has been remarkably successful in describing a wide range of homogeneous reactions, it fundamentally neglects important aspects of molecular dynamics, such as spatial heterogeneity and memory effects ³.

An alternative approach to addressing this type of problem is stochastic simulation ⁴⁰. The stochastic nature of chemical reactions was rigorously formalized through the Chemical Master Equation (CME) ^{4), (24}, which describes the time evolution of the probability distribution over discrete molecular states. In the 1970s, Gillespie’s pioneering work provided an exact computational realization of the CME through the Stochastic Simulation Algorithm (SSA), also known as the event-driven Monte Carlo simulation ^{11), (12), (13}.

Unlike deterministic models based on ordinary differential equations, the SSA explicitly considers the probabilistic nature of chemical reactions, computing the time to the next reaction event at each step. This approach has established a paradigm for linking microscopic molecular events to macroscopic kinetic behaviour, and has become a cornerstone of computational chemical physics. In this regime, the temporal evolution of the system must be described by a stochastic Markov process. Classical chemical kinetics, which is based on ordinary differential equations, describes the evolution of average concentrations. Although numerous variants have been proposed to improve efficiency and address multiscale complexity, the foundational Markovian assumption remains central to the method.

According to the Markovian assumption, the probability of a reaction occurring within the next infinitesimal time interval depends solely on the system’s current state, independent of its past history. This memoryless property is encoded mathematically in the Poisson statistics of reaction events, which is reflected in the exponential waiting-time distribution underlying the classical SSA. At a microscopic level, this assumption is valid when molecular collisions arise from uncorrelated diffusive encounters in a uniform medium.

However, a growing body of experimental evidence demonstrates systematic deviations from this idealized behavior across diverse chemical and physical systems. Examples include subdiffusive transport in biological environments ^{29), (36), (43}, non-exponential relaxation in polymer networks ³⁸, anomalous reaction kinetics on heterogeneous catalytic surfaces ^{9), (20}, and history-dependent rate coefficients in viscoelastic media ¹⁵. All of these phenomena exhibit temporal correlations that are incompatible with Markovian dynamics. In these systems, past molecular configurations continue to influence present reaction probabilities through mechanisms such as structural rearrangements, geometric constraints or slow environmental relaxation.

In complex chemical environments like these, the Markovian assumption on which classical kinetic models are based may be inadequate. Such cases often exhibit wide distributions of waiting times between reaction events. In such situations, the exponential waiting-time distribution associated with Markovian dynamics is replaced by heavy-tailed distributions. Within the framework of Continuous Time Random Walk (CTRW) theory ^{33), (37}, waiting-time densities of the form ψ(t) ∼ t ^−1−α , with 0 < α < 1, naturally lead, in the hydrodynamic limit, to evolution equations involving fractional time derivatives ^{3), (19), (33), (37}. These equations replace integer-order time derivatives with non-integer (fractional) derivatives ^{6), (7), (18), (35}. From the perspective of chemical kinetics, the fractional derivative therefore encodes non-local temporal correlations arising from environmental heterogeneity or transport-limited reaction mechanisms.

Recent work by Angstmann, Donnelly, and Henry ^{1), (2), (44} has demonstrated that exact event-driven stochastic algorithms can be constructed for fractional master equations by sampling waiting times from Mittag-Leffler distributions ^{10), (14}, thereby reproducing time-fractional dynamics at the trajectory level. However, extending these methods to fractional kinetics requires replacing exponential clocks with Mittag-Leffler distributed waiting times, which is a nontrivial modification with significant computational and theoretical implications ^{5), (22}.

An alternative approach is to develop a generalization of the stochastic method based on a time-driven Monte Carlo algorithm ^{8), (16), (17), (27), (28}. This method has been used, particularly in simplified kinetic models and pedagogical contexts ^{26), (31), (34), (41}. In these approaches, time is divided into fixed intervals and reaction events are tested probabilistically at each stage. Although these methods are conceptually transparent and straightforward to implement, they are inherently approximate and only reproduce Poissonian kinetics in the limit of a vanishing time step. Despite their methodological differences, both event-driven SSA and time-driven Monte Carlo algorithms belong to the broader class of stochastic simulation techniques based on probabilistic descriptions of microscopic dynamics.

Taking this alternative approach to stochastic simulations into account, this work suggests generalising the time-driven Monte Carlo method in order to simulate processes that are described by fractional differential equations. Our approach builds on classical algorithms by rigorously deriving fractional transition probabilities based on Mittag-Leffler statistics. This provides an alternative to event-driven methods that is conceptually transparent and computationally efficient for moderate system sizes. It is also well-suited to pedagogical exposition ^{25), (42}. To the best of our knowledge, this work is the first to systematically investigate discrete-time Monte Carlo methods for fractional chemical kinetics, providing detailed error analyses, stability criteria and direct comparisons with exact fractional solutions.

The remainder of this paper is organised as follows: Section 2 introduces the theoretical foundations of fractional chemical kinetics, emphasising the connection between the CTRW formalism and time-fractional kinetic equations. This section also presents the generalised fractional Gillespie algorithm and the construction of the time-driven Monte Carlo framework. Section 3 reports on extensive numerical simulations of unimolecular fractional decay processes. These include systematic convergence studies, error analysis, and direct comparisons with exact analytical solutions and the classical Markovian limit. Finally, Section 4 summarises the main results and outlines perspectives for future work.

2 METHODOLOGY

2.1 Continuous-Time Random Walk model

In 1965, the classical article by Montroll and Weiss ³³ introduced the Continuous-Time Random Walk (CTRW) model as a conceptual generalization of the classical random walk. In this formulation, two fundamental assumptions are relaxed: (i) the waiting times between successive events become random variables, and (ii) the spatial jumps and the waiting times are statistically independent. Later, this framework was interpreted as showing that the average dynamics are governed by fractional differential equations rather than integer-order ones ²⁹.

Consider initially an irreversible transition process: A → B. Let ψ(t)dt denote the probability that the transition occurs within the time interval between t and t + dt. The survival function, Ψ(t), defined as the probability that the transition has not yet occurred up to time t, is given by the expression:

$\Psi \left(t\right)={\int }_t^{\infty }\psi \left(t\text{'}\right)dt\text{'}=1-{\int }_0^t\psi \left(t\text{'}\right)dt\text{'}.$ (2.1)

The function Ψ(t) is the survival probability (or complementary cumulative distribution). In CTRW, it is often assumed that ψ(t) ∼ t ^−1−α , and therefore Ψ(t) ∼ t ^−α by integration.

Since the process is irreversible and species A cannot be replenished, the probability of the system remaining in state A at time t is exactly equal to the survival function. Therefore, assuming that the function is normalized, P _A (t) = Ψ(t). Passing to the Laplace domain, the following relation is obtained directly:

${\hat{P}}_A\left(s\right)\hat{\Psi }\left(s\right)=\frac{1-\hat{\Psi }\left(s\right)}{s},$ (2.2)

which expresses the survival probability in terms of the Laplace transform of the waiting-time density.

It is convenient to introduce the memory kernel in the Laplace domain through the following definition

$\hat{K}\left(s\right):=\frac{s\hat{\Psi }\left(s\right)}{1-\hat{\psi }\left(s\right)}.$ (2.3)

Using this definition, the expression for ${\hat{P}}_A\left(s\right)$ can be rewritten in the compact form

${\hat{P}}_A\left(s\right)=\frac{1}{s+\hat{K}\left(s\right)}.$ (2.4)

Assuming the initial condition P _A (0) = 1, this relation is equivalent to the algebraic equation

$s{\hat{P}}_A\left(s\right)-P_A\left(0\right)+\hat{K}\left(s\right){\hat{P}}_A\left(s\right)=0.$ (2.5)

Upon returning to the time domain by inverse Laplace transform, one obtains the integro-differential master equation characteristic of the CTRW formalism,

$\frac{d{P}_A\left(t\right)}{dt}=-{\int }_0^{t}K\left(t-t\text{'}\right)P_A\left(t\text{'}\right)dt\text{'}.$ (2.6)

This equation explicitly expresses the nonlocal character in time of the dynamics, since the rate of change of P _A (t) depends on the entire past history of the system weighted by the memory kernel K(t).

In the classical Markovian limit, the waiting times are exponentially distributed. Considering

$\psi \left(t\right)=k{e}^{-kt},$ (2.7)

one has

$\hat{\psi }\left(s\right)=\frac{k}{k+s}.$ (2.8)

Substituting this expression into the definition of the kernel yields

$\hat{K}\left(s\right)=\frac{s\hat{\psi }\left(s\right)}{1-\hat{\psi }\left(s\right)}=k,$ (2.9)

whose inverse transform is

$K\left(t\right)=k\delta \left(t\right).$ (2.10)

Inserting this result into the master equation immediately gives

$\frac{d{P}_A\left(t\right)}{dt}=-k{\int }_0^t\delta \left(t-t\text{'}\right)P_A\left(t\text{'}\right)dt\text{'}=-k{P}_A\left(t\right),$ (2.11)

which corresponds to the usual first-order kinetic law characteristic of memoryless processes.

Consider now the non-Markovian case, in which the waiting-time density exhibits a long tail. In this context, one assumes an asymptotic form compatible with power-law distributions, characterized in the Laplace domain by

$\hat{\psi }\left(s\right)=\frac{1}{1+{\left(\tau s\right)}^{\alpha }},0<\alpha <1.$ (2.12)

The choice of $\hat{\psi }\left(s\right)$ corresponds in the time domain to a waiting-time statistics governed by a Mittag-Leffler function ¹⁴,

$\psi \left(t\right)=\frac{t^{\alpha -1}}{{\tau }^{\alpha }}{E}_{\alpha ,\alpha }\left(-{\left(t/\tau \right)}^{\alpha }\right).$ (2.13)

The two-parameter Mittag-Leffler function E _α,β (z) is defined by the following power series:

$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\left(\alpha k+\beta \right)},$

where z ∈ ℂ is a complex variable, α and β are complex parameters constrained by Re(α) > 0 for convergence, and Γ(·) denotes the classical Gamma function. In this case, the survival function associated with the waiting-time density ψ(t) is

$\psi \left(t\right)=E_{\alpha ,1}\left(-{\left(t/\tau \right)}^{\alpha }\right),$ (2.14)

with the asymptotic limit for t → ∞ given by

$\psi \left(t\right)~\frac{{\tau }^{\alpha }}{\Gamma \left(1-\alpha \right)}{t}^{-\alpha }.$ (2.15)

From the expression of $\hat{\psi }\left(t\right)$, the memory kernel in the Laplace domain becomes

$\hat{K}\left(s\right):=\frac{s\hat{\psi }\left(s\right)}{1-\hat{\psi }\left(s\right)}={\tau }^{-\alpha }{s}^{1-\alpha }.$ (2.16)

For 0 < α < 1, one has 1 − α ∈ (0, 1). The inverse transform yields the temporal kernel

$K\left(t\right)=\frac{{\tau }^{-\alpha }}{\Gamma \left(\alpha -1\right)}{t}^{\alpha -2},$ (2.17)

which decays slowly in time, reflecting the presence of long-range memory. When 0 < α < 1, one has α − 2 < −1, so that t ^α−2 is singular at t = 0, requiring interpretation in the sense of distributions.

Substituting this kernel into the master equation, one obtains

$\frac{d{P}_A\left(t\right)}{dt}=-\frac{{\tau }^{-\alpha }}{\Gamma \left(\alpha -1\right)}{\int }_0^t{\left(t-t\text{'}\right)}^{\alpha -2}{P}_A\left(t\text{'}\right)dt\text{'}.$ (2.18)

In the Laplace domain, this equation is equivalent to

$s{\hat{P}}_A\left(s\right)-P_A\left(0\right)=-{\tau }^{-\alpha }\frac{{\hat{P}}_A\left(s\right)}{s^{\alpha -1}},$ (2.19)

or, equivalently,

$s^{\alpha }{\hat{P}}_A\left(s\right)-s^{\alpha -1}{P}_A\left(0\right)=-{\tau }^{-\alpha }{\hat{P}}_A\left(s\right).$ (2.20)

In the present work, the temporal memory of the system is modeled through the Caputo fractional derivative, defined as

$D_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\int }_0^t\frac{f\text{'}\left(s\right)}{{\left(t-s\right)}^{\alpha }}ds,0<\alpha <1.$

The Caputo formulation is particularly suitable for physical applications because it allows the use of classical initial conditions expressed in terms of integer-order derivatives ³⁵. In the context of stochastic chemical kinetics, this operator introduces a history-dependent contribution to the reaction dynamics, effectively accounting for memory effects that arise from non-exponential waiting-time statistics. If f(t) is a sufficiently regular function for t ≥ 0, then the Laplace transform of the Caputo derivative is given by

$\mathcal{L}\left\{D_t^{\alpha }f\left(t\right)\right\}\left(s\right)=s^{\alpha }\hat{f}\left(s\right)-s^{\alpha -1}f\left(0\right),$ (2.21)

for 0 ≤ α ≤ 1, where $\hat{f}\left(s\right)=\mathcal{L}\left\{f\left(t\right)\right\}\left(s\right)$ ^{6), (7), (18), (35}. Returning to the time domain from Equation (2.20), one finally obtains the fractional kinetic equation:

$D_t^{\alpha }{P}_A\left(t\right)=-{\tau }^{-\alpha }{P}_A\left(t\right),$ (2.22)

where $D_t^{\alpha }$ denotes the Caputo fractional derivative of order α.

Fractional kinetics therefore does not arise as an ad-hoc phenomenological hypothesis, but rather as a direct mathematical consequence of the fundamental probabilistic assumptions of the CTRW model. In a broader modeling context, this constructive, non-ad-hoc framework has recently been corroborated by Monteiro et al. ³². From a physicochemical standpoint, this derivation provides a microscopic interpretation of the fractional order α as a parameter quantifying the degree of temporal heterogeneity in the underlying stochastic dynamics. Values α < 1 indicate the presence of long-time correlations, which may arise from transport limitations, trapping processes, or structural disorder in the medium. In heterogeneous catalytic systems, for instance, such effects may originate from diffusion through porous structures or adsorption-desorption cycles occurring over a wide range of characteristic times. The classical Markovian limit is recovered when α = 1, corresponding to exponential waiting-time statistics and memoryless reaction dynamics.

2.2 Gillespie algorithm and fractional stochastic simulation method

The Gillespie stochastic simulation algorithm (SSA) ^{11), (12), (13} provides an exact, event-driven Monte Carlo framework for simulating chemically reacting systems under the assumption of Markovian dynamics, namely exponentially distributed waiting times between reaction events. Within this framework, the temporal evolution of the system is constructed by explicitly sampling both the time to the next reaction and the identity of the reaction channel. This yields trajectories that are statistically consistent with the Chemical Master Equation ^{26), (41}.

Consider a set of elementary reactions, R _j , occurring in a homogeneous medium. The corresponding propensity function for each reaction channel is defined as follows:

$a_0=\sum _j{a}_j,$ (2.23)

where a _j = k _j h _j (X), h _j = ∏_s X _s denotes the number of distinct reactant combinations compatible with reaction R _j , and k _j is the associated stochastic rate constant. Here, X = (X ₁ , X ₂ , . . . , X _S )^T represents the state vector of the system, where each component X _s denotes the current population (number of molecules) of the chemical species s. The propensity to react is derived from the probability of a molecular collision, based on statistical mechanics arguments. Two independent random numbers, r ₁ , r ₂, uniformly distributed on U(0, 1), are generated. The waiting time to the next reaction event is then sampled according to the following equation:

$\tau =\frac{1}{a_0}\ln \left(\frac{1}{r_1}\right).$ (2.24)

The probability that the selected reaction channel is R _j is given by a _j /a ₀, and operationally the reaction index j is determined as the smallest integer satisfying

$\sum _{k=1}^j{a}_k>r_2{a}_0.$ (2.25)

The system time is updated as t ← t + τ, and the state vector is modified according to X ← X + ν _j . This procedure is iterated until a prescribed final time is reached. The vector ν _j = (ν _1j , ν _2j , . . . , ν _Sj )^T defines the net changes, where each component ν _sj represents the stoichiometric coefficient of species s when a single reaction event j occurs. In the classical SSA, the exponential distribution of τ directly reflects the Markovian assumption underlying Poisson reaction statistics.

As demonstrated by Angstmann, Donnelly, and Henry ^{1), (44}, a Fractional Stochastic Algorithm can be developed whose stochastic trajectories converge, in the ensemble sense, to the solution of the time-fractional master equation. The key modification relative to the classical Gillespie algorithm lies in the sampling of the waiting times τ, which now follow a Mittag-Leffler probability density:

${\psi }_{\alpha }\left(\tau \right)=a_0{\tau }^{\alpha -1}{E}_{\alpha ,\alpha }\left(-a_0{\tau }^{\alpha }\right),$ (2.26)

where a ₀ denotes the total propensity of all reaction channels available in the current system state. Within the CTRW formalism, these non-exponential waiting times lead, in the macroscopic limit, to time-fractional master equations involving nonlocal-in-time operators, most commonly Caputo fractional derivatives.

Equivalently, introducing the characteristic timescale ${\tau }_s=a_0^{-1/\alpha }$, the survival probability can be written in the rescaled form

${\Psi }_{\alpha }\left(\tau \right)=E_{\alpha ,1}\left[-{\left(\frac{\tau }{{\tau }_s}\right)}^{\alpha }\right].$ (2.27)

The efficient generation of random variables with this distribution is non-trivial, since the Mittag-Leffler function does not, in general, admit an elementary inverse.

From a computational perspective, it is sufficient to generate the dimensionless random variable:

$X=\frac{\tau }{{\tau }_s},$ (2.28)

and then recover the physical waiting time via τ = τ _s X. A fundamental result by Kozubowski and Rachev ^{21), (22} shows that the random variable X, characterized by the survival function

$Pr\left(X>x\right)=E_{\alpha ,1}\left(-x^{\alpha }\right),$ (2.29)

admits the exact stochastic representation

$X=-\ln \left(U\right){\left(\frac{\sin \left(\alpha \pi \right)}{\tan \left(\alpha \pi V\right)}-\cos \left(\alpha \pi \right)\right)}^{1/\alpha },$ (2.30)

where U and V are independent uniform random variables with uniform distribution on the interval [0, 1]. This closed-form transformation provides an exact sampling procedure, avoiding the need for numerical approximation of the Mittag-Leffler function.

The generalized Gillespie algorithm for irreversible fractional processes of the form $A\stackrel{k}{\to }B$ is constructed as follows: At the initial time t = 0, the system contains N _A0 molecules of species A. Unlike the classical Gillespie algorithm, no single global rate is associated with the system. Instead, each molecule is treated as an independent stochastic entity endowed with its own internal clock. Each molecule j is assigned a random waiting time τ _j , drawn independently from a Mittag-Leffler distribution with survival function ^{30), (39}

${\Psi }_{\alpha }\left(t\right)=E_{\alpha ,1}\left(-k^{\alpha }{t}^{\alpha }\right),$ (2.31)

where 0 < α < 1.

Once the waiting times {τ _j } have been assigned, the next reaction event is determined by a competition between clocks. The time increment to the next event is given by:

$\Delta t=\underset{j=1,\dots ,N_A}{min}\left\{{\tau }_j\right\},$ (2.32)

The molecule associated with this minimum is the next to react. The system time is updated to t ← t + ∆t, and the reacting molecule is removed from the population: N _A ← N _A − 1. The residual waiting times of all remaining molecules that have not yet reacted are updated according to τ _i ← τ _i − ∆t, where i ≠ j. This update step is essential, as it preserves the memory inherent in the non-Markovian dynamics. Unlike in the classical SSA, the waiting-time distribution is not reset after each event.

These steps are repeated until no more molecules of type A remain in the system. Because the reaction is unimolecular, individual molecules undergo transitions independently. Consequently, the survival probability of a single molecule governed by this heavy-tailed waiting-time distribution is given exactly by the Mittag-Leffler function. Furthermore, under the assumption that statistical fluctuations become negligible when the number of independent realizations in an ensemble is sufficiently large, averaging over many independent realizations yields:

(2.33)

which satisfies the fractional derivative equation

(2.34)

As α approaches 1, the Mittag-Leffler distribution converges to an exponential distribution, residual waiting times lose memory and the algorithm reverts to the classical Gillespie scheme. However, for 0 < α < 1, the dynamics are intrinsically non-Markovian, and no reformulation based on global propensities is possible without compromising theoretical consistency.

2.3 Time-driven Monte Carlo approach and proposed generalized framework

The Monte Carlo method emerged as a general computational paradigm in the 1940s, in the context of the Manhattan Project. Mathematicians such as Stanislaw Ulam and John von Neumann advocated the systematic use of random sampling to solve problems in neutron transport, diffusion and stochastic dynamics that could not be solved using traditional methods. This line of development was consolidated through the contributions of Nicholas Metropolis and Enrico Fermi, with the analogy to games of chance leading to the name ’Monte Carlo’. The first widely recognised publication to introduce the method under this name was published in 1949 by Metropolis and Ulam ^{8), (16), (17), (27), (28}. Since then, the term has come to denote a broad family of stochastic simulation techniques that are widely used in physics, chemistry and physical chemistry to study equilibrium properties, diffusion phenomena, elementary reactions and decay processes.

In the field of chemical kinetics, McQuarrie et al. provided a systematic formulation of stochastic reaction dynamics based on Poisson processes ^{26), (41}, while also discussing temporal discretisation as a practical approximation. This work laid the theoretical foundations for the time-driven Monte Carlo method, which is now widely used. In this approach, time is divided into equal intervals, and the state of the system is updated probabilistically at each step. Each entity (e.g. molecule, particle, radioactive nucleus or reactive site) is assigned a probability p, which is typically proportional to a rate constant λ , of undergoing an event within a given time interval. For each entity, a uniform random number r in the range (0, 1) is generated; if r < p, the event occurs; otherwise, the entity remains unchanged during that step. After all entities have been evaluated, the system state is updated and time advances by ∆t. This makes the algorithm intrinsically time-driven, as opposed to the event-driven nature of the Gillespie method.

This discretised strategy exhibits well-known structural limitations. The introduction of a finite time step introduces an approximation error, and for the method to be valid, it is necessary that λ∆t is small; otherwise, the possibility of multiple events occurring within a single interval is neglected. The method also becomes inefficient when events are rare, since many time steps result in no change to the system. Despite these limitations, discrete-time Monte Carlo methods remain important historically, transparent conceptually, and useful practically, particularly for exploratory studies, pedagogical purposes, or situations in which temporal exactness is not essential.

Now consider a unimolecular process governed by the fractional differential equation ^{23), (45}

$D_t^{\alpha }N\left(t\right)=-\lambda N\left(t\right),$ (2.35)

where $D_t^{\alpha }$ is the Caputo fractional derivative and λ is a constant. We consider the associated initial value problem with n − 1 < α ≤ n and initial conditions [D ^l N(t)]_t=0 = b _l , with l = 1, 2, . . . , n.

Applying the Laplace transform to both sides of Equation (2.35) and using the standard property of the Caputo derivative ^{6), (7), (18), (35},

$\mathcal{L}\left[\left(D_t^{\alpha }N\right)\left(t\right)\right]=s^{\alpha }L\left[N\left(t\right)\right]-\sum _{k=0}^{n-1}{s}^{\alpha -k-1}{\left[D^{k}N\left(t\right)\right]}_{t=0},$ (2.36)

one obtains, after straightforward algebraic manipulations, the equation:

$\mathcal{L}\left[N\left(t\right)\right]=\sum _{l=1}^n\frac{b_l{s}^{\alpha -l}}{s^{\alpha }+\lambda },$ (2.37)

with l = k + 1 and b _l = [D ^l N(t)]_t=0 .

Each term of Equation (2.37) is inverted by the two-parameter Mittag-Leffler function E _α,β (−λt ^α ), defined by ^{10), (14}

$\mathcal{L}\left[t^{\beta -1}{E}_{\alpha ,\beta }\left(-\lambda t^{\alpha }\right)\right]=\frac{s^{\alpha -\beta }}{s^{\alpha }+\lambda }.$ (2.38)

The solution of Equation (2.35) is therefore

$N\left(t\right)=\sum _{l=1}^n{b}_l{t}^{l-1}{E}_{\alpha ,l}\left(-\lambda t^{\alpha }\right).$ (2.39)

For 0 < α ≤ 1, one has n = 1 and b ₁ = N(0), yielding

$N\left(t\right)=N\left(0\right)E_{\alpha ,1}\left(-\lambda t^{\alpha }\right),$ (2.40)

which recovers the standard exponential decay in the limit α = 1.

In the present work, the time-driven Monte Carlo framework is adapted to dynamics governed by the fractional equation (2.35).

For a unidirectional elementary reaction A → B, each molecule of A is assigned a probability of conversion per unit time. In the classical Poisson case, this probability is given by P = k∆t over a sufficiently small interval ∆t. Under fractional dynamics, the probability P = dN/N is obtained by differentiating the analytical solution,

$dN=-\lambda \alpha t^{\alpha -1}dt\frac{d{E}_{\alpha ,1}\left(x\right)}{dx},$ (2.41)

with x = −λt ^α and

$\frac{d{E}_{\alpha ,1}\left(x\right)}{dx}=\frac{E_{\alpha ,0}\left(x\right)}{\alpha x}.$ (2.42)

This leads to the effective transition probability

$P\left(t\right)=\frac{dN}{N}=t^{-1}dt\left[\frac{E_{\alpha ,0}\left(-\lambda t^{\alpha }\right)}{E_{\alpha ,1}\left(-\lambda t^{\alpha }\right)}\right],$ (2.43)

which naturally reduces to P = λ dt when α = 1. In the suggested implementation, the singularity is circumvented by avoiding calculations at t = 0.

The Ehrenfest urn model, introduced by Paul and Tatiana Ehrenfest in 1907, provides a paradigmatic discrete stochastic representation of relaxation towards equilibrium ⁸. In its simplest formulation, a particle is randomly selected at each time step and transferred between two urns, generating a discrete Markov chain whose state corresponds to the number of particles in one urn. In the context of chemical kinetics, the urns can be interpreted as chemical species and the particles as molecules, with each transfer representing a stochastic reaction event.

Figure 1a illustrates the effective transition probability defined by Equation (2.43) for the fractional order α = 0.8 and compares it with the standard result corresponding to the classical case α = 1. For α = 1, Eq. (2.43) reduces identically to the constant Poissonian rate P = λdt, reflecting the memoryless character of Markovian first-order kinetics and the exponential distribution of inter-event times. In contrast, the fractional case α = 0.8 exhibits a time-dependent transition probability that decays algebraically with time and directly reflects the non-Markovian nature of the underlying dynamics.

Figure 1

In the Appendix, the pseudocode for the proposed algorithm is presented.

3 RESULTS AND DISCUSSION

In this section, we present the implementation of the stochastic simulation of the kinetic process described by a fractional differential equation using event-driven and time-driven Monte Carlo schemes, as described in Section 2. The numerical results are compared with the exact analytical solution. As a case study, we consider the irreversible process $A\stackrel{k}{\to }B$.

We begin by presenting the results obtained using the generalized fractional Gillespie algorithm, as presented in ^{11), (12), (13}. The implementation of this method requires an efficient procedure for generating random numbers distributed according to Mittag-Leffler statistics. As discussed in Section 2, we employ the algorithm proposed by Kozubowski and Rachev ²² for this purpose. Figure 1b compares the empirical survival function obtained from N = 10⁶ independent samples using the Kozubowski and Rachev algorithm with the exact theoretical expression, Equation (2.14). The excellent agreement between the two curves provides strong numerical evidence of the correctness of the sampling procedure adopted for generating waiting times with a Mittag-Leffler distribution.

Figure 2a, obtained from the numerical experiment, compares three distinct dynamical descriptions for the irreversible reaction: (i) the ensemble-averaged result computed using the generalized fractional Gillespie algorithm (circle symbols with error bars) with α = 0.8, (ii) the exact analytical solution of the fractional kinetic equation, expressed in terms of the Mittag-Leffler function (continuous line) with α = 0.8, and (iii) the classical Markovian solution characterized by exponential decay (dashed line). The simulation parameters are k = 0.1 and initial population N _A0 = 1000, with the ensemble average computed from M = 50 independent realizations.

Figure 2

Figure 2a shows the excellent agreement between the Gillespie ensemble average, ⟨N _A (t)⟩, and the exact theoretical prediction, Equation (2.40), throughout the entire time interval considered. Details of this algorithm and analysis of the convergence can be found in references ^{11), (12), (13}.

We now present the results obtained using the generalized time-driven Monte Carlo method applied to the same problem.

The distinction between the two methods, Gillespie and time-driven Monte Carlo, lies primarily in how time is treated and how the system transitions are selected. In the time-driven Monte Carlo approach, we adopt the standard Metropolis-style algorithm, in which time is measured in discrete units. At each step, a particle is selected randomly and an attempted transition is performed. Even if the transition is rejected, the “time” variable still advances by one unit.

Figure 2b shows the results of a Monte Carlo simulation of fractional unimolecular decay corresponding to a single stochastic trajectory, compared with the analytical solutions based on the Mittag-Leffler function (α = 0.8) and with the classical exponential limit (α = 1). As expected for a single-trajectory simulation, the stochastic process presents a highly noisy character (see Figure 3a). A clear qualitative agreement is observed in Figure 2b between the simulated Monte Carlo trajectory and the analytical fractional solution. At short times, the decay is relatively rapid and remains close to exponential behavior, reflecting the fact that memory effects do not yet dominate the dynamics. As time evolves, the fractional solution progressively deviates from the classical exponential solution, exhibiting a slower decay that unequivocally highlights the physical impact of long-range memory incorporated through the fractional order α < 1.

Figure 3

The fluctuations visible in the Monte Carlo simulation points are statistically expected, given the intrinsically noisy nature of the stochastic process and the fact that a single trajectory is considered. Figure 3b presents the average over an ensemble of M = 50 independent trajectories of the fractional unimolecular decay process obtained by time-driven Monte Carlo simulation, together with error bars corresponding to the standard deviation σ _N (t). These results are compared with the analytical fractional solution (α = 0.8) based on the Mittag-Leffler function and with the classical exponential solution (α = 1). This comparison constitutes a test of the consistency of the stochastic algorithm in the non-Markovian regime.

3.1 Numerical Convergence with Respect to the Time Step (∆t → 0)

The evaluation of numerical convergence with respect to the time step ∆t is fundamental to ensure that the discrete solution asymptotically approaches the continuous limit of the underlying fractional or Markovian kinetic model as ∆t → 0. A numerical scheme is considered convergent if the global discretization error ε _M (∆t) vanishes in the limit of infinitesimal time steps:

$\underset{\Delta t\to 0}{\lim }{\epsilon }_M\left(\Delta t\right)=0.$ (3.1)

In this study, the discretization error is quantified by the maximum absolute deviation over the observation interval [0, T]:

(3.2)

where ⟨N ^(m) (t)⟩_M represents the ensemble average obtained from M Monte Carlo realizations and N _exact(t) denotes the analytical solution of the corresponding fractional differential equation.

For a numerical method of order p, the asymptotic behavior of the error follows the power law ε _M (∆t) ∼ C(∆t)^p , where C is a system-dependent constant independent of the discretization parameter. In a logarithmic scale, this relationship linearizes as:

$\log {\epsilon }_M=\log C+p\log \Delta t.$ (3.3)

Consequently, the slope of the linear fit in a log-log plot provides the empirical convergence order p. This analysis is of particular importance for non-Markovian dynamics governed by Continuous Time Random Walks (CTRW), where the nonlocal temporal coupling and the heavy-tailed nature of the waiting-time distributions often lead to slower convergence rates compared to classical Markovian kinetics (α = 1).

Table 1 summarizes the convergence results for the classical regime (α = 1) and the subdiffusive fractional regime (α = 0.8). Both cases exhibit a monotonic decay of the error as ∆t decreases, confirming the internal consistency and stability of the proposed stochastic algorithm. Linear regression of the log-log data yields an empirical convergence order of p ≈ 1 (R ² = 0.9899) for the Markovian case, whereas for the fractional case (α = 0.8), the order reduces to p ≈ 0.565 (R ² = 0.9986). The reduction in the empirical convergence order indicates a deceleration in the algorithm’s convergence rate.

Table 1:
Global discretization error ε _M (∆t) for the Markovian (α = 1) and non-Markovian (α = 0.8) regimes.

∆t	0.5000	0.2000	0.1000	0.0500	0.0100
α = 1	9.4220	3.4006	2.2987	1.4323	0.9619
α = 0.8	16.6112	10.1059	7.0030	4.3546	1.8596

This occurs because the heavy-tailed probability density function ψ(t) ∼ t ^−1−α implies a wider divergence in realization behaviors; while some experience extended periods of stasis without reactive events, others undergo rapid successive steps. Consequently, this significantly increases the ensemble variance σ ²(t) compared to the Poissonian statistics characteristic of the α = 1 case. For α = 0.8, the system explores the state space more slowly, requiring a larger ensemble M to adequately sample the tails of the Mittag-Leffler distribution. This point suggests that standard Monte Carlo sampling becomes less efficient as α decreases.

3.2 Statistical Convergence with Respect to the Ensemble Size (M → ∞)

In the Monte Carlo formulation of stochastic chemical kinetics, the physically relevant observable is the ensemble-averaged trajectory rather than an individual realization. Let N ^(m) (t) denote the m-th independent stochastic realization of the process. The empirical ensemble mean is defined as the estimator:

(3.4)

For a fixed time t, the set ${\left\{N^{\left(m\right)}\left(t\right)\right\}}_{m=1}^M$ consists of independent and identically distributed (i.i.d.) random variables with finite variance σ ²(t). According to the Strong Law of Large Numbers, the empirical mean converges to the theoretical expectation:

(3.5)

which guarantees the statistical consistency of the Monte Carlo estimator. Furthermore, the Central Limit Theorem (CLT) dictates the asymptotic distribution of the estimation error:

(3.6)

Consequently, the statistical fluctuation of the ensemble average scales as 𝒪(M ^−1/2 ). This M ^−1/2 scaling is an intrinsic property of Monte Carlo estimators, remaining theoretically independent of the fractional order α or the temporal discretization step ∆t. Increasing M reduces the statistical noise (variance) of the estimator but does not mitigate the systematic discretization error (bias) associated with the finite time-step ∆t. Therefore, a slope near −1/2 in a log-log plot of the error versus M provides quantitative evidence of robust statistical convergence.

Table 2 presents the evolution of the maximum absolute error,

(3.7)

as a function of the ensemble size M for the Markovian (α = 1) and non-Markovian (α = 0.8) regimes. The data present an overall downward trend, subject to stochastic fluctuations intrinsic to the process. Empirical fits yield convergence slopes of β ≈ 0.619 (R ² = 0.9678) for α = 1 and β ≈ 0.3421 (R ² = 0.7766) for α = 0.8. The deviation from the ideal 0.5 slope in the fractional case, coupled with the lower R ² value, reflects the higher intrinsic variance and the ”heavy-tailed” influence of the waiting-time distributions on the sampling efficiency.

Table 2:
Statistical convergence parameter ε _∆t (M) as a function of the ensemble size M for α = 1 and α = 0.8.

M	10	20	50	100	250	500	750	1000	1500	2000
α = 1	12.4449	9.5330	3.7934	3.5878	2.7735	1.1470	1.0130	0.9619	0.5232	0.3995
α = 0.8	15.9084	6.9266	3.1386	2.2092	2.4035	2.9014	1.9800	1.8596	1.8140	1.5789

3.3 Structural Convergence in the Limit α → 1

A fundamental consistency requirement for any fractional kinetic model is the recovery of classical Markovian dynamics as the anomalous exponent α approaches unity. For a unimolecular irreversible reaction A → B, the fractional framework should collapse into the first-order rate equation, where the Caputo fractional derivative reduces to a standard temporal derivative, yielding the exponential survival probability:

$P\left(t\right)=exp\left(-\lambda t\right).$ (3.8)

In the fractional regime, the analytical solution is governed by the Mittag-Leffler function:

$P_{\alpha }\left(t\right)=E_{\alpha }\left(-\lambda t^{\alpha }\right).$ (3.9)

It is mathematically established that $li{m}_{\alpha \to 1}-E_{\alpha }\left(-\lambda t^{\alpha }\right)=exp\left(-\lambda t\right)$ pointwise for t ≥ 0. However, for a rigorous validation of the numerical scheme, pointwise convergence is insufficient. We define a stronger criterion based on the uniform norm over the observation interval [0, T]:

$\mu \left(\alpha \right)=\underset{t\in \left[0,T\right]}{sup}\left|P_{\alpha }\left(t\right)-exp\left(-\lambda t\right)\right|.$ (3.10)

Structural consistency is verified if lim_α→1 − µ(α) = 0, ensuring that the maximum deviation between the non-Markovian and classical trajectories vanishes globally. Physically, this limit signifies the extinction of long-memory effects. As α → 1, the heavy-tailed waiting-time distribution characteristic of Continuous Time Random Walks (CTRW) loses its power-law scaling and converges to an exponential distribution. This transition restores the Markov property, transforming the nonlocal temporal operator into a local one.

Table 3 quantifies this structural deviation. The results demonstrate a strictly monotonic decay of µ(α) as α increases from 0.4 to 1.0. The residual value at α = 1.0 (µ ≈ 3.6647) represents the combined contribution of the finite temporal discretization (∆t = 0.1) and the statistical variance of the ensemble (M = 50). Despite these numerical constraints, the clear trend provides robust evidence that the fractional model is a legitimate generalization of classical kinetics, smoothly recovering the exponential regime.

Table 3:
Structural convergence toward the Markovian limit (α → 1) quantified by the uniform norm µ(α). Parameters: λ = 0.1, M = 50, ∆t = 0.1.

α	0.4	0.5	0.6	0.7	0.8	0.9	α → 1.0
µ(α)	313.5582	257.4353	199.4208	150.9512	101.3750	52.4843	3.6647

3.4 Computational Cost Analysis

A comparative analysis of computational performance reveals a significant disparity in efficiency between the two methodologies. Utilizing a 13th Gen Intel(R) Core(TM) i7-1335U processor with 16 GB of RAM, the Event-Driven Monte Carlo (EDMC) approach completed the simulation in a mere 3 s for parameters M = 1000, α = 0.8, and T _max = 50, as measured by the MATLAB tic/toc routine. In contrast, the Time-Driven Monte Carlo (TDMC) method exhibited an increase in computational cost as the temporal resolution was refined; execution times surged from 119 s at ∆t = 0.5 to 1199 s at ∆t = 0.05. These results underscore the inherent advantage of the event-driven scheme, which achieves superior precision by leaping directly between state transitions, whereas the time-driven approach incurs heavy overhead by processing redundant steps during periods of system inactivity.

4 CONCLUSION

In this study, we present a detailed numerical and conceptual analysis of a discrete-time Monte Carlo framework for fractional chemical kinetics. We focus particularly on comparing it with the generalised Gillespie algorithm and the exact solution in the hydrodynamic limit (N → ∞).

Introducing fractional derivatives into the governing chemical process introduces a memory effect. In a stochastic process, this implies that the probability of a transition occurring depends on the entire past history or the ’waiting time’ since the last event. In a fixed time-step scheme, the process becomes Markovian (memoryless) if the transition probability is treated as constant within each interval. In order to capture memory effects in discrete time, an instantaneous fractional transition probability p(t) must be employed that depends explicitly on the elapsed time.

The numerical results demonstrate that, for moderate fractional orders (α ≳ 0.5), the ensemble-averaged Monte Carlo trajectories converge reliably to the exact Mittag-Leffler solutions of the corresponding time-fractional kinetic equation.

However, a systematic error analysis revealed that the performance of discrete-time schemes deteriorates significantly for small fractional orders (α ≪ 1). In this strongly non-Markovian regime, long waiting times and the absence of a finite mean waiting time exacerbate discretisation errors and restrict the effectiveness of ensemble averaging. The time-driven Monte Carlo method is inefficient for systems in which the acceptance probability of a move is very low, as substantial time is spent processing rejected attempts. By contrast, event-driven algorithms (SSA/Gillespie) select only events that will occur, rather than attempting moves that may fail. These findings suggest that discrete-time Monte Carlo methods require further refinement.

This framework provides a robust basis for expanding discrete-time fractional Monte Carlo methods to encompass more intricate reaction networks, reversible processes and systems coupled to anomalous transport. The proposed approach provides a unified and comprehensive view of non-Markovian chemical dynamics.

Data availability

Do not apply.

REFERENCES

• ¹ C.N. Angstmann, S.J.M. Burney, B.I. Henry, D.S. Han, B.A. Jacobs & Z. Xu. A Stochastic Simulation Method for Fractional Order Compartment Models (2024). URL https://arxiv.org/abs/2312.05268
  » https://arxiv.org/abs/2312.05268
• ² C.N. Angstmann, A.M. Erickson, B.I. Henry, A.V. McGann, J.M. Murray & J.A. Nichols. Fractional order compartment models. SIAM J. Appl. Math., 77(2) (2017), 430-446.
• ³ T. Aquino & M. Dentz. Chemical Continuous Time Random Walks. Phys. Rev. Lett., 119 (2017). doi:10.1103/PhysRevLett.119.230601. Article 230601.
  » https://doi.org/10.1103/PhysRevLett.119.230601
• ⁴ D.A. Beard & H. Qian. “Cambridge texts in biomedical engineering: Chemical biophysics: Quantitative analysis of cellular systems”. Cambridge University Press, Cambridge, England (2008).
• ⁵ D.O. Cahoy, V.V. Uchaikin & W.A. Woyczynski. Parameter estimation for fractional Poisson processes. J. Stat. Plan. Inference, 140(11) (2010), 3106-3120.
• ⁶ M. Caputo. Linear Models of Dissipation whose Q is almost Frequency Independent-II. Geophys. J. Int., 13(5) (1967), 529-539.
• ⁷ K. Diethelm. “The analysis of fractional differential equations”. Lecture Notes in Mathematics. Springer, Berlin, Germany, 2010 ed. (2010).
• ⁸ P. Ehrenfest & T. Ehrenfest. “The conceptual foundations of the statistical approach in mechanics”. Dover Books on Physics & Chemistry. Dover Publications (2014).
• ⁹ G. Ertl. Reactions at surfaces: from atoms to complexity (Nobel Lecture). Angew. Chem. Int. Ed Engl., 47(19) (2008), 3524-3535.
• ¹⁰ R. Garrappa. Numerical evaluation of two and three parameter Mittag-Leffler functions. SIAM J. Numer. Anal., 53(3) (2015), 1350-1369.
• ¹¹ D.T. Gillespie. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22(4) (1976), 403-434.
• ¹² D.T. Gillespie. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem., 81(25) (1977), 2340-2361.
• ¹³ D.T. Gillespie. Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem., 58(1) (2007), 35-55.
• ¹⁴ R. Gorenflo, A.A. Kilbas, F. Mainardi & S.V. Rogosin. “Mittag-Leffler functions, related topics and applications”. Springer Monographs in Mathematics. Springer, New York, NY, 2014 ed. (2014).
• ¹⁵ I. Goychuk. “Viscoelastic Subdiffusion: Generalized Langevin Equation Approach”. John Wiley & Sons, Ltd (2012), p. 187-253. doi:10.1002/9781118197714.ch5.
  » https://doi.org/10.1002/9781118197714.ch5
• ¹⁶ J.M. Hammersley. “Monte Carlo Methods”. Monographs on Statistics and Applied Probability. Springer, Dordrecht, Netherlands (2011).
• ¹⁷ M. Kac. Random walk and the theory of Brownian motion. Am. Math. Mon., 54(7) (1947), 369.
• ¹⁸ A.A. Kilbas, H.M. Srivastava & J.J. Trujillo. “Theory and applications of fractional differential equations: Volume 204”. North-Holland Mathematics Studies. Elsevier Science, London, England, 204 ed. (2006).
• ¹⁹ J. Klafter & I.M. Sokolov. “First steps in random walks”. Oxford University Press, London, England (2011).
• ²⁰ R. Kopelman. Fractal reaction kinetics. Science, 241(4873) (1988), 1620-1626.
• ²¹ T.J. Kozubowski. Computer simulation of geometric stable distributions. Journal of Computational and Applied Mathematics, 116(2) (2000), 221-229. doi:10.1016/S0377-0427(99)00318-0.
  » https://doi.org/10.1016/S0377-0427(99)00318-0
• ²² T.J. Kozubowski & S.T. Rachev. Univariate Geometric Stable Laws. J. Comput. Anal. Appl., 1(2) (1999), 177-217.
• ²³ N.H.T. Lemes, J.P.C. dos Santos & J.P. Braga. A generalized Mittag-Leffler function to describe nonexponential chemical effects. Appl. Math. Model., 40(17-18) (2016), 7971-7976.
• ²⁴ J. Liang & H. Qian. Computational cellular dynamics based on the chemical master equation: A challenge for understanding complexity. J. Comput. Sci. Technol., 25(1) (2010), 154-168.
• ²⁵ J. Liang & H. Qian. Computational cellular dynamics based on the chemical master equation: A challenge for understanding complexity. J. Comput. Sci. Technol. , 25(1) (2010), 154-168.
• ²⁶ D.A. McQuarrie. Stochastic approach to chemical kinetics. J. Appl. Probab., 4(3) (1967), 413-478.
• ²⁷ N. Metropolis & S. Ulam. The Monte Carlo method. J. Am. Stat. Assoc., 44(247) (1949), 335-341.
• ²⁸ N.C. Metropolis. The Beginning of the Monte Carlo Method. Los Alamos Science, (15) (1987), 125-130. URL https://mcnp.lanl.gov/pdf_files/Article_1987_LAS_Metropolis_125-130.pdf Special Issue on Stanislaw Ulam; also published as LA-UR-88-9067.
  » https://mcnp.lanl.gov/pdf_files/Article_1987_LAS_Metropolis_125-130.pdf
• ²⁹ R. Metzler & J. Klafter. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep., 339(1) (2000), 1-77.
• ³⁰ R. Metzler & J. Klafter. From stretched exponential to inverse power-law: fractional dynamics, Cole-Cole relaxation processes, and beyond. J. Non Cryst. Solids, 305(1-3) (2002), 81-87.
• ³¹ W.D. Moebs & E.A. Haglund. A simple Monte Carlo method for teaching chemical kinetics. J. Chem. Educ., 53(8) (1976), 506.
• ³² N.Z. Monteiro, S.R. Mazorche & R.W. dos Santos. Constructive fractional models through Mittag-Leffler functions. Computational and Applied Mathematics, 43(4) (2024), 177. doi:10.1007/s40314-024-02680-z.
  » https://doi.org/10.1007/s40314-024-02680-z
• ³³ E.W. Montroll & G.H. Weiss. Random walks on lattices. II. J. Math. Phys., 6(2) (1965), 167-181.
• ³⁴ F. Nascimento, B. Nascimento Junior, L. Cardoso, R. Albuquerquer & N. Oliveira-Neto. An approach to the kinetics and thermodynamics of elementary chemical reactions using a stochastic model. Quim. Nova, (2018).
• ³⁵ I. Podlubny. “Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications”, volume 198 of Mathematics in Science and Engineering. Academic Press, 1 ed. (1998). doi:10.1016/S0076-5392(98) 80001-8.
  » https://doi.org/10.1016/S0076-5392(98) 80001-8
• ³⁶ M.J. Saxton. Anomalous diffusion due to obstacles: a Monte Carlo study. Biophys. J., 66(2 Pt 1) (1994), 394-401.
• ³⁷ H. Scher & E.W. Montroll. Anomalous transit-time dispersion in amorphous solids. Phys. Rev., 12(6) (1975), 2455-2477.
• ³⁸ H. Schiessel & A. Blumen. Hierarchical analogues to fractional relaxation equations. J. Phys. A Math. Gen., 26(19) (1993), 5057-5069.
• ³⁹ I.M. Sokolov, J. Klafter & A. Blumen. Fractional kinetics. Phys. Today, 55(11) (2002), 48-54.
• ⁴⁰ D. Stoyan & G. Döge. 50 Jahre Metropolis et al. Mitteilungen der Deutschen Mathematiker-Vereinigung, 11(1) (2003), 17-22. doi:10.1515/dmvm-2003-0008.
  » https://doi.org/10.1515/dmvm-2003-0008
• ⁴¹ N.G. Van Kampen. “Stochastic processes in physics and chemistry”. North-Holland Personal Library. North-Holland, Oxford, England, 3 ed. (2007).
• ⁴² G. Varshney & T.K. Gupta. Recent advances in mathematical modeling of reaction kinetics: A critical review. J. Mountain Res., 19(1) (2024).
• ⁴³ M. Weiss, M. Elsner, F. Kartberg & T. Nilsson. Anomalous subdiffusion is a measure for cytoplasmic crowding in living cells. Biophys. J. , 87(5) (2004), 3518-3524.
• ⁴⁴ Z. Xu, C.N. Angstmann, D. Han, B.I. Henry, S.J.M. Burney & B.A. Jacobs. An exact stochastic simulation method for fractional order compartment models. SIAM J. Appl. Math. , 84(5) (2024), 2132-2151.
• ⁴⁵ A. Yadav & W. Horsthemke. Kinetic equations for reaction-subdiffusion systems: derivation and stability analysis. Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 74(6 Pt 2) (2006). Article 066118.

Funding

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) - Finance Code 001, and Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG) - APQ-02561-25.

APPENDIX

Algorithm 1
Time-Step Monte Carlo Simulation of Fractional Dynamics.

1:	Input:N0 (initial number of particles), λ (decay parameter), ∆t (time step), t max (maximum simulation time), α (fractional order, 0 < α ≤ 1)
2:	Initialization:
3:	N ← N 0
4:	t0 ← 0
5:	𝒯 ← {t i = i∆t | i = 0, 1, . . . , ⌊t max /∆t⌋}
6:	NMC(t0) ← N 0
7:	fori = 1 to |𝒯| − 1 do
8:	t ← 𝒯[i]
9:	x ← −λt α
10:	Compute the instantaneous fractional transition rate
	$p\left(t\right)\leftarrow \frac{1}{t}\left|\frac{E_{\alpha ,0}\left(x\right)}{E_{\alpha ,1}\left(x\right)}\right|$
11:	ndec ← 0
12:	forj = 1 to N do
13:	Generate r ∼ 𝒰(0, 1)
14:	ifr < p(t) ∆t then
15:	ndec ← n dec + 1
16:	end if
17:	end for
18:	N ← N − n dec
19:	NMC(t) ← N
20:	end for
21:	Theoretical reference solution:
22:	for allt ∈ 𝒯 do
23:	Nth(t) ← N 0 E α,1 (−λt α )
24:	end for
25:	Output: Time evolution N MC(t) and comparison with N th(t)

*

Corresponding author: Nelson Henrique Teixeira Lemes - E-mail: [e-mail]

Associate editor:

Lisandro Armando Parente

Publication Dates

• Publication in this collection
  26 Oct 2026
• Date of issue
  2026

History

• Received
  28 Mar 2026
• Accepted
  03 Aug 2026