A Deterministic and Stochastic Fractional-Order Model for Computer Virus Propagation with Caputo-Fabrizio Derivative: Analysis, Numerics, and Dynamics

1  Introduction

Rapid computer network expansion in the digital age has made systems more vulnerable to malware, especially viruses. Antivirus techniques and robust network topologies depend on understanding viral transmission mechanisms. Mathematics provides a formal framework for virus dynamics, outbreak thresholds, and intervention evaluation. Classical models, based on integer-order differential equations, focus on virus spread but cannot account for memory and hereditary effects in complex digital and biological systems [1–3].

Fractional calculus has developed into a robust extension of conventional modelling methods, enabling derivatives of non-integer order to include non-locality, memory, and hereditary characteristics [4]. Initial formulations like the Riemann–Liouville and Caputo derivatives are defined by singular kernels, limiting their use for modelling smooth memory processes in practical scenarios [5]. To address these constraints, Caputo and Fabrizio proposed a fractional derivative including a nonsingular exponential kernel [6–8].

Fractional calculus has emerged as a powerful framework for modeling complex systems with memory and hereditary properties that cannot be adequately captured by classical integer-order derivatives. A comprehensive overview of its real-world relevance is provided by Baleanu et al. [9], who demonstrated the effectiveness of fractional operators in diverse applications spanning physics, engineering, and applied sciences. Building on these foundations, Atangana [10] introduced the concept of fractal–fractional differentiation, which unifies fractal geometry with fractional calculus and enables a more faithful representation of multi-scale and irregular dynamics. This formulation has proven particularly useful in the analysis of chaotic systems. For instance, Alsulami et al. [11] investigated the controlled chaos of the Newton–Leipnik system using fractal–fractional operators, revealing enhanced control performance and richer dynamical behavior compared to classical models. Similarly, Yan et al. [12] applied Adomian decomposition techniques in conjunction with fractal–fractional operators to a Lorenz-type system, highlighting the significant impact of memory effects on chaotic trajectories and convergence properties. In parallel, Alhazmi et al. [13] employed the Caputo–Fabrizio operator to study chaotic systems with nonsingular kernels, demonstrating improved numerical stability and an accurate depiction of fading memory effects. Collectively, these studies underscore the growing importance of fractal–fractional and nonsingular fractional operators in the analysis, control, and numerical approximation of nonlinear and chaotic dynamical systems. The use of fractional-order differential equations in dynamical system modeling has a long and well-established history, particularly in the context of mechanical systems exhibiting damping and memory effects [14, 15].

More recently, comprehensive investigations of fractional glucose–insulin systems have been conducted, addressing both theoretical and numerical aspects. Saber and Alahmari [16] analyzed existence, stability, and numerical behavior using residual power series and generalized Runge–Kutta methods, offering detailed comparisons among various fractional operators. The application of fractal–fractional operators to nonlinear glucose–insulin models was further advanced by Althubyani et al. [17], who derived semi-analytical solutions via the Adomian decomposition method and demonstrated the strong influence of fractal memory on system dynamics. In parallel, Saber and Mirgani [18] developed and analyzed several numerical schemes based on Laplace residual power series methods and Atangana–Baleanu as well as Caputo derivatives, revealing rich dynamical behaviors, including stability transitions and chaos control, particularly under variable-order formulations. Saber et al. [19] employed the Jumarie–Stancu collocation series method alongside a multi-step generalized differential transform approach to construct efficient numerical solutions. Furthermore, the robustness of fractional glucose–insulin systems was rigorously examined through Hyers–Ulam stability analysis and control strategies in [20], ensuring reliability of approximate solutions under perturbations. Collectively, these studies demonstrate that fractional, fractal–fractional, and variable-order modeling frameworks, together with advanced numerical and stability analyses, provide a powerful and flexible foundation for capturing the complex dynamics of glucose–insulin regulation in diabetic systems. Further contributions encompass the creation of novel numerical approximations utilising non-local and non-singular kernels [21], and the application of fractal–fractional operators in chaotic and biological systems [22,23]. Additionally, numerical techniques tailored for Caputo–Fabrizio (CF) operators, including multi-step Adams–Bashforth (AB) methods, Laplace–Adomian decomposition methods (LADM), and Newton polynomial-based methods (NPM), have been suggested to guarantee stability and precision in addressing these systems [24,25].

Choice of Caputo–Fabrizio Operator. The selection of the Caputo–Fabrizio derivative for this cyber-epidemiological model is motivated by several key considerations specific to computer virus dynamics. First, the exponential kernel of the CF operator exp⁡[−α(t−τ)1−α] naturally represents fading memory effects, which align with the gradual decay of influence from past infection states in digital systems—such as delayed patch deployments, gradual antivirus signature updates, and persistent latent infections. This contrasts with the Mittag-Leffler kernel of the Atangana–Baleanu (AB) operator, which exhibits power-law decay that may overemphasize distant past states in cyber contexts where memory effects typically diminish exponentially. Second, the CF derivative’s non-singular kernel avoids computational singularities at t=τ, enabling more stable numerical simulations crucial for real-time cybersecurity applications. Third, the CF formulation facilitates the development of efficient numerical schemes (like our adapted Adams–Bashforth method) with proven convergence properties, whereas AB operators often require more complex numerical handling of the Mittag-Leffler function. Finally, the CF parameter α provides a clear interpretation as a memory strength index: α→1 recovers memoryless classical dynamics, while α→0 corresponds to strong hereditary effects—an intuitive mapping for cybersecurity practitioners designing memory-aware defense strategies.

The modelling of computer viruses has progressed from deterministic compartmental frameworks to stochastic and fractional-order models. Classical models classify the computer population into susceptible, infected, and recovered compartments [26,27], whereas more sophisticated systems incorporate latent stages, antivirus activation, and network feedback effects [28]. The integration of fractional derivatives facilitates a more accurate depiction of memory-dependent processes in digital contexts, enhancing predictions regarding virus survival and control efficacy. Furthermore, fractional-order models have lately been utilised in various applied sciences, including thermoelasticity and nanofluid dynamics, underscoring its transdisciplinary potential [29–32]. Mathematical modeling has long played a pivotal role in understanding the dynamics of computer virus propagation and the effectiveness of control strategies in networked systems. Classical integer-order epidemic models have been widely employed to describe virus transmission mechanisms; however, such models often fail to account for inherent delay effects, long-term memory, and hereditary characteristics observed in complex cyber environments. To address these limitations, nonstandard computational approaches incorporating delay and memory effects have been proposed, demonstrating improved realism and analytical depth in cyber-epidemiological modeling [33]. In recent years, fractional calculus has emerged as a powerful extension of classical modeling frameworks, offering enhanced capability to represent nonlocal and memory-dependent dynamics through derivatives of non-integer order. Efficient numerical methods for fractional differential equations have further strengthened their applicability across physics, engineering, and applied sciences by enabling accurate simulation of systems governed by memory effects [34]. Moreover, nonlinear and chaotic behaviors have been identified in virus-related and biological models, highlighting the necessity of advanced mathematical tools to capture complex dynamical phenomena that arise in both biological and cyber systems [35]. From a theoretical perspective, ensuring the well-posedness of fractional-order models—particularly the existence and uniqueness of solutions—is essential for validating their mathematical and practical relevance, and such foundational properties have been rigorously established for various classes of fractional and nonlocal problems [36,37].

Model Scope and Limitations. This study presents a novel fractional-order computer virus propagation model based on the Caputo–Fabrizio derivative. The primary objectives are to establish a theoretical framework for incorporating memory effects into cyber-epidemiological dynamics and to develop stable numerical methods for simulating such systems. The model is validated through comprehensive numerical simulations and comparative analyses with classical integer-order and alternative fractional formulations. However, it is important to note the following limitations: (1) The model assumes homogeneous mixing within the computer population and does not explicitly account for network topology or contact heterogeneity, which are known to significantly influence real-world infection dynamics; (2) The validation is conducted numerically without calibration against empirical virus spread datasets, which limits direct applicability to specific real-world scenarios; (3) Additional factors such as patching delays, human behavior patterns, and adaptive defense mechanisms are not incorporated in the current formulation. These simplifications enable tractable mathematical analysis and clear isolation of memory effects, but future extensions should address these aspects for enhanced practical relevance.

Advancements Beyond Existing CF-Based Cyber-Epidemiological Models. While recent works have employed Caputo–Fabrizio derivatives in virus and worm propagation modeling [38], this study introduces several key theoretical and methodological advances. First, the proposed compartmental structure {S, L, B, R} explicitly separates latently infected (L) from actively breaking-out (B) computers, providing cyber-specific granularity that more accurately represents digital infection chains—unlike the generic {S, E, I, R} frameworks commonly adopted in prior CF-based cyber models. Second, we develop a novel two-step Adams–Bashforth scheme specifically adapted for CF operators, achieving proven second-order convergence with computational efficiency comparable to integer-order schemes, thereby advancing beyond the first-order discretizations typically employed in CF implementations. Third, the complete stochastic formulation, together with the derived SCF–AB scheme, enables probabilistic risk assessment and ensemble simulations, whereas most existing CF-based cyber models remain deterministic. Fourth, the integrated stability framework—including Hyers–Ulam stability analysis and comprehensive sensitivity assessment—provides robustness guarantees and quantitatively evaluates parameter influences relevant to cybersecurity policy (e.g., the transmission rate β exhibits a sensitivity of +0.82, while the memory parameter α shows a sensitivity of −0.29). These combined advances in cyber-specific modeling, numerical methodology, and practical cybersecurity applicability clearly distinguish the proposed framework from existing CF-based epidemiological approaches.

Practical Implications and Future Directions. Despite these limitations, the theoretical insights from this fractional modeling approach can inform the design of more resilient cybersecurity systems. The memory effects captured by the Caputo–Fabrizio derivative provide a mathematical basis for understanding delayed responses in cyber-epidemics, which aligns with observed behaviors in real networks. Future research could integrate this framework with network-aware epidemic models [39,40], human behavior factors [41], and adaptive defense strategies [42] to develop more accurate predictive tools. Additionally, the stochastic extension of the model supports uncertainty quantification in dynamic threat environments, which can be leveraged by machine learning models for anomaly detection [43].

Recent advances in computer virus modeling have increasingly focused on capturing memory effects, network heterogeneity, and complex dynamical behaviors that cannot be adequately represented by classical integer-order formulations [44]. Fractional-order models, in particular, have proven effective in describing nonlocal interactions and hereditary characteristics inherent in cyber systems. Notably, the Atangana–Baleanu fractional derivative with a nonsingular Mittag–Leffler kernel has been successfully employed to investigate computer virus spread in networked environments, demonstrating how fractional memory significantly influences virus persistence and control efficiency [45]. In parallel, modified compartmental virus models incorporating network attributes and key system parameters have provided deeper insights into stability regions, outbreak thresholds, and the role of connectivity in shaping virus dynamics [46]. From a theoretical standpoint, the stability analysis of fractional-order nonlinear systems remains a central concern, and Lyapunov-based approaches combined with generalized Mittag–Leffler stability theory offer rigorous criteria for assessing asymptotic behavior and robustness of equilibria in fractional dynamical systems [47]. These developments collectively motivate the use of fractional derivatives and advanced stability tools in the present study to achieve a more realistic and mathematically sound representation of computer virus propagation.

This study presents a novel computer virus propagation model utilising the Caputo–Fabrizio derivative. The model categorises the network into four compartments: susceptible, latently infected, actively spreading, and antivirus-capable systems. The use of the CF derivative facilitates the explicit incorporation of memory and hereditary effects, thereby accurately representing the dynamics of virus spread and antivirus response. We demonstrate the existence and uniqueness of solutions through fixed-point theory, contingent upon Lipschitz continuity and boundedness conditions, and we develop a two-step Adams-Bashforth numerical approach tailored for CF operators. The scheme’s convergence and stability are thoroughly examined, and comprehensive simulations are conducted for different fractional orders to evaluate the influence of memory on virus dynamics. The expanded framework integrates stochastic influences via a stochastic CF formulation and encompasses a thorough sensitivity analysis to pinpoint critical parameters affecting viral dynamics. The findings indicate that fractional-order dynamics substantially influence outbreak intensity and timing, providing insights into effective antivirus techniques and cyber defence policy.

The key contributions of this study are:

(i)   The formulation of a novel CF-based computer virus propagation model with four network compartments;

(ii)   The theoretical proof of existence and uniqueness of solutions under standard assumptions;

(iii)   The development of a stable and convergent two-step Adams–Bashforth method adapted to CF operators;

(iv)   Comprehensive numerical investigations illustrating how memory effects influence infection dynamics and antivirus efficiency;

(v)   Extension to stochastic framework and sensitivity analysis for robust cybersecurity policy insights.

The remaining parts of the paper are structured as follows: Section 2 presents the mathematical formulation of the suggested model and its characteristics. It establishes results pertaining to existence and uniqueness. Also delineates the numerical methodology and convergence assessment. It ultimately provides numerical simulations, a comparative study with classical and alternative fractional models, and a sensitivity analysis. Section 3 expands the approach to encompass stochastic situations. Section 4 offers an extensive analysis of the results, while Section 5 ends with final observations and possible extensions.

2  Model Properties of the Caputo–Fabrizio Fractional Computer Virus Model

2.1 Model Description and Preliminaries

This section presents a refined formulation of a fractional-order computer virus propagation model using the Caputo–Fabrizio (CF) fractional derivative. The CF operator provides a non-singular exponential kernel, making it well-suited for capturing memory-dependent dynamics in cyber-epidemiological processes such as latency, delayed activation, and antivirus response.

We divide the total computer population into the following four compartments:

•   S⁡(t): Susceptible computers not yet infected,

•   L⁡(t): Latently infected computers,

•   B⁡(t): Breaking-out computers actively transmitting the virus,

•   R⁡(t): Antivirus-protected computers.

The classical (integer-order) virus transmission model is

dSdt=μ−βS(L+B)−μS,dLdt=βS(L+B)−(ϵ+μ)L,dBdt=ϵL−(γ+μ)B,dRdt=γB−(ϱ+μ)R.(1)

The Caputo fractional model of order 0<α<1 is

DtαCS(t)=μ−βS(L+B)−μS,DtαCL(t)=βS(L+B)−(ϵ+μ)L,DtαCB(t)=ϵL−(γ+μ)B,DtαCR(t)=γB−(ϱ+μ)R,(2)

where the Caputo derivative is defined as

DtαCf(t)=1Γ(1−α)∫0tf′(τ)(t−τ)αdτ.

For 0<α<1, the Caputo–Fabrizio (CF) fractional derivative is defined as [6,7]

𝒟tαCFf(t)=11−α∫0tf′(τ)exp[−α(t−τ)1−α]dτ.

Using the CF operator, the virus propagation model becomes

𝒟tαCFS(t)=μ−βS(t)(L(t)+B(t))−μS(t),𝒟tαCFL(t)=βS(t)(L(t)+B(t))−(ϵ+μ)L(t),𝒟tαCFB(t)=ϵL(t)−(γ+μ)B(t),𝒟tαCFR(t)=γB(t)−(ϱ+μ)R(t).(3)

System (3) is the Caputo-Fabrizio Fractional Computer Virus Model that will be analyzed in this work. The model parameters and their interpretations are summarized in Table 1.

The fractional order α controls the strength of memory in the system: α→1 recovers memoryless classical dynamics, whereas α→0 corresponds to strong memory influence. The parameters, particularly β, allow investigation of Virus-Free Steady States (VFSS) and Virus-Persistent Steady States (VPSS), such as β=0.5 (VFSS) and β=0.8 (VPSS). This generalized framework preserves the essential structure of epidemiological models while incorporating realistic memory effects inherent in cyber-epidemics.

2.2 Boundedness of Solutions

Before establishing Lipschitz continuity, we first demonstrate that the state variables remain bounded within a biologically feasible region.

Lemma 1 (Boundedness): Let Ω={(S,L,B,R)∈R+4:S+L+B+R≤1} denote the feasible region. For any initial condition (S(0),L(0),B(0),R(0))∈Ω and any fractional order α∈(0,1], the solution of the Caputo-Fabrizio fractional-order system (3) remains in Ω for all t≥0. Consequently, there exists a positive constant M=1 such that

0≤S(t),L(t),B(t),R(t)≤Mfor all t∈[0,T].

Proof: Consider the total computer population N(t)=S(t)+L(t)+B(t)+R(t). From system (3), we have:

𝒟tαCFN(t)=μ−μS(t)−(ϵ+μ)L(t)+ϵL(t)−(γ+μ)B(t)+γB(t)−(ϱ+μ)R(t).

Simplifying this expression yields:

𝒟tαCFN(t)=μ−μN(t)−ϱR(t)≤μ−μN(t).

This inequality holds because ϱR(t)≥0. Now consider the comparison fractional differential equation:

𝒟tαCFy(t)=μ−μy(t),y(0)=N(0).(4)

Applying the Laplace transform method for Caputo-Fabrizio derivatives [7], the solution of (4) is:

y(t)=1+(N(0)−1)exp⁡(−αμt1−α+αμ).

By the comparison principle for fractional differential equations [47], it follows that N(t)≤y(t) for all t≥0. Since N(0)≤1 (by the initial condition in Ω), we obtain:

N(t)≤1+(N(0)−1)exp⁡(−αμt1−α+αμ)≤1.

Therefore, S(t)+L(t)+B(t)+R(t)≤1 for all t≥0. Since all state variables are nonnegative by their biological interpretation, we conclude that 0≤S(t),L(t),B(t),R(t)≤1 for all t∈[0,T]. Thus, we can take M=1. ◻

Remark 1: The boundedness result holds for both the integer-order case (α=1) and the fractional-order case (0<α<1). The exponential term in the solution of (4) reduces to the standard exponential e−μt when α=1, consistent with classical epidemic models.

2.3 Equilibrium and Stability Analysis

The equilibrium states of both the integer-order system (1) and the fractional-order system (2) are identical, as setting the Caputo–Fabrizio derivative to zero yields the same algebraic equations. This follows from the property that 𝒟tαCFf(t)=0 if and only if f(t) is constant. Let (S∗,L∗,B∗,R∗) denote an equilibrium point. Then:

μ−βS∗(L∗+B∗)−μS∗=0,(5)

βS∗(L∗+B∗)−(ϵ+μ)L∗=0,(6)

ϵL∗−(γ+μ)B∗=0,(7)

γB∗−(ϱ+μ)R∗=0.(8)

From (7), we have B∗=ϵγ+μL∗. Let k=ϵγ+μ; then B∗=kL∗. Substituting into (6) gives:

βS∗L∗(1+k)=(ϵ+μ)L∗.

If L∗=0, then B∗=0, and from (5) we obtain S∗=1. Eq. (8) then yields R∗=0. This defines the disease-free equilibrium (DFE):

E0=(1,0,0,0).

If L∗≠0, we obtain from (6):

S∗=ϵ+μβ(1+k).

Substituting into (5) and solving for L∗ gives:

L∗=μϵ+μ(1−1ℛ0),

where the basic reproduction number ℛ0 is defined as:

ℛ0=β(1+k)ϵ+μ=β(ϵ+γ+μ)(ϵ+μ)(γ+μ).

Consequently, the endemic equilibrium (EE) is:

S∗=1ℛ0,L∗=μϵ+μ(1−1ℛ0),B∗=kL∗=ϵγ+μ⋅μϵ+μ(1−1ℛ0),R∗=γϱ+μB∗=γϱ+μ⋅ϵγ+μ⋅μϵ+μ(1−1ℛ0).

The endemic equilibrium exists and is biologically meaningful (positive) if and only if ℛ0>1.

2.4 Stability Analysis for Integer-Order Model

For the integer-order system (1), the stability properties are well-established in epidemic modeling theory [26]:

1.   Local Asymptotic Stability (LAS): The DFE E0 is LAS when ℛ0<1 and unstable when ℛ0>1. The endemic equilibrium E∗ is LAS when ℛ0>1, provided all parameters are positive.

2.   Global Asymptotic Stability (GAS): For the integer-order model, stronger stability results can be established:

•   The DFE is GAS on the feasible region Ω={(S,L,B,R)∈R+4:S+L+B+R≤1} when ℛ0≤1. This can be proven using the Lyapunov function

V(S,L,B,R)=L+ϵ+μϵB,

which satisfies V˙≤(ϵ+μ)(ℛ0S−1)L≤0 when ℛ0≤1.

•   The endemic equilibrium E∗ is GAS on Ω∖{E0} when ℛ0>1. This can be established constructing a suitable Lyapunov function of the form

V=S−S∗ln⁡S+L−L∗ln⁡L+c1(B−B∗ln⁡B)+c2(R−R∗ln⁡R),

with appropriate constants c1,c2>0.

3.   Global Exponential Asymptotic Stability (GEAS): Under the condition ℛ0>1, the endemic equilibrium exhibits exponential convergence rates. Specifically, there exist constants M>0 and δ>0 such that

‖X(t)−E∗‖≤Me−δt‖X(0)−E∗‖,

for all initial conditions in Ω, where X(t)=(S(t),L(t),B(t),R(t)).

Global Exponential Asymptotic Stability (GEAS):

Theorem 1: For the integer-order system (1), under the condition ℛ0>1, the endemic equilibrium E∗=(S∗,L∗,B∗,R∗) is globally exponentially asymptotically stable on the feasible region Ω∖{E0}. This means there exist constants M>0 and δ>0 such that:

‖X(t)−E∗‖≤Me−δt‖X(0)−E∗‖,

for all t≥0 and for any initial condition X(0)∈Ω, where X(t)=(S(t),L(t),B(t),R(t)).

Proof: We employ the Lyapunov direct method. Consider the candidate Lyapunov function:

V(X)=S−S∗ln⁡S+L−L∗ln⁡L+ϵ+μϵ(B−B∗ln⁡B)+γ(ϵ+μ)ϵ(ϱ+μ)(R−R∗ln⁡R).

Define the deviations:

x=S−S∗,y=L−L∗,z=B−B∗,w=R−R∗.

The time derivative of V along solutions of (1) is:

V˙=(1−S∗S)S˙+(1−L∗L)L˙+ϵ+μϵ(1−B∗B)B˙+γ(ϵ+μ)ϵ(ϱ+μ)(1−R∗R)R˙.

Substituting the expressions for S˙,L˙,B˙,R˙ from (1) and using the equilibrium conditions:

μ=βS∗(L∗+B∗)+μS∗,βS∗(L∗+B∗)=(ϵ+μ)L∗,ϵL∗=(γ+μ)B∗,γB∗=(ϱ+μ)R∗,

we obtain after algebraic simplification:

V˙=−μ(S−S∗)2S−βS∗(L∗+B∗)[S∗S+S(L+B)S∗(L∗+B∗)−2]−(ϵ+μ)(L−L∗L)2L∗L.

Using the inequality a+b≥2ab for a,b>0, we have:

S∗S+S(L+B)S∗(L∗+B∗)≥2L+BL∗+B∗.

Thus,

V˙≤−μ(S−S∗)2S−(ϵ+μ)(L−L∗)2L.

Since S(t),L(t)∈(0,1] for all t>0, there exists a constant c>0 such that:

V˙≤−c((S−S∗)2+(L−L∗)2)≤−c‖X−E∗‖2,

where the last inequality follows from the fact that B and R are linearly dependent on L at equilibrium.

Now, note that V is positive definite and radially unbounded in Ω. Moreover, there exist constants k1,k2>0 such that:

k1‖X−E∗‖2≤V(X)≤k2‖X−E∗‖2.

Combining with V˙≤−c‖X−E∗‖2, we obtain:

V˙≤−ck2V.

By the Gronwall inequality, this implies:

V(X(t))≤V(X(0))e−δt,

with δ=c/k2. Consequently,

‖X(t)−E∗‖≤V(X(t))k1≤V(X(0))k1e−δt/2≤Me−δt/2‖X(0)−E∗‖,

where M=k2/k1. This completes the proof of global exponential asymptotic stability. ◻

This proof follows the Lyapunov function approach commonly used in epidemic models with mass-action incidence. For similar stability analyses in integer-order compartmental models, see [26,47].

For the Caputo–Fabrizio fractional-order system (3), stability analysis requires special consideration due to the memory effects inherent in the fractional derivative. The linearized stability conditions differ from the integer-order case:

1.   Local Asymptotic Stability: For a fractional-order system with Caputo–Fabrizio derivative of order α∈(0,1), an equilibrium point is LAS if all eigenvalues λi of the Jacobian matrix evaluated at the equilibrium satisfy [7]:

Re(λi)<0and|arg⁡(λi)|>απ2.

For the DFE E0, the Jacobian matrix is:

J(E0)=(−μ−β−β00β−(ϵ+μ)β00ϵ−(γ+μ)000γ−(ϱ+μ)).

The eigenvalues are λ1=−μ, λ2=−(ϱ+μ), and the remaining two are determined by the submatrix:

M=(β−(ϵ+μ)βϵ−(γ+μ)).

The characteristic equation of M is λ2+a1λ+a2=0 with

a1=(ϵ+μ)+(γ+μ)−β,a2=(ϵ+μ)(γ+μ)−β(ϵ+γ+μ)=(ϵ+μ)(γ+μ)(1−ℛ0).

When ℛ0<1, both a1>0 and a2>0, ensuring all eigenvalues have negative real parts. Moreover, since the fractional order α∈(0,1), the condition |arg⁡(λi)|>απ2 is automatically satisfied for eigenvalues with negative real parts. Thus, the DFE is LAS for the fractional-order system when ℛ0<1.

2.   Global Stability: Global stability analysis for fractional-order systems is more challenging. However, using fractional Lyapunov direct methods [47], it can be shown that the DFE is globally asymptotically stable when ℛ0≤1. For the endemic equilibrium, numerical evidence suggests global stability when ℛ0>1, though rigorous analytical proofs remain an open research problem for this specific model.

3.   Memory Effects on Stability: The fractional order α influences the convergence rate to equilibrium. Smaller values of α (stronger memory effects) typically lead to slower convergence, as evidenced by our numerical simulations in Section 2.9. This contrasts with the integer-order case where convergence is exponential. The fractional-order system exhibits algebraic decay toward equilibrium when memory effects are strong (α close to 0), transitioning to exponential decay as α→1.

For comparative analysis, see Table 2:

Remark 2: The fractional-order model preserves the same epidemic threshold ℛ0 as the integer-order model, but the memory effects encoded by α modify the transient dynamics and convergence rates. This has practical implications for cybersecurity: while the threshold for virus extinction remains unchanged, the time to achieve eradication or control is prolonged under strong memory effects (small α), necessitating longer-term intervention strategies.

Remark 3: For the Caputo–Fabrizio derivative, the stability condition |arg⁡(λi)|>απ2 is less restrictive than for the classical Caputo derivative when α is small, due to the non-singular kernel. This mathematical property aligns with the physical interpretation that memory effects in cyber-epidemics are typically finite and fading rather than singular.

2.5 Existence and Uniqueness of Solutions

To establish the well-posedness of the Caputo–Fabrizio fractional computer virus model, we reformulate the system in terms of nonlinear mappings and apply fixed-point theory.

Let us define the nonlinearities associated with the model:

ω1(S,L,B,R,t)=μ−βS(L+B)−μS,ω2(S,L,B,R,t)=βS(L+B)−(ϵ+μ)L,ω3(S,L,B,R,t)=ϵL−(γ+μ)B,ω4(S,L,B,R,t)=γB−(ϱ+μ)R.

Let Λ(α) denote the normalization constant:

Λ(α)=1−α+αΓ(1+α).

The associated integral operator is:

𝒟tαCFf(t)=1−αΛ(α)∫0tf(ϕ)dϕ+αΛ(α)∫0tf(ϕ)(t−ϕ)dϕ.

Theorem 2 (Lipschitz Continuity): Consider the nonlinear functions ωi defined in (3). Under the boundedness established in Lemma 1, each ωi is Lipschitz continuous on the domain Ω with respect to its corresponding state variable.

Proof: By Lemma 1, all state variables are bounded by M=1 on [0,T]. Let S1,S2,L1,L2,B1,B2,R1,R2∈[0,1]. We examine each nonlinearity separately.

For ω1:

‖ω1(S1,L,B,R,t)−ω1(S2,L,B,R,t)‖=|β(L+B)(S2−S1)+μ(S2−S1)|≤(β⋅2M+μ)‖S1−S2‖=(β⋅2+μ)‖S1−S2‖=L1‖S1−S2‖,

where L1=2β+μ.

For ω2:

‖ω2(S,L1,B,R,t)−ω2(S,L2,B,R,t)‖=|βS(L2−L1)−(ϵ+μ)(L2−L1)|≤(βM+ϵ+μ)‖L1−L2‖=(β+ϵ+μ)‖L1−L2‖=L2‖L1−L2‖,

where L2=β+ϵ+μ.

For ω3:

‖ω3(L,B1,t)−ω3(L,B2,t)‖=|(γ+μ)(B2−B1)|=(γ+μ)‖B1−B2‖=L3‖B1−B2‖,

where L3=γ+μ.

For ω4:

‖ω4(B,R1,t)−ω4(B,R2,t)‖=|(ϱ+μ)(R2−R1)|=(ϱ+μ)‖R1−R2‖=L4‖R1−R2‖,

where L4=ϱ+μ.

Thus, each ωi is Lipschitz continuous with Lipschitz constants L1=2β+μ, L2=β+ϵ+μ, L3=γ+μ, and L4=ϱ+μ, respectively. ◻

Theorem 3 (Existence and Uniqueness): Let the initial conditions satisfy (S(0),L(0),B(0),R(0))∈Ω, and let each nonlinearity ωi satisfy the Lipschitz condition established in Theorem 2. Then the Caputo-Fabrizio fractional-order virus model (3) admits a unique solution on [0,T] for any finite T>0.

Proof: We reformulate the system (3) in integral form using the associated CF integral operator defined in Section 2.1:

𝒫S(t)=S(0)+ItαCFω1(S,L,B,R,t),𝒬L(t)=L(0)+ItαCFω2(S,L,B,R,t),ℛB(t)=B(0)+ItαCFω3(S,L,B,R,t),𝒮R(t)=R(0)+ItαCFω4(S,L,B,R,t).

By Lemma 1, all state variables remain in [0,M] with M=1. The Lipschitz continuity (Theorem 2) ensures that each operator is well-defined. We now establish that these operators are contractions.

Consider the operator 𝒫. For any S1,S2∈C([0,T],[0,M]), we have:

‖𝒫S1(t)−𝒫S2(t)‖≤1−αΛ(α)L1‖S1−S2‖∞+αΛ(α)∫0tL1‖S1−S2‖∞dϕ≤(1−αΛ(α)L1+αTΛ(α)L1)‖S1−S2‖∞=L1(1−α+αTΛ(α))‖S1−S2‖∞.

Define C1=L1(1−α+αTΛ(α)). If C1<1, then 𝒫 is a contraction. Similar calculations yield constants C2,C3,C4 for operators 𝒬,ℛ,𝒮, respectively. By choosing T sufficiently small, we can ensure all Ci<1. Consequently, by Banach’s Fixed Point Theorem, each operator admits a unique fixed point, which corresponds to a unique local solution on [0,T]. Since T is arbitrary and the boundedness holds for all time, the solution can be extended uniquely to any finite interval [0,T] using standard continuation arguments [2].

For uniqueness, assume two solutions S1(t) and S2(t) exist. Then:

‖S1(t)−S2(t)‖≤C1‖S1−S2‖∞.

Since C1<1 for sufficiently small t, we conclude S1(t)=S2(t). Similar reasoning applies to L, B, R. Hence, the system admits a unique solution on [0,T]. ◻

2.6 Hyers-Ulam Stability

The concept of Hyers–Ulam (HU) stability offers a valuable framework for evaluating the resilience of solutions to fractional differential equations in the presence of minor perturbations. In the framework of the Caputo–Fabrizio (CF) fractional computer virus model, HU stability guarantees that the approximate trajectories of the system remain proximal to the exact solutions when the initial conditions or model functions experience minor perturbations. This trait is crucial in digital epidemic modelling, as uncertainties in parameter values or measurement inaccuracies might influence system dynamics.

Definition 1 (Hyers–Ulam Stability): ([48]) Consider the CF system

{DtαCFS(t)=ω1(S,L,B,R,t),DtαCFL(t)=ω2(S,L,B,R,t),DtαCFB(t)=ω3(S,L,B,R,t),DtαCFR(t)=ω4(S,L,B,R,t),

with 0<α<1. Let X(t)=(S(t),L(t),B(t),R(t)) and define F(X(t),t)=(ω1,ω2,ω3,ω4)⊤. The system is said to be Hyers–Ulam stable on [0,T] if for every ε>0 and for every approximate solution X~(t)=(S~(t),L~(t),B~(t),R~(t)) satisfying

‖ItαCFX~(t)−F(X~(t),t)‖≤ε,

there exists an exact solution X(t) of the system such that

‖X~(t)−X(t)‖≤Cε,∀t∈[0,T],

where C>0 is a constant independent of ε.

Theorem 4 (Hyers–Ulam Stability of the CF Virus Model): Let ωi, i=1,…,4, be Lipschitz continuous with Lipschitz constants Li as established in Theorem 1. Then the CF fractional computer virus model (1) is Hyers–Ulam stable on [0,T].

Proof: From Theorem 2, the CF system admits a unique solution under the Lipschitz condition. Suppose X~(t) is an approximate solution with perturbation bounded by ε. Then for the exact solution X(t) we have

ItαCF(X~(t)−X(t))=F(X~(t),t)−F(X(t),t)+δ(t),

where ‖δ(t)‖≤ε. Using the Lipschitz continuity of F, it follows that

‖ItαCF(X~(t)−X(t))‖≤L‖X~(t)−X(t)‖+ε,

with L=max{L1,L2,L3,L4}. By applying the fractional Grönwall inequality adapted to the CF operator, one obtains

‖X~(t)−X(t)‖≤ε1−L(eLt−1)≤Cε,t∈[0,T],

for some constant C>0. This proves the Hyers–Ulam stability of the system. ◻

Remark 4: The above result guarantees that the CF fractional computer virus model is robust with respect to perturbations. In practical terms, even if parameter estimation or data acquisition introduces small errors, the model’s trajectories remain close to the true dynamics. This provides additional confidence in the applicability of the model for real-world cybersecurity analysis.

2.7 Numerical Scheme

To achieve this, we transform the Equations in (1) by applying the fundamental theorem of calculus, following the method outlined in [9] to

S⁡(ϕ)−S⁡(0)=1−αM(α)ℱ1(S,L,B,R,ϕ)+αM(α)∫0tℱ1(S,L,B,R,σ)dσ,L⁡(ϕ)−L⁡(0)=1−αM(α)ℱ2(S,L,B,R,ϕ)+αM(α)∫0tℱ2(S,L,B,R,σ)dσ,B⁡(ϕ)−B⁡(0)=1−αM(α)ℱ3(S,L,B,R,ϕ)+αM(α)∫0tℱ3(S,L,B,R,σ)dσ,R⁡(ϕ)−R⁡(0)=1−αM(α)ℱ4(S,L,B,R,ϕ)+αM(α)∫0tℱ4(S,L,B,R,σ)dσ.

At ϕn+1, we obtain:

S⁡(ϕn+1)−S⁡(0)=1−αM(α)ℱ1(Sn,Ln,Bn,Rn,ϕn)+αM(α)∫0ϕn+1ℱ1(S,L,B,R,ϕ)dϕ,L⁡(ϕn+1)−L⁡(0)=1−αM(α)ℱ2(Sn,Ln,Bn,Rn,ϕn)+αM(α)∫0ϕn+1ℱ2(S,L,B,R,ϕ)dϕ,B⁡(ϕn+1)−B⁡(0)=1−αM(α)ℱ3(Sn,Ln,Bn,Rn,ϕn)+αM(α)∫0ϕn+1ℱ3(S,L,B,R,ϕ)dϕ,R⁡(ϕn+1)−R⁡(0)=1−αM(α)ℱ4(Sn,Ln,Bn,Rn,ϕn)+αM(α)∫0ϕn+1ℱ4(S,L,B,R,ϕ)dϕ.

Similarly, at ϕn, we have:

S⁡(ϕn)−S⁡(0)=1−αM(α)ℱ1(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)+αM(α)∫0ϕnℱ1(S,L,B,R,ϕ)dϕ,L⁡(ϕn)−L⁡(0)=1−αM(α)ℱ2(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)+αM(α)∫0ϕnℱ2(S,L,B,R,ϕ)dϕ,B⁡(ϕn)−B⁡(0)=1−αM(α)ℱ3(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)+αM(α)∫0ϕnℱ3(S,L,B,R,ϕ)dϕ,R⁡(ϕn)−R⁡(0)=1−αM(α)ℱ4(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)+αM(α)∫0ϕnℱ4(S,L,B,R,ϕ)dϕ.

Subtracting the two equations results:

S⁡(ϕn+1)−S⁡(ϕn)=1−αM(α)[ℱ1(Sn,Ln,Bn,Rn,ϕn)−ℱ1(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)]+αM(α)∫ϕnϕn+1ℱ1(S,L,B,R,ϕ)dϕ,L⁡(ϕn+1)−L⁡(ϕn)=1−αM(α)[ℱ2(Sn,Ln,Bn,Rn,ϕn)−ℱ2(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)]+αM(α)∫ϕnϕn+1ℱ2(S,L,B,R,ϕ)dϕ,B⁡(ϕn+1)−B⁡(ϕn)=1−αM(α)[ℱ3(Sn,Ln,Bn,Rn,ϕn)−ℱ3(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)]+αM(α)∫ϕnϕn+1ℱ3(S,L,B,R,ϕ)dϕ,R⁡(ϕn+1)−R⁡(ϕn)=1−αM(α)[ℱ4(Sn,Ln,Bn,Rn,ϕn)−ℱ4(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)]+αM(α)∫ϕnϕn+1ℱ4(S,L,B,R,ϕ)dϕ,

where the integral terms are:

∫ϕnϕn+1ℱ1(S,L,B,R,ϕ)dϕ=3h2ℱ1(Sn,Ln,Bn,Rn,ϕn)−h2ℱ1(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1),∫ϕnϕn+1ℱ2(S,L,B,R,ϕ)dϕ=3h2ℱ2(Sn,Ln,Bn,Rn,ϕn)−h2ℱ2(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1),∫ϕnϕn+1ℱ3(S,L,B,R,ϕ)dϕ=3h2ℱ3(Sn,Ln,Bn,Rn,ϕn)−h2ℱ3(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1),∫ϕnϕn+1ℱ4(S,L,B,R,ϕ)dϕ=3h2ℱ4(Sn,Ln,Bn,Rn,ϕn)−h2ℱ4(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1).

Thus, we obtain:

S⁡(ϕn+1)−S⁡(ϕn)=1−αM(α)[ℱ1(Sn,Ln,Bn,Rn,ϕn)−ℱ1(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)]+3αh2M(α)ℱ1(Sn,Ln,Bn,Rn,ϕn)−αh2M(α)ℱ1(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1),L⁡(ϕn+1)−L⁡(ϕn)=1−αM(α)[ℱ2(Sn,Ln,Bn,Rn,ϕn)−ℱ2(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)]+3αh2M(α)ℱ2(Sn,Ln,Bn,Rn,ϕn)−αh2M(α)ℱ2(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1),B⁡(ϕn+1)−B⁡(ϕn)=1−αM(α)[ℱ3(Sn,Ln,Bn,Rn,ϕn)−ℱ3(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)]+3αh2M(α)ℱ3(Sn,Ln,Bn,Rn,ϕn)−αh2M(α)ℱ3(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1),R⁡(ϕn+1)−R⁡(ϕn)=1−αM(α)[ℱ4(Sn,Ln,Bn,Rn,ϕn)−ℱ4(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1)]+3αh2M(α)ℱ4(Sn,Ln,Bn,Rn,ϕn)−αh2M(α)ℱ4(Sn−1,Ln−1,Bn−1,Rn−1,ϕn−1).