Time-fractional cocaine–heroin dynamics with spatially heterogeneous diffusion: analysis and fractional Parker–Sochacki approximation

Abstract

Fractional reaction–diffusion models provide a powerful framework for describing dynamical systems in which memory effects and spatial interactions are significant. In this study, we develop and analyze a spatially heterogeneous time-fractional reaction–diffusion model for the coupled dynamics of cocaine and heroin abuse. Memory effects associated with addiction persistence, delayed behavioral responses, and relapse are incorporated via Caputo time-fractional derivatives, whereas spatial heterogeneity is represented by space-dependent diffusion coefficients, including both smoothly varying diffusion profiles and multi-patch configurations that capture regional disparities. To efficiently approximate the resulting nonlinear fractional system, we propose a fractional extension of the Parker–Sochacki method for time-fractional reaction–diffusion equations with heterogeneous diffusion. The proposed scheme provides an explicit and computationally efficient numerical algorithm that significantly reduces the substantial memory requirements typically associated with fractional-order models. The mathematical properties of the model are investigated rigorously. The existence, uniqueness, positivity, and boundedness of solutions are established using the theory of sectorial operators, and the equilibrium states are characterized. Furthermore, global stability conditions are derived using appropriate Lyapunov functionals. Numerical simulations validate the proposed approach and demonstrate the combined effects of memory and spatial heterogeneity on cocaine–heroin dynamics. In particular, lower fractional orders are shown to promote prolonged persistence of substance abuse, while heterogeneous diffusion generates pronounced nonuniform spatial distributions. The proposed framework integrates methodological advances with epidemiological insight, offering a robust tool for studying substance-use dynamics in spatially heterogeneous environments.

1 Introduction

The misuse of psychoactive substances, including cocaine and heroin, remains a major global public-health concern, driven by complex interactions among behavioral dependence, relapse, treatment engagement, and population-level social processes [1, 2]. From a theoretical standpoint, these dynamics generate nonlinear feedback mechanisms, threshold effects, and multiple long-term outcomes that are not readily accessible through empirical observation alone. Mathematical modeling provides a principled framework for analyzing such systems, enabling rigorous investigation of asymptotic behavior, stability properties, and the emergence of collective patterns from individual-level interactions. In this context, mechanistic population models formulated through dynamical systems theory have become indispensable tools for studying addiction dynamics, complementing empirical approaches by revealing structural features that govern persistence, extinction, and coexistence of substance-use behaviors [3–5].

Classical population models of substance abuse are often formulated as systems of ordinary differential equations, implicitly assuming homogeneous mixing and neglecting spatial structure. For example, White and Comiskey [6] formulated an integer-order compartmental model for heroin epidemics, dividing the population into susceptible, untreated users, and treated users, and analyzed the basic reproduction number and the stability of equilibria. Similarly, Nyabadza and Hove-Musekwa [7] developed a deterministic model of substance abuse dynamics that incorporates treatment processes and establishes threshold conditions for persistence and control. In the context of methamphetamine abuse, Nyabadza, Njagarah, and Smith [8] proposed an integer-order model that incorporated drug-supply dynamics and examined intervention strategies through equilibrium, stability, and numerical analyses. These models, while insightful, are based on classical ODE frameworks and therefore assume spatial homogeneity and instantaneous local interactions. While such models yield important insight into average population behavior, they are inherently unable to represent geographic variability in human mobility, treatment availability, or localized clustering of drug use. From a theoretical perspective, this limitation is significant, as spatial structure can fundamentally alter stability properties, persistence conditions, and long-term outcomes. Reaction–diffusion frameworks provide a natural extension by coupling local interaction dynamics with spatial dispersal, thereby enabling the study of spatial pattern formation, wave propagation, and hotspot emergence. In epidemiological and ecological contexts, spatially extended reaction–diffusion models with constant diffusion coefficients have demonstrated how spatial structure alone can influence long-term dynamics and pattern formation [3–5, 9].

In recent years, fractional calculus has gained increasing attention as a modeling paradigm for biological and social systems in which memory and hereditary effects play a central role, see Akindeinde et al. [10] and Madani et al. [11]. Time-fractional differential equations, typically formulated using Caputo derivatives, provide a natural framework for describing long-term dependence and nonlocal temporal behavior that cannot be captured by integer-order models [12–14]. In population and epidemic dynamics, fractional models have revealed qualitatively different persistence and stability properties, often associated with slower relaxation and enhanced memory effects [15, 16]. For example, Mittag–Leffler kernel-based models have been employed to investigate e-cigarette smoking dynamics in Uçar et al. [17], and in other related behavioral processes [18–20], demonstrating the importance of nonlocal memory effects in characterizing addiction persistence and intervention outcomes. Despite these advances, most spatially extended models of drug abuse remain restricted to integer-order temporal dynamics, whereas fractional models are predominantly spatially homogeneous, limiting their ability to capture geographic variability and localized persistence.

From both analytical and computational perspectives, combining fractional time derivatives with spatially heterogeneous diffusion poses substantial challenges. Although well-posedness and stability of fractional evolution equations have been investigated using semigroup and Lyapunov-based techniques [21–23], relatively little attention has been given to heterogeneous fractional reaction–diffusion systems motivated by substance-use dynamics. Moreover, standard numerical methods for time-fractional partial differential equations often incur high computational cost due to global memory requirements, making them impractical for long-time or large-scale simulations. Considerable effort has therefore been devoted to developing efficient numerical schemes for fractional differential equations, including predictor–corrector methods [24], finite-difference approaches [25], spectral methods [26], decomposition techniques [27], and structure-preserving algorithms. Nevertheless, many existing approaches remain computationally demanding due to the global memory requirements of fractional operators, motivating the development of explicit, memory-efficient alternatives for heterogeneous fractional reaction–diffusion systems. Recent contributions in this direction further underscore the need for scalable numerical methods tailored to complex biological applications. To the best of our knowledge, no existing study simultaneously develops a time-fractional, spatially heterogeneous reaction–diffusion model for cocaine–heroin dynamics and an efficient numerical strategy tailored to such systems.

Motivated by these gaps, we develop a time-fractional reaction–diffusion framework for the coupled dynamics of cocaine and heroin abuse in spatially heterogeneous environments. Memory effects associated with addiction persistence and relapse are incorporated through Caputo fractional derivatives, while spatial variability is modeled via space-dependent diffusion coefficients. In addition to smoothly varying diffusion profiles, a multi-patch formulation with piecewise-constant diffusion is introduced to represent abrupt regional contrasts within a unified modeling structure.

The present work contributes to the literature in four main ways. First, it integrates fractional-order memory effects and heterogeneous diffusion within a unified cocaine–heroin modeling framework, thereby capturing both temporal non-locality and spatial variability in substance-use dynamics. Second, the proposed model accommodates both smoothly varying and multi-patch diffusion structures, enabling the investigation of a broad range of mobility scenarios. Third, the resulting heterogeneous fractional system is analyzed rigorously: well-posedness and positivity of solutions are established using the theory of sectorial operators; equilibrium states are derived; and the global stability of the drug-free and endemic equilibria is investigated via appropriate Lyapunov functionals. Finally, we develop a fractional Parker–Sochacki method for time-fractional reaction–diffusion systems with heterogeneous diffusion. The proposed method extends classical power-series techniques to fractional partial differential equations, accommodates nonlinear reaction mechanisms and diffusion operators in divergence form, and avoids the global memory requirements typical of standard fractional schemes while remaining explicit and computationally efficient. To the best of our knowledge, this is the first application of the fractional Parker–Sochacki method to a cocaine–heroin reaction–diffusion model. Numerical simulations illustrate the analytical findings and examine how memory effects and spatial heterogeneity jointly influence cocaine–heroin dynamics, demonstrating convergence to drug-free and endemic equilibria and revealing enhanced persistence at lower fractional orders.

Taken together, the proposed framework advances existing substance-abuse models by unifying fractional memory effects, spatially heterogeneous diffusion, rigorous qualitative analysis, and a computationally efficient fractional Parker–Sochacki scheme to a single cocaine–heroin modeling framework.

The remainder of the manuscript is organized as follows. Section 2 summarizes the fractional calculus preliminaries required for the subsequent analysis. The time'90fractional reaction–diffusion model and its underlying assumptions are formulated in Section 3. Well' posedness of the resulting heterogeneous fractional system is established in Section 4, followed by the derivation of equilibrium states and a rigorous stability analysis in Section 5. The proposed fractional Parker–Sochacki numerical method is developed in Section 6. Section 7 presents numerical simulations that illustrate the theoretical results and examine the effects of memory and spatial heterogeneity. Finally, concluding remarks and perspectives on future work are presented in Section 8.

2 Fractional calculus preliminaries

This section compiles the specific concepts from fractional calculus required for the subsequent analytical and numerical developments. The presentation is limited to those identities and algebraic properties that play a direct role in formulating the heterogeneous fractional reaction–diffusion model and in constructing the fractional Parker–Sochacki method. Throughout, 0 < α ≤ 1 denotes the fractional order of differentiation in the Caputo sense, and Γ(·) denotes Euler's Gamma function.

2.1 Caputo fractional derivative

Let u:[0, T] → ℝ be sufficiently smooth. The Caputo fractional derivative of order α∈(0, 1] is defined by

The Caputo formulation is adopted since it allows the use of classical initial conditions and is well-suited for modeling systems with temporal memory. For α = 1, the operator reduces to the standard first-order derivative.

A key identity used in the derivation of the fractional Parker–Sochacki recurrence relations is the Caputo derivative of fractional power functions. For k≥1,

This identity underpins the explicit construction of fractional power-series solutions and is repeatedly invoked in both the analytical and numerical sections.

2.2 Nonlinear convolution structure

The reaction terms in the cocaine–heroin model include polynomial nonlinearities, such as bilinear products of state variables. Within the fractional power-series framework, products of the form

are represented through the discrete convolution of the corresponding series coefficients,

This convolution structure enables nonlinear terms to be handled systematically within the Parker–Sochacki framework and ensures closure of the fractional series expansions. The stability of such nonlinear convolutions in fractional power-series representations has been analyzed in Jornet [28].

3 Model formulation

Let U ⊂ ℝ² be a bounded spatial domain with smooth boundary ∂U, and let b>0 denote a fixed final time. Motivated by the spatio-temporal cocaine–heroin SCHR model introduced in Zinihi et al. [5] and its subsequent extension, we formulate a spatially heterogeneous cocaine–heroin abuse model that incorporates both space-dependent diffusion and a Caputo fractional time derivative of order 0 < α ≤ 1. Fractional-order dynamics allows the model to account for memory effects and delayed behavioral responses commonly observed in substance abuse and treatment processes.

The total population is partitioned into six interacting compartments,

representing, respectively, densities of individuals susceptible to drug initiation, active cocaine users, individuals undergoing cocaine treatment, active heroin users, individuals undergoing heroin treatment, and recovered individuals who have ceased drug use. Each compartment is defined by the spatial location z∈U and the time t∈[0, b]. Transitions between compartments are governed by drug initiation, progression to more severe use, treatment and recovery, and natural mortality, reflecting the complex interplay between cocaine and heroin abuse, as shown in the transfer diagram in Figure 1.

Figure 1

The transition structure captures key behavioral features of substance-use dynamics. Susceptible individuals enter the active cocaine-user class via the bilinear incidence term βSC, which models initiation driven by social exposure and interaction with current users. Active cocaine users may enter treatment at rate μ₁, progress to heroin use at rate σ, recover naturally at rate γ₁, or leave the population due to natural mortality. Individuals undergoing treatment may relapse to active use at rates μ₂ and κ₂ for cocaine and heroin, respectively, reflecting the recurrent nature of substance dependence. The heroin-related compartments follow an analogous progression, with treatment uptake occurring at a rate κ₁ and recovery at rates γ₂ and γ₄. Recovered individuals are assumed to remain permanently abstinent during the study period.

The Caputo fractional derivative is used to capture the long-term memory and hereditary effects associated with substance use. In particular, an individual's current propensity to initiate, continue, or relapse into drug use is influenced not only by present conditions but also by accumulated past experiences, treatment history, and behavioral dependence. The fractional order 0 < α ≤ 1 quantifies the strength of memory effects, with smaller values of α corresponding to stronger memory and slower behavioral adaptation.

To capture heterogeneity in spatial mobility, individuals in each compartment are allowed to diffuse across the domain with distinct space-dependent diffusion coefficients d_i = d_i(z), i∈{S, C, U_c, H, U_h, R}, where and satisfy uniform bounds 0 < d_min ≤ d_i(z) ≤ d_max. Spatial diffusion is assumed to occur in the absence of population flux across the boundary, consistent with a closed geographical region.

The resulting dynamics are governed by the following Caputo time–fractional reaction–diffusion system:

The model is supplemented with homogeneous Neumann boundary conditions,

which impose zero flux across the boundary ∂U, where ν denotes the outward unit normal vector. Initial conditions are prescribed by

where all initial population distributions are assumed to be nonnegative and sufficiently regular. The model parameters are described in Table 1.

4 Well-posedness of the fractional cocaine–heroin model

We now establish well–posedness for the fractional, space–heterogeneous reaction–diffusion system (Equations 4–6). Our strategy follows the theory of fractional evolution equations generated by sectorial operators, together with Lipschitz continuity of the nonlinearities.

4.1 Functional setting

Let W(U): =(L²(U))⁶ and write We define A:D(A)⊂W(U) → W(U) componentwise by

with domain

We assume and d_i(z) ≥ d_min > 0. Under these conditions, each scalar operator −∇·(d_i∇·) with Neumann boundary conditions is self–adjoint, positive, and sectorial on L²(U) [14]. Hence, the block–diagonal operator A is also sectorial on W(U).

Define X: =(L²(U)∩L^∞(U))⁶, and ξ:X→W(U) by

Because ξ consists only of linear and quadratic componentwise polynomials, it defines a locally Lipschitz mapping ξ:X→W(U) (see Lemma 4.1).

With these definitions, the fractional cocaine–heroin system (Equations 4–6) can be written in the compact abstract form

Let X: =(L²(U)∩L^∞(U))⁶ be endowed with the norm

In what follows, we still denote by W(U) = (L²(U))⁶ the natural energy space, and ξ:X→W(U) by the reaction terms in Equation 8.

Lemma 4.1. The nonlinear operator ξ:X→W(U) is locally Lipschitz continuous. More precisely, for every K>0 there exists L_K>0 such that

for all ϑ¹, ϑ²∈X with , k = 1, 2.

Proof. Fix K>0 and set

If ϑ∈B_K, then by definition for all i = 1, …, 6. Let us denote , k = 1, 2 and write

Then Hölder's inequality yields

Hence, using also the linear terms in ξ(ϑ), we obtain

for some constant c₁(K) depending on K and the parameters β, η₁. A similar estimate holds for the second component, with an extra linear dependence on . The remaining components ξ₃, …, ξ₆ are linear combinations of C, U_c, H, U_h, R and thus satisfy estimates of the form

with constants c_i>0 depending only on the model parameters.

Summing the six componentwise estimates and using the definition of the W(U)–norm, we obtain

for some L_K>0 depending on K and the parameters but not on . This is the desired local Lipschitz estimate on B_K, hence ξ is locally Lipschitz on X.

4.2 Mild formulation

We consider the abstract fractional Cauchy problem (Equation 8)

where A is a sectorial diffusion operator on X and ξ:X→W(U) is locally Lipschitz (see Lemma 4.1). Following the standard theory of fractional evolution equations, the Cauchy problem is equivalent to the Volterra–Mittag–Leffler integral form

where E_α and E_{α, α} denote the Mittag–Leffler operator families generated by A.

Equation 9 is obtained by applying the Laplace transform to the Caputo derivative and using the functional calculus for sectorial operators; see [12–14, 21]. We therefore adopt (Equation 9) as the definition of a mild solution to the cocaine–heroin system. A function ϑ∈C([0, T];X) is a mild solution on [0, T] if and only if it satisfies (Equation 9) for all t∈[0, T].

4.2.1 Local existence and uniqueness of mild solutions

Since A is sectorial on W(U), the corresponding Mittag–Leffler families and are bounded on W(U); more precisely, there exist constants M₁, M₂>0 such that

for all t > 0, see Podlubny [12], Kilbas et al. [13], Gal and Warma [14], and de Andrade et al. [21]. We also note that and preserve X and are bounded on X with the same type of estimates.

Theorem 4.2 (Local existence and uniqueness of mild solutions). Assume that A is sectorial on W(U) and generates the operator families and satisfying (Equation 10), and that ξ:X→W(U) is locally Lipschitz. Then, for every ϑ₀∈X there exists T^*>0 and a unique mild solution

of the abstract problem

that is, ϑ satisfies (Equation 9) for all t∈[0, T^*]. Moreover, the solution depends continuously on the initial datum ϑ₀ in the X-norm.

Proof. By Lemma 4.1, ξ:X→W(U) is locally Lipschitz. Let ϑ₀∈X be fixed. Choose K>0 such that ||ϑ₀|| _X ≤ K/2 and consider the closed ball

where T>0 will be chosen later. Define the operator Φ:B_K→C([0, T];X) by

Step 1: Φ mapsB_Kinto itself. By Equation 10 and boundedness of E_α and E_{α, α} on X,

For ϑ ∈ B_K, Lemma 4.1 implies that ξ is Lipschitz on the set {ϑ (t) : t ∈ [0, T], ϑ ∈B_K}, so there exist constants L_K, C_K > 0 such that

for all t ∈ [0, T] and . Using Equation 10 and Fubini's theorem, we obtain

Hence,

Choosing T > 0 such that

we obtain Φ(B_K)⊂B_K.

Step 2: Φ is a contraction onB_KforT > 0 small enough. Let . Then

so by Equation 10 and the Lipschitz estimate for ξ,

Taking the supremum over t∈[0, T] gives

Choosing T > 0 so small that , the operator Φ is a contraction on B_K.

By the Banach fixed-point theorem, Φ has a unique fixed point ϑ∈B_K, which satisfies (Equation 9). This proves the existence and uniqueness of a mild solution in C([0, T];X), with T: = T^* > 0 chosen as above. Continuous dependence on ϑ₀ follows from the same contraction argument applied to the difference of the two solutions with different initial data.

Remark 4.3 (Strong solutions under additional regularity). The analysis above establishes the existence of mild solutions for initial data . If, in addition, the initial state satisfies

then the mild solution becomes a strong solution for all t > 0, in the sense that

and Equation 8 holds in W(U) pointwise in time. Thus, higher regularity of the initial distribution leads to enhanced regularity of the evolving cocaine–heroin population profile.

4.3 Positivity and biological invariance

Proposition 4.4 (Positivity of solutions). Assume that the model parameters Λ, β, η_i, μ_i, σ, γ_i, κ_i≥0 and that the initial condition satisfies

Then the unique mild solution ϑ∈C([0, T^*);X) given by Theorem 4.2 satisfies

Proof. The Neumann diffusion operator A generates a positive Mittag–Leffler family , in the sense that whenever ψ≥0 a.e. in U; see, for example, de Andrade et al. [21] and Gal and Warma [14].

On the other hand, the reaction operator ξ is quasi-positive: if satisfies ϑ_j≥0 for all j and ϑ_i = 0 for some i, then a direct inspection of the components shows that ξ_i(ϑ)≥0. For instance, if S = 0 then ξ₁(ϑ) = Λ≥0; if C = 0 then ξ₂(ϑ) = μ₂U_c≥0; and similarly for the other components.

Using the mild formulation (Equation 9) and these two facts, one obtains that the positive cone in X is invariant under the solution operator: starting from nonnegative initial data, the term remains nonnegative for all t≥0, and the integrand in Equation 9 remains nonnegative whenever ϑ(·, s)≥0. A standard approximation argument (e.g., Picard iteration or time-stepping constructing monotone sequences of approximate solutions in C([0, T];X)) ensures that positivity is preserved in the limit mild solution; see, for example, de Andrade et al. [21], Gal and Warma [14], and Podlubny [12].

5 Equilibria and stability analysis

In this section, we characterize the steady states of the fractional space–heterogeneous SCHR system (Equation 4) and derive the basic reproduction number R₀ using the next–generation approach. Because the Caputo derivative and the space–dependent diffusion vanish at equilibrium, the equilibrium points coincide with those of the corresponding spatially homogeneous ODE system studied in Zinihi et al. [5].

5.1 Equilibria and the basic reproduction number

An equilibrium satisfies the algebraic system obtained by setting all time derivatives and spatial derivatives to zero in Equation 4:

5.1.1 The drug-free equilibrium

The drug-free equilibrium corresponds to the absence of cocaine and heroin consumption, i.e. Substituting into Equation 11 gives yielding the drug-free equilibrium

5.1.2 The basic reproduction number

The infected subsystem corresponding to the drug-user compartments is given by

Let denote the infected-state vector. Following the next-generation matrix method, the infected subsystem is expressed as

where ℱ(X) contains all new drug-use generation terms, while 𝒱(X) comprises all remaining transition terms. The corresponding next-generation matrices are obtained from the Jacobians

evaluated at the drug-free equilibrium

Accordingly, the matrices of new infections and transition terms are given by

and

where

Therefore, the basic reproduction number is defined as the spectral radius of the next-generation matrix FV⁻¹, namely which yields

5.1.3 Endemic equilibrium

For the endemic equilibrium, observe that if C^* > 0 or H^* > 0, the system (Equation 11) admits a non-trivial equilibrium corresponding to persistent cocaine/heroin use. To determine the endemic equilibrium , we continue from Equation 11 and solve the resulting algebraic system sequentially.

From the third equation, we obtain

Substituting this expression into the second equation yields

Hence,

Using the expression for the basic reproduction number

it follows that

Substituting S^* into the first equation gives

and therefore

Furthermore, from the fifth equation,

Substituting this into the fourth equation yields

Consequently,

and from the sixth equation,

5.2 Stability analysis

We analyze the stability of equilibria using a fractional Lyapunov framework. Consider the autonomous Caputo system

Since the Caputo derivative has no classical chain rule, the stability arguments will rely on a derivative inequality rather than regular chain rule identities.

Lemma 5.1 (Fractional Lyapunov Inequality). Let x(t) solve , and let V:ℝ^m → ℝ be convex and C¹ with V(0) = 0. Then

Under this inequality, fractional LaSalle invariance follows: ifalong trajectories and the sublevel set {V ≤ V(x₀)} is compact and positively invariant, then V(x(t)) is nonincreasing and every solution converges to the largest invariant set contained in. When this set consists of a single equilibrium, it is globally asymptotically stable.

The above result was proved for a general convex V in Tuan and Trinh [23], and for the quadratic Lyapunov case earlier in Aguila-Camacho et al. [22].

5.2.1 Global stability of the drug–free equilibrium

Recall the drug-free equilibrium We choose the Lyapunov function

where the positive parameters α₁, α₂ and α₃ are to be determined. Further, we define the Lyapunov functional Apply Lemma 5.1 to obtain

Because the boundary term vanishes under homogeneous Neumann conditions, and , the diffusion operator produces a strictly dissipative quadratic form, ensuring

For the reaction contribution, we compute

where

Choose

so that the coefficients of U_c, H, and U_h vanish identically. Consequently,

Using

we obtain

where

Hence,

Moreover,

if and only if

Substituting these conditions into the equilibrium Equation 11 yields

Therefore, by the fractional LaSalle invariance principle, every solution converges to the drug-free equilibrium . We, therefore obtain the following result.

Theorem 5.2. If ℛ₀ < 1, then the drug-free equilibrium

is globally asymptotically stable for all 0 < α ≤ 1.

5.2.2 Global stability of the endemic equilibrium

Let

denote the endemic equilibrium defined in Section 5.1.3. Assume that ℛ₀ > 1, so that the endemic equilibrium exists and is positive.

Consider the Volterra-type Lyapunov function

where

and ω₁, ω₂, ω₃ are positive parameters to be chosen appropriately.

Define the Lyapunov functional

Since

it follows that

with equality if and only if

Applying Lemma 5.1, we obtain

We compute

and

Using homogeneous Neumann boundary conditions together with

integration by parts yields

Next, we simplify the reaction terms by recalling the endemic equilibrium identities

where

For v = S,

Using the endemic equilibrium identity we obtain

Hence, the susceptible reaction term becomes

Similarly, the remaining ξ_v(ϑ) are rewritten using the corresponding equilibrium identities.

Altogether, the reaction contribution becomes

We decompose I_react into four parts and analyze each group of terms separately.

[A.] Dissipative terms:Recall the identity: for all v > 0,

Consequently,

and

[B.] Tranfer terms:

The transfer terms combine pairwise using the arithmetic–geometric mean inequality: for X > 0, Indeed, with we obtain

by taking

Similarly, with we obtain

[C.] Incidence terms:

Next, we consider the incidence terms

Using

we obtain

Similarly,

and therefore

Combining the two incidence contributions yields

Clearly, the term

Furthermore, observe from the dissipative term for C in Equation 17:

which can be paired with the last expression in Equation 21 as

Using the equilibrium condition

it follows that

and thus

Recall Young's inequality

The mixed term is controlled via Young's inequality:

In a similar manner, the other mixed term satisfies

and is absorbed into the dissipative S and C terms as explained in the sequel.

[D.] Other term:

The remaining mixed term involving C and H in Equation 15 is

by Young's inequality.

Since the system admits a positively invariant bounded region established earlier, there exists C_max > 0 such that 0 ≤ C(t, z) ≤ C_max. Moreover, assuming uniform persistence of the endemic state, there exists H_min > 0 such that H(t, z)≥H_min for (t, z)∈[0, ∞) × U. Thus,

Recall that the dissipative H-term satisfies

Thus, choosing yields

Furthermore, it holds

Hence,

Therefore, if

the positive contribution

in Equation 26 is absorbed by the negative dissipative C-term. A similar argument shows that the terms in Equation 22 will be absorbed into the dissipative S and C terms.

Consequently, the integrands in Equation 15 are negative, and we conclude that

Also,

so 𝒢(t) is nonincreasing for all t≥0. Furthermore, the derivative vanishes only if S = S^*, C = C^*, , H = H^*, , R = R^* everywhere in U. Hence, the only invariant state satisfying is . By the fractional LaSalle's principle, every solution converges to . We obtain the following result.

Theorem 5.3. If ℛ₀ > 1, then the endemic equilibriumis globally asymptotically stable for all 0 < α ≤ 1.

6 Novel numerical method for the fractional heterogeneous PDE model

The goal of this section is to develop and implement a fractional Parker–Sochacki method for the numerical approximation of the fractional SCHR reaction–diffusion system (Equation 4). The proposed approach combines a method-of-lines discretization in space with a fractional Parker–Sochacki time integrator.

The development of efficient numerical methods for time-fractional partial differential equations remains challenging because the nonlocal nature of fractional derivatives typically requires storing and repeatedly evaluating the entire solution history. Consequently, many standard finite-difference [25] and finite-element approaches [29] incur substantial computational and memory costs, particularly in long-time simulations and heterogeneous spatial settings.

To address these challenges, the proposed scheme extends the classical Parker–Sochacki method [30–33] to fractional reaction–diffusion dynamics. By exploiting local power-series representations following spatial semi-discretization, the method avoids the global memory requirements typical of many conventional fractional solvers while remaining fully explicit and computationally efficient. Furthermore, the proposed framework naturally accommodates nonlinear reaction terms, space-dependent diffusion coefficients, and interface conditions expressed in divergence form, making it particularly suitable for the heterogeneous fractional model considered in this study.

The accuracy and reliability of the proposed scheme will be evaluated by comparing it with an established numerical method to validate its applicability to heterogeneous time-fractional reaction–diffusion systems.

6.1 Spatial discretization

The spatial diffusion operator ∇·(d(z)∇u(z, t)) is discretized in conservative divergence form to preserve positivity, mass balance and stability when the diffusion coefficient varies in space. Let be a uniform partition of (0, L) with Δz = L/(M−1), and denote by u_i(t)≈u(z_i, t) the grid values. Since the term d(z)∂_zu represents flux, we evaluate it at cell interfaces by defining midpoint diffusion values

so that flux continuity across cell boundaries is preserved even when d(z) is heterogeneous.

We then define the discrete diffusion operator 𝒟:ℝ^M → ℝ^M component-wise by

which is the standard second-order finite-difference approximation of ∂_z(d∂_zu).

Zero-flux boundary conditions ∂_νu = 0 are enforced by introducing ghost points u₀ = u₁, u_M+1 = u_M, which give

ensuring no artificial flux enters or leaves the domain. Thus 𝒟 is the discrete diffusion operator associated with the heterogeneous Laplace operator ∂_z(d(z)∂_z·).

Applying this discretization to every compartment of Equation 4 yields a coupled system of 6M fractional ODEs

where , 𝒟 is block–diagonal and F contains all nonlinear reaction terms.

6.2 Fractional time integration

Time integration of Equation 28 proceeds via a fractional Parker–Sochacki expansion. On a step t∈[t_j, t_j+Δt] each nodal value is approximated by a fractional Taylor–Mittag–Leffler series

Recall that the Caputo fractional derivative satisfies

Let G_{k, i} be the coefficient of in the expansion of the right-hand side of Equation 28. Matching coefficients of Equation 28 gives the fractional PSM recurrence

where each coefficient consists of diffusion and nonlinear contributions:

Thus, advancing one time step requires computing G_{0, i}, …, G_{N−1, i} and then evaluating (Equation 29) at t_j+1. In essence, after computing U_{k, i} for k = 0, …, N, the solution is updated via

This produces a fully explicit time step: no history convolution and no nonlocal kernel storage are required. The only memory cost is storing (N+1) coefficient vectors per state variable, making the method substantially lighter than the Grünwald–Letnikov or predictor–corrector schemes. Gamma ratios Γ(αk+1)/Γ(α(k+1)+1) are pre-computed, and spatial nodes update independently within each coefficient loop, enabling efficient CPU parallelization or GPU acceleration.

7 Numerical simulations and discussion of results

In this section, we present numerical simulations to explore the model's dynamics across different epidemiological scenarios. The governing fractional–order system is solved using the novel numerical scheme developed in Section 6, which is designed to handle the nonlocal memory effects associated with the Caputo fractional derivative. All computations are implemented in Matlab, ensuring numerical stability and efficiency over long time intervals.

The primary objective of the simulations is to investigate how the fractional order α∈(0, 1] influences the transient dynamics of the cocaine–heroin abuse system and the rate at which solutions approach equilibrium. In particular, the proposed numerical method is first validated in the classical integer-order case (α = 1) by direct comparison with the standard Matlab solver ode15s. The implementation is then extended to the fractional setting to examine dynamics in both the drug-free and endemic regimes and to quantify how memory effects alter the qualitative and quantitative behavior relative to the integer-order model.

The parameter values used in this study are summarized in Tables 2, 3. Due to the limited availability of quantitative data describing transitions between cocaine use, heroin use, treatment uptake, and relapse, the adopted parameter values are taken from the substance-abuse model developed in Zinihi et al. [5], where they were selected for biological plausibility and consistency with the available literature. These values are used to illustrate the qualitative behavior of the proposed model and to facilitate numerical investigations. The same parameter set is used to facilitate comparison with the corresponding integer-order model while isolating the effects of memory and spatial heterogeneity (see Section 7.1 below). We emphasize that the analytical results presented in this work are independent of the specific parameter values employed in the simulations. We run the simulation for z∈[0, 2] and t∈[0, 500].

7.1 Method validation

To validate the proposed numerical scheme, we first consider the classical integer-order model (α = 1) with constant diffusion coefficients d_S(z) = d_C(z) = d_H(z) = d_Uc(z) = d_Uc(z) = d_R(z) = 0.1. In this setting, we benchmark the numerical solutions obtained using the proposed method against those from the standard Matlab solver ode15s. The results in Figure 2 show excellent agreement between the two approaches across all state variables, confirming the accuracy and consistency of the developed scheme in the absence of fractional memory effects. This validation provides a reliable baseline for subsequent fractional-order simulations.

Figure 2

In what follows, we consider two representative choices of the diffusion coefficients: smoothly varying spatial diffusion and piecewise-constant (patch-based) diffusion.

7.2 Smooth heterogeneous diffusion

The heterogeneous diffusion coefficients are chosen as where L = 2. Initial population of the compartments is chosen as S₀ = 30, C₀ = 10, H₀ = 5,_Uc0 = 3,_Uh0 = 3, R₀ = 0.

7.2.1 Drug-free regime

We consider the drug-free scenario, corresponding to parameter values for which the basic reproduction number satisfies ℛ₀ = 0.694 < 1. To provide a comprehensive description of the system behavior, we report in Figures 3–5 the spatio-temporal surface plots of all state variables, together with their corresponding spatial mean profiles. This dual representation allows us to capture both the full space–time evolution of the solution and its averaged temporal dynamics.

Figure 3

Figure 4

The surface plots show that the solution trajectories remain nonnegative and evolve smoothly toward the drug-free equilibrium. Figure 6 shows that irrespective of initial conditions, the populations tend toward the drug-free equilibrium, consistent with the analytical stability results in Section 5.2. The system approaches equilibrium relatively rapidly, both in the pointwise spatial dynamics and in the mean behavior. Figure 7 shows that as the fractional order decreases, the convergence becomes progressively slower, with the mean trajectories exhibiting extended transients before stabilization.

Figure 5

Figure 6

These results suggest that the fractional order plays a biologically meaningful role in the model dynamics by encoding the long-term influence of past drug use and treatment histories. Smaller values of α correspond to stronger memory effects, which manifest as prolonged persistence of drug-use behaviors and delayed population-level recovery, even in parameter regimes where elimination is theoretically guaranteed (see Figure 7). This behavior reflects the well-documented inertia in substance-use dynamics in which behavioral change and treatment outcomes depend not only on current conditions but also on accumulated past exposure.

Figure 7

7.2.2 Endemic regime

We next examine the endemic case, for which ℛ₀ = 1.604 > 1 and the system admits a unique endemic equilibrium. Figures 8–10 illustrate the corresponding solution profiles for various fractional orders. As shown in Figure 11, the fractional order plays a crucial role in shaping the system dynamics. While the equilibrium levels remain largely unchanged across different values of α, the pathways by which the system approaches the endemic equilibrium vary significantly. Lower fractional orders lead to slower convergence and smoother trajectories, thereby delaying the attainment of endemic levels. Furthermore, Figure 12 confirms the asymptotic stability of the endemic equilibrium as proved in Theorem 5.3.

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

7.3 Patch-based diffusion

In the preceding sections, spatial periodic heterogeneity was introduced, allowing mobility to vary smoothly across the domain. While this is mathematically convenient and well-suited to environments with regular spatial repetition, the mobility patterns associated with illicit drug use, treatment access, and social mixing are rarely periodic. Instead, they are often shaped by sharp socioeconomic or infrastructural contrasts between neighboring regions, such as urban versus peri-urban areas, zones with unequal access to treatment facilities, or regions separated by administrative boundaries. It should be emphasized that the purpose of the multi-patch formulation is not to investigate additional fractional-order effects but to examine a distinct form of spatial heterogeneity arising from abrupt regional differences in mobility and treatment accessibility.

Motivated by these considerations, we introduce a multi-patch diffusion framework as a parsimonious and ecologically standard representation of spatial heterogeneity. In this setting, the spatial domain is decomposed into subregions (patches), each characterized by homogeneous mobility properties, while abrupt changes in diffusion occur across patch interfaces. Such models are widely used in reaction–diffusion systems to capture discontinuous spatial structure and regional differences in movement intensity, and they rely on classical transmission conditions enforcing continuity of population density and diffusive flux across interfaces [3, 9].

We restrict attention to the simplest nontrivial configuration, namely, a two-patch, piecewise-constant diffusion model. Let the one-dimensional spatial domain be U = (0, L), 0 < ξ < L, partitioned into two sub-domains U₁ = (0, ξ), U₂ = (ξ, L), with interface Γ = {ξ}. For each population compartment X∈{S, C, U_c, H, U_h, R}, the diffusion coefficient is assumed to be constant within each patch, but may differ between patches:

We assume that U₁ denotes a hotspot region for drug users, and that z∈(0, ξ) exhibits lower diffusion. On the other hand, z∈(ξ, L) in U₂ exhibits higher diffusion. All reaction parameters are taken to be spatially homogeneous, so that The heterogeneity enters the system exclusively through differential mobility.

Let X_i(t, z) denote the restriction of the state variable X to patch U_i (i = 1, 2). The resulting two-patch fractional reaction–diffusion model is given by

Consistent with the assumption adopted in earlier sections, we impose zero-flux (Neumann) boundary conditions at the outer boundaries:

for all compartments X. At the interface z = ξ, classical transmission conditions are enforced to guarantee physical consistency of the model. Specifically, we require

• continuity of the state variables, ensuring no artificial population jump across the interface, and

• continuity of diffusive flux, ensuring conservation of mass.

That is,

for each X∈{S, C, U_c, H, U_h, R}.

Initial data are also inherited from the baseline spatial model (Equation 4) and are assumed to be consistent across patches:

The heterogeneous diffusion coefficients are chosen as follows:

and d_Uc(z) = d_Uc(z) = d_R(z) = 0.1.

We solved the two-patch configuration with ξ = 0.9 using the space-dependent diffusion coefficients described above, applying the fractional numerical method introduced in Section 7. The resulting spatio–temporal dynamics are shown in Figures 13–16. Across all fractional orders considered, the solutions converge to the corresponding drug-free and endemic equilibria, and the spatial profiles remain consistent across patch interfaces. These simulations demonstrate the robustness and computational efficiency of the proposed numerical framework for fractional reaction–diffusion systems with piecewise-heterogeneous diffusion coefficients.

Figure 13

Figure 14

Figure 15

Figure 16

8 Conclusion

We have developed and analyzed a spatially heterogeneous, time-fractional reaction–diffusion framework for modeling the coupled dynamics of cocaine and heroin abuse. Fractional temporal dynamics were used to capture memory effects associated with addiction persistence and relapse, whereas space-dependent diffusion accounted for regional variability in mobility and intervention resources, including abrupt contrasts, using a multi-patch formulation.

The heterogeneous fractional system was shown to be well-posed, and positivity of solutions was established via sectorial operator theory. Equilibrium states corresponding to drug-free and endemic regimes were derived, and their global stability was established using appropriate Lyapunov functionals, thereby providing rigorous insight into the model's long-term behavior under memory and spatial heterogeneity.

A central methodological contribution of the study is the development of a fractional Parker–Sochacki method for time-fractional reaction–diffusion systems with heterogeneous diffusion. The method extends classical power-series techniques to fractional partial differential equations, accommodates nonlinear reaction mechanisms and interface conditions in divergence form, and avoids the global memory requirements characteristic of many standard fractional schemes. Numerical simulations demonstrated both the robustness of the method and its effectiveness in capturing the combined influence of memory and spatial heterogeneity, revealing slower relaxation and enhanced persistence for lower fractional orders, as well as nonuniform spatial patterns induced by heterogeneous diffusion.

Overall, the results highlight the importance of jointly accounting for memory effects and spatial structure in mathematical models of substance-use dynamics. Beyond the specific cocaine–heroin context, the analytical and computational framework developed here applies to a broader class of fractional reaction–diffusion systems arising in biological, epidemiological, and social settings. Future studies may consider higher-dimensional domains, additional behavioral mechanisms, and interactions among fractional dynamics, spatial heterogeneity, and control strategies.

References

• 1.

  AkereleEOluponaT. Drugs of Abuse. Psych Clin North Am. (2017) 40:501–17. doi: 10.1016/j.psc.2017.05.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 2.

  Carmona AraújoACasalRJGoulãoJMartinsAP. Misuse of psychoactive medicines and its consequences in the European Union – a scoping review. J Subst Use. (2023) 29:629–40. doi: 10.1080/14659891.2023.2213325

  • CrossRef
  • Google Scholar
  • View reference in article

• 3.

  CosnerC. Reaction–diffusion equations and ecological modeling. Tutor Mathem Biosci. (2008) 4:77–115. doi: 10.1007/978-3-540-74331-6_3

  • CrossRef
  • Google Scholar
  • View reference in article

• 4.

  LamKYLiuSLouY. Reaction–diffusion models in spatial ecology and epidemiology: theory, methods and applications. arXiv preprint arXiv:2004.07978 (2021).

  • Google Scholar
  • View reference in article

• 5.

  ZinihiASidi AmmiMREhrhardtMBachirA. Dynamical analysis of a cocaine–heroin epidemiological model with spatial distributions. Adv Cont Discr Models. (2025) 2025:1–22. doi: 10.1186/s13662-025-03924-w

  • CrossRef
  • Google Scholar
  • View reference in article

• 6.

  WhiteEComiskeyC. Heroin epidemics, treatment and ODE modelling. Math Biosci. (2007) 208:312–24. doi: 10.1016/j.mbs.2006.10.008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 7.

  NyabadzaFHove-MusekwaSD. From heroin epidemics to methamphetamine epidemics: modelling substance abuse in a South African province. Math Biosci. (2010) 225:132–40. doi: 10.1016/j.mbs.2010.03.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8.

  NyabadzaFNjagarahJBHSmithRJ. Modelling the dynamics of crystal meth (Tik) Abuse in the presence of drug-supply chains in South Africa. Bull Math Biol. (2012) 75:24–48. doi: 10.1007/s11538-012-9790-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9.

  RousselMR. Reaction-diffusion equations. Technical Report, University of Lethbridge. (2005).

  • Google Scholar
  • View reference in article

• 10.

  AkindeindeSOOkyereEAdewumiAOLebeloRSFabelurinOOMooreSE. Caputo fractional-order SEIRP model for COVID-19 Pandemic. Alexandr Eng J. (2022) 61:829–45. doi: 10.1016/j.aej.2021.04.097

  • CrossRef
  • Google Scholar
  • View reference in article

• 11.

  MadaniNHammouchZJafariH. Fractal–fractional modeling of drug addiction dynamics: capturing memory'driven effects. Math Methods Appl Sci. (2025) 49:2114–28. doi: 10.1002/mma.70232

  • CrossRef
  • Google Scholar
  • View reference in article

• 12.

  PodlubnyI. Fractional Differential Equations. London: Academic Press. (1999).

  • Google Scholar
  • View reference in article

• 13.

  KilbasAASrivastavaHMTrujilloJJ. Theory and Applications of Fractional Differential Equations. New York: Elsevier. (2006).

  • Google Scholar
  • View reference in article

• 14.

  GalCGWarmaM. The semilinear parabolic problem. In: Fractional-in-Time Semilinear Parabolic Equations and Applications, (2020). p. 63–124. doi: 10.1007/978-3-030-45043-4_3

  • CrossRef
  • Google Scholar
  • View reference in article

• 15.

  LiYChenYPodlubnyI. Mittag–Leffler stability of fractional order nonlinear dynamic systems. Automatica. (2009) 45:1965–9. doi: 10.1016/j.automatica.2009.04.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 16.

  LiYChenYPodlubnyI. Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag–Leffler stability. Comput Mathem Applic. (2010) 59:1810–21. doi: 10.1016/j.camwa.2009.08.019

  • CrossRef
  • Google Scholar
  • View reference in article

• 17.

  UçarEUçarSEvirgenFÖzdemirN. Investigation of e-cigarette smoking model with Mittag-Leffler kernel. Found Comput Dec Sci. (2021) 46:97–109. doi: 10.2478/fcds-2021-0007

  • CrossRef
  • Google Scholar
  • View reference in article

• 18.

  UçarSÖzdemirNKocaİAltunE. Novel analysis of the fractional glucose–insulin regulatory system with non-singular kernel derivative. Eur Phys J Plus. (2020) 135:414. doi: 10.1140/epjp/s13360-020-00420-w

  • CrossRef
  • Google Scholar
  • View reference in article

• 19.

  HristovJ. Sigmoids based on Mittag-Leffler functions: ideas, modeling, and computational experiments. Trans Comput Model Intell Syst. (2026) 3:10073. doi: 10.65112/tcmis.10073

  • CrossRef
  • Google Scholar
  • View reference in article

• 20.

  EvirgenFUçarSÖzdemirN. Mathematical analysis and optimal control of a Caputo fractional diabetes system with parameter identification. J Comput Appl Mathem. (2026) 477:117151. doi: 10.1016/j.cam.2025.117151

  • CrossRef
  • Google Scholar
  • View reference in article

• 21.

  de AndradeBCarvalhoANCarvalho-NetoPMMarín-RubioP. Semilinear fractional differential equations: global solutions, critical nonlinearities and comparison results. Topol Methods Nonl Anal. (2015) 45:439. doi: 10.12775/TMNA.2015.022

  • CrossRef
  • Google Scholar
  • View reference in article

• 22.

  Aguila-CamachoNDuarte-MermoudMAGallegosJA. Lyapunov functions for fractional order systems. Commun Nonl Sci Numer Simul. (2014) 19:2951–7. doi: 10.1016/j.cnsns.2014.01.022

  • CrossRef
  • Google Scholar
  • View reference in article

• 23.

  TuanHTTrinhH. Stability of fractional' order nonlinear systems by Lyapunov direct method. IET Control Theory Applic. (2018) 12:2417–22. doi: 10.1049/iet-cta.2018.5233

  • CrossRef
  • Google Scholar
  • View reference in article

• 24.

  DiethelmKFordNJFreedAD. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. (2002) 29:3–22. doi: 10.1023/A:1016592219341

  • CrossRef
  • Google Scholar
  • View reference in article

• 25.

  LiCN'GboNSuF. Finite difference methods for nonlinear fractional differential equation with ψ-Caputo derivative. Physica D. (2024) 460:134103. doi: 10.1016/j.physd.2024.134103

  • CrossRef
  • Google Scholar
  • View reference in article

• 26.

  HadaadJMAllameMAal-RkhaisHAKTavassoli-KajaniM. A Müntz spectral method for solving fractional Fredholm integro–differential equations with convergence analysis. Alexandr Eng J. (2026) 134:585–95. doi: 10.1016/j.aej.2025.12.040

  • CrossRef
  • Google Scholar
  • View reference in article

• 27.

  MomaniSShawagfehN. Decomposition method for solving fractional Riccati differential equations. Appl Math Comput. (2006) 182:1083–92. doi: 10.1016/j.amc.2006.05.008

  • CrossRef
  • Google Scholar
  • View reference in article

• 28.

  JornetM. Power-series solutions of fractional-order compartmental models. Comput Appl Mathem. (2024) 43:67. doi: 10.1007/s40314-023-02579-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 29.

  MarosGIzsákF. Finite element methods for fractional-order diffusion problems with optimal convergence order. Comput Mathem Applic. (2020) 80:2105–14. doi: 10.1016/j.camwa.2020.09.006

  • CrossRef
  • Google Scholar
  • View reference in article

• 30.

  AkindeindeSOBelloKAAdewumiAO. Multistage Parker-Sochacki method for fractional ODE and PDE models: application to the Brusselator system. Boletim da Soc Paran Matem. (2025) 43:1–14.

  • Google Scholar
  • View reference in article

• 31.

  AkindeindeSAdesanyaSLebeloRSMoloiKC. A new multistage Parker-Sochacki method for solving the Troesch's problem. Int J Eng Technol. (2020) 9:592. doi: 10.14419/ijet.v9i2.13231

  • CrossRef
  • Google Scholar
  • View reference in article

• 32.

  AkindeindeSOAdewumiOALebeloRS. Approximate analytical solution of a power-law pseudoplastic fluid in a boundary layer over a porous plate: a new application of multi-stage parker-sochacki method. Int J Eng Res Africa. (2021) 55:1–14. doi: 10.4028/www.scientific.net/JERA.55.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 33.

  AkindeindeSO. A new multistage technique for approximate analytical solution of nonlinear differential equations. Heliyon. (2020) 6:e05188. doi: 10.1016/j.heliyon.2020.e05188

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article