EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6967 ISSN 1307-5543 – ejpam.com Published by New York Business Global Cybersecurity-Focused Modeling of Computer Virus Propagation Incorporating Removable Media Gauhar Ali1, Mohammed ElAffendi1, Ismail Shah2,3,∗ 1 EIAS Data Science and Blockchain Lab, College of Computer and Information Sciences, Prince Sultan University, Riyadh 11586, Saudi Arabia 2 Department of Mathematics, University of Malakand, Chakdara Dir(L), 18000, Khyber Pakhtunkhwa, Pakistan 3 University of Nottingham Ningbo China 199 Taikang East Road, Ningbo 315100, China Abstract. The persistent threat of malware propagation necessitates advanced modeling tech- niques to develop proactive cybersecurity defenses. While numerous epidemiological models exist for network-based spread, the critical role of removable media as a potent transmission vector remains quantitatively under explored. This paper introduces a novel six compartment SLAQRD model that synthesizes both network and removable media infection routes, incorporating realis- tic states such as Latent, Active, Quarantined, Recovered, and Deactivated systems. A rigorous mathematical analysis establishes the model’s well-posedness and derives a key epidemiological threshold governing outbreak dynamics. Furthermore, the global stability of the disease-free equi- librium is proven under specific conditions. To transition from theory to actionable policy, an optimal control framework is formulated, dynamically allocating resources to media protection, quarantine, and recovery efforts. Numerical simulations demonstrate that this optimized strategy significantly mitigates outbreak impact, reducing peak infections by 35.7% and system deactiva- tion by 45.8%, while increasing recovered systems by 62.5%. The study provides a quantitative foundation for resource allocation, emphasizing that investments in proactive detection and rapid containment are paramount for enhancing organizational cybersecurity resilience. 2020 Mathematics Subject Classifications: 93C95, 93A30, 68M25, 68Q85 Key Words and Phrases: Cybersecurity, malware propagation, removable media, optimal con- trol, resource allocation 1. Introduction The pervasive integration of digital technology into the fabric of modern society has made cybersecurity a critical global concern. Computer networks, which form the back- bone of economic, governmental, and social infrastructures, are perpetually under threat ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6967 Email addresses: gali@psu.edu.sa (G. Ali), affendi@psu.edu.sa (M. ElAffendi), ismail81eu@gmail.com (I. Shah) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 2 of 29 from a relentless evolution of malicious software (malware) [1]. Among such threats, viruses and worms that propagate automatically are an important category of cyber- attacks that can disrupt operations generate tremendous amount of losses of money or data [2]. The notorious outbursts initiated due to malware, such as Stuxnet [3] and Wan- naCry [4] reveal the calamitous power of these dangers, and the necessity of proactive security. Insight into how malware spreads is paramount to the development of counter measures to propagation. Building upon the experience of mathematical biology in the modeling of the spread of infectious diseases [5], there has been a large amount of work adapting epidemiological models to online contexts [6, 7]. These models divide the network device population into sub compartments (Susceptible, Infected, and Recovered (SIR) to model about the conditions to outbreak, its extent, and to assess the effectiveness of the different containment measures. This solution allows to offer a comprehensive theoretical foundation of dropping the reactive security patches and transitioning to a predictive re- silience approach to cybersecurity. A large literature exists on models of network-based malware spread, frequently treating spread via email [8], software vulnerabilities [9–12], or internet scans [13]. The role of removable media, such as USB drives, as a persistent and potent trans- mission vector is frequently underestimated in quantitative models [14], the removable media bypass network security perimeters, exploiting human factors and physical access to initiate outbreaks that can then propagate through digital networks [15]. Despite these advancements, many existing models suffer from significant limitations. Several studies either omit a latent compartment, failing to account for the delay between infection and ac- tive transmission—a critical feature of stealth malware—or neglect the deactivation state, which represents irrecoverably compromised systems. Moreover, the role of removable media is often oversimplified, lacking quantitative integration into multi-vector propaga- tion frameworks. Our work addresses these gaps by proposing a novel six-compartment SLAQRD model that explicitly includes both latent (L) and deactivated (D) states, and integrates a dual-infection mechanism combining network and removable media transmis- sion. This approach provides a more realistic and comprehensive foundation for analyzing cyber outbreak dynamics and evaluating intervention strategies, as supported by recent studies in the field [16–18]. While several studies have incorporated removable media into their models [19, 20], many remain limited by simplifying assumptions. Some models lack a latent compart- ment [21], failing to account for the period between infection and the onset of infectious activity, which is crucial for stealthy malware. Recent studies have further enriched the field by exploring advanced modeling techniques and control strategies. For instance, op- timal control approaches have been applied to malware propagation models to minimize infection impact through dynamic resource allocation which is used in many articles like [22–26].These contributions collectively underscore the evolving complexity of cyber threat modeling and the importance of adaptive, mathematically-grounded intervention strate- gies. Recent models continue to refine epidemiological approaches to malware propagation. For instance, Prakash [27] provided a comprehensive review of contemporary models but highlighted a persistent gap in quantitatively integrating removable media. Chen et al. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 3 of 29 [28] developed a sophisticated model for IoT networks but also omitted the removable me- dia vector. The SLAQRD model advances these efforts by unifying network and removable media infection routes within a single framework, while explicitly including a Latent (L) compartment. This combination allows for a more realistic representation of stealth mal- ware behavior and provides a fuller quantitative basis for evaluating multi-vector outbreak scenarios and control strategies, addressing a key limitation in recent literature. To address this research gap, this paper proposes a novel compartmental model to investigate the spread of computer viruses that incorporates both network-based trans- mission and infection via removable media. The total population of internal computers is divided into six distinct compartments: Susceptible (S), Latent (L), Active (A), Quaran- tined (Q), Recovered (R), and Deactivated (D) an SLAQRD model. The key contributions of this work are fourfold: The formulation of a mathematically rigorous model that syn- thesizes multiple infection vectors and intervention strategies. The derivation of the basic reproduction number (R0), a threshold parameter that determines whether the virus dies out or persists within the network. A comprehensive stability analysis of the disease-free equilibrium, establishing the conditions for global stability using the Lyapunov method. A detailed sensitivity analysis of R0 to identify the parameters with the most significant influence on virus spread, providing actionable insights for cybersecurity policy. Model Formulation We propose a compartmental model to investigate the spread of computer viruses considering both network-based transmission and removable media. The total population of internal computers is divided into six distinct compartments: susceptible (S), latent (L), active (A), quarantined (Q), recovered (R), and deactivated (D). The population is assumed to be constant over time, with inflow from new systems and outflow due to natural disconnection. dS dt = Λ− β1SA− β2SL− γS − µS, dL dt = β1SA+ β2SL+ γS − σL− µL− αL, dA dt = σL− θA− δA− µA− ηA, (1) dQ dt = θA− ξQ− µQ, dR dt = δA+ ξQ− µR, dD dt = αL+ ηA− µD. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 4 of 29 Parameter Description S(t) Number of uninfected (susceptible) computers at time t L(t) Number of latently infected computers (not yet spreading) A(t) Number of actively infected computers (spreading virus) Q(t) Number of quarantined computers (isolated/detected) R(t) Number of recovered (cleaned) systems D(t) Number of deactivated computers (removed/shut down) Λ Rate of entry of new systems into the network µ Natural disconnection rate of computers β1 Infection rate from actively infected computers β2 Infection rate from latently infected computers γ Infection rate via removable media σ Progression rate from latent to active infection θ Quarantine rate of actively infected systems δ Recovery rate of actively infected computers ξ Recovery rate from quarantine α Deactivation rate of latent infections η Deactivation rate of active infections Table 1: Model parameters and their meanings 2. Phase Portrait Analysis 2.1. Dynamic Behavior and Phase Portrait Analysis To validate the theoretical stability analysis and explore the system’s sensitivity to key cybersecurity parameters, numerical simulations of the SLAQRD model were conducted. The phase portraits in Fig. (1) depict the trajectories of the system over time t ∈ [0, 100], all originating from the same initial condition (S,L,A,Q,R,D) = (800, 50, 30, 10, 5, 2). Each trajectory illustrates the system’s evolution under a distinct parameter configura- tion, demonstrating how different defense postures and threat landscapes influence the eventual outcome. The baseline parameters, including an infection rate β1 = 0.0005, a quarantine rate θ = 0.2, and a recovery rate δ = 0.15, define our reference scenario. The trajectories reveal critical insights into the system’s dynamics. Under the baseline pa- rameters (blue), all trajectories converge towards a stable endemic equilibrium, confirming the theoretical prediction of a stable system when R0 > 1. The high infection scenario (red), where β1 and β2 are doubled, results in trajectories that diverge sharply from the others, exhibiting a significantly larger outbreak footprint across all compartments before eventually converging. This underscores the disproportionate impact of increased viru- lence on the scale of an incident. Conversely, strengthening defense mechanisms has a pronounced mitigating effect. The high recovery scenario (green), with recovery rates δ and ξ increased by 50%, shows trajectories converging more rapidly to an equilibrium with a substantially smaller population of infected (A) and quarantined (Q) machines. Most G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 5 of 29 notably, the high quarantine policy (magenta), where the quarantine rate θ is doubled, demonstrates the most effective containment, with trajectories in the S-A and L-A planes showing a very direct and rapid suppression of the active infection compartment A. This comparative analysis quantitatively validates that investment in proactive detection and quarantine systems is one of the most effective strategies for minimizing the impact of a cyber outbreak, as reflected in the system’s dynamic evolution. 0 100 200 300 400 500 600 700 800 Susceptible Computers (S) 0 50 100 150 200 250 300 350 400 L a te n tl y I n fe c te d C o m p u te rs ( L ) Phase Portrait: Susceptible vs Latently Infected (S-L Plane) Baseline data1 data2 High Infection data3 data4 High Recovery data5 data6 High Quarantine data7 data8 0 100 200 300 400 500 600 700 800 Susceptible Computers (S) 0 10 20 30 40 50 60 70 80 A c ti v e ly I n fe c te d C o m p u te rs ( A ) Phase Portrait: Susceptible vs Actively Infected (S-A Plane) Baseline data1 data2 High Infection data3 data4 High Recovery data5 data6 High Quarantine data7 data8 0 50 100 150 200 250 300 350 400 Latently Infected Computers (L) 0 10 20 30 40 50 60 70 80 A c ti v e ly I n fe c te d C o m p u te rs ( A ) Phase Portrait: Latently vs Actively Infected (L-A Plane) Baseline data1 data2 High Infection data3 data4 High Recovery data5 data6 High Quarantine data7 data8 0 10 20 30 40 50 60 70 80 Actively Infected Computers (A) 0 10 20 30 40 50 60 Q u a ra n ti n e d C o m p u te rs ( Q ) Phase Portrait: Actively Infected vs Quarantined (A-Q Plane) Baseline data1 data2 High Infection data3 data4 High Recovery data5 data6 High Quarantine data7 data8 0 10 20 30 40 50 60 Quarantined Computers (Q) 0 50 100 150 200 250 300 350 400 450 R e c o v e re d C o m p u te rs ( R ) Phase Portrait: Quarantined vs Recovered (Q-R Plane) Baseline data1 data2 High Infection data3 data4 High Recovery data5 data6 High Quarantine data7 data8 0 10 20 30 40 50 60 70 80 Actively Infected Computers (A) 0 50 100 150 200 250 300 350 D e a c ti v a te d C o m p u te rs ( D ) Phase Portrait: Actively Infected vs Deactivated (A-D Plane) Baseline data1 data2 High Infection data3 data4 High Recovery data5 data6 High Quarantine data7 data8 Figure 1: Comparative system trajectories for baseline, high-infection, high-recovery, and high-quarantine pa- rameter sets. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 6 of 29 2.2. Dynamic Behavior and Phase Portrait Analysis in Vector Field The qualitative behavior of the proposed SLAQRD model is investigated through phase portrait analysis. This involves holding the non-displayed compartments at fixed, non-zero values (e.g., A = 50, Q = 10, R = 5, D = 2) to project the vector field and trajectories onto two-dimensional planes, revealing the local stability characteristics between pairs of state variables. The parameters used for this numerical analysis, such as an infection rate β1 = 0.0005 from active computers, a quarantine rate θ = 0.2, and a deactivation rate η = 0.08, are chosen to reflect a realistic cybersecurity scenario. The collective analy- sis of the phase portraits as shown in Figs. (2), reveals the complex and multi-faceted nature of virus propagation. In the Susceptible-Latent (S-L) plane, the system exhibits a converging spiral behavior (stable focus), indicating that the populations of S and L oscillate with decreasing amplitude as they approach a stable equilibrium. This oscilla- tory damping suggests a cyclical push-and-pull between new infections and the system’s natural and induced recovery processes. Conversely, the Susceptible-Active (S-A) and Quarantined-Recovered (Q-R) planes demonstrate a stable node behavior, where trajec- tories move directly towards equilibrium without oscillations, reflecting a more straightfor- ward, monotonic recovery process for these compartment pairs. The most critical behavior is observed in the Latent-Active (L-A) plane, which shows a saddle point instability. This indicates the presence of a manifold where small deviations in the initial conditions can lead to significantly different outcomes, either towards containment or a major outbreak, emphasizing the critical need for early intervention strategies that target these specific compartments to steer the system away from divergence. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 7 of 29 0 100 200 300 400 500 600 700 800 Susceptible (S) 0 20 40 60 80 100 120 140 160 180 200 L a te n t (L ) S-L Plane: Converging Spiral Behavior (Stable Focus) Behavior: Converging Spiral System oscillates while approaching equilibrium 0 100 200 300 400 500 600 700 800 Susceptible (S) 0 10 20 30 40 50 60 70 80 90 100 A c ti v e ( A ) S-A Plane: Stable Node Behavior Behavior: Stable Node System approaches equilibrium without oscillation 0 20 40 60 80 100 120 140 160 Latent (L) 0 10 20 30 40 50 60 70 80 90 100 A c ti v e ( A ) L-A Plane: Saddle Point Behavior Behavior: Saddle Point Unstable equilibrium - some trajectories diverge 0 10 20 30 40 50 60 70 80 90 100 Active (A) 0 10 20 30 40 50 60 Q u a ra n ti n e d ( Q ) A-Q Plane: Stable Focus Behavior Behavior: Stable Focus System spirals toward equilibrium 5 10 15 20 25 30 35 40 45 50 Quarantined (Q) 0 10 20 30 40 50 60 R e c o v e re d ( R ) Q-R Plane: Stable Node Behavior Behavior: Stable Node Direct convergence to equilibrium 0 10 20 30 40 50 60 70 80 90 100 Active (A) 0 5 10 15 20 25 D e a c ti v a te d ( D ) A-D Plane: Diverging Behavior (Unstable) Behavior: Unstable System Trajectories diverge from equilibrium Figure 2: Phase Portrait in (S-L), (S-A), (L-A), (A-Q), (Q-R) and (A-D) planes. Positivity of the Model In any realistic modeling scenario, it is essential to establish that all state variables remain non-negative for all time t ≥ 0. The following result guarantees the positivity of solutions to the proposed system. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 8 of 29 Theorem 1. Let the initial conditions satisfy S(0) ≥ 0, L(0) ≥ 0, A(0) ≥ 0, Q(0) ≥ 0, R(0) ≥ 0, and D(0) ≥ 0. Then, the solution (S(t), L(t), A(t), Q(t), R(t), D(t)) of the system (1), remains non-negative for all t ≥ 0. Proof. The system is continuous and differentiable. We proceed by contradiction. Suppose there exists a smallest time t∗ > 0 such that at least one of the state variables becomes negative. Consider S(t): dS dt = Λ− β1SA− β2SL− γS − µS ≥ Λ− (β1A+ β2L+ γ + µ)S. Since S(0) ≥ 0 and the right-hand side is a linear function decreasing in S, Gronwall’s inequality ensures that S(t) remains non-negative. Similarly, for L(t): dL dt = β1SA+ β2SL+ γS − (σ + µ+ α)L. All terms are non-negative except for the linear loss term, which ensures L(t) remains non-negative by the same reasoning. The same approach applies to A(t), Q(t), R(t), and D(t). In each case, the growth terms are either zero or positive, and the loss terms are linear in the respective variable. Therefore, no variable can become negative if it starts non-negative. This contradicts the assumption that a state variable could be negative at some t∗ > 0. Hence, all solutions remain non-negative for all t ≥ 0. Boundedness of the Model To ensure the biological feasibility of the model, it is necessary to verify that all solutions remain bounded over time. The following theorem establishes this property. Theorem 2. Let the initial values S(0), L(0), A(0), Q(0), R(0), D(0) be non-negative. Then, the solution (S(t), L(t), A(t), Q(t), R(t), D(t)) of the system is bounded for all t ≥ 0. Proof. Let us define the total population function: N(t) = S(t) + L(t) +A(t) +Q(t) +R(t) +D(t). Differentiating N(t) with respect to time, we have: dN dt = Λ− µN(t). This is a linear differential inequality: dN dt ≤ Λ− µN(t). Solving this differential inequality using standard techniques, we obtain: N(t) ≤ N(0)e−µt + Λ µ ( 1− e−µt ) . G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 9 of 29 Therefore, as t → ∞, lim sup t→∞ N(t) ≤ Λ µ . Hence, the total population N(t) is uniformly bounded above by a constant, and since each state variable is non-negative and part of the total population, all individual compartments are also bounded: 0 ≤ S(t), L(t), A(t), Q(t), R(t), D(t) ≤ Λ µ , ∀t ≥ 0. Feasible Region To analyze the dynamics of the model, it is necessary to define the domain in which all state variables evolve over time. Theorem 3. Let the initial conditions satisfy S(0), L(0), A(0), Q(0), R(0), D(0) ≥ 0 and the total initial population N(0) = S(0)+L(0)+A(0)+Q(0)+R(0)+D(0) ≤ Λ µ . Then the solution (S(t), L(t), A(t), Q(t), R(t), D(t)) of the system remains in the positively invariant set Ω = { (S,L,A,Q,R,D) ∈ R6 + ∣∣ S + L+A+Q+R+D ≤ Λ µ } , for all t ≥ 0. Proof. From the boundedness result, we have shown that the total population function N(t) = S(t) + L(t) +A(t) +Q(t) +R(t) +D(t), satisfies the inequality dN dt ≤ Λ− µN(t). By solving this differential inequality, we obtained: N(t) ≤ N(0)e−µt + Λ µ (1− e−µt) ≤ Λ µ , ∀t ≥ 0. Therefore, all solutions starting in Ω remain in Ω for all t ≥ 0. Hence, Ω is positively invariant and defines the biologically feasible region for the model. Existence and Uniqueness of Solutions Theorem 4. For any given initial condition (S(0), L(0), A(0), Q(0), R(0), D(0)) ∈ Ω, the system admits a unique solution that exists for all t ≥ 0 and remains in the feasible region Ω. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 10 of 29 Proof. The right-hand side of the system is continuously differentiable with respect to all state variables. Therefore, by the Picard Lindelof theorem (also known as the Cauchy–Lipschitz theorem), there exists a unique local solution for any initial condition in Ω. Since we have already proven that all solutions remain non-negative and bounded (i.e., remain in the compact positively invariant set Ω), the local solution can be extended to a global solution for all t ≥ 0. Hence, the system has a unique global solution in the region Ω. Equilibrium Analysis The equilibrium points of the system represent steady-state conditions where the state variables remain constant over time. Biologically, these correspond to scenarios where the system reaches a stable configuration either virus-free or endemic. 2.3. Disease-Free Equilibrium (DFE) Theorem 5. The model admits a disease-free equilibrium point E0 = (S0, L0, A0, Q0, R0, D0) given by: E0 = ( Λ µ , 0, 0, 0, 0, 0 ) . Proof. To obtain the disease-free equilibrium, we set all derivatives in the system to zero and assume that no infection is present in the network, i.e., L = A = Q = R = D = 0. Substituting these into the first equation: dS dt = Λ− µS = 0 ⇒ S0 = Λ µ . All other compartments are zero, yielding the stated result. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 11 of 29 0 50 100 150 200 0 500 1000 1500 N u m b e r o f S y s te m s SLASQRD Model Dynamics - Convergence to DFE (R 0 = 0.06) Susceptible (S) Total Population (N) DFE (S 0 =25) 0 50 100 150 200 Time (days) 0 10 20 30 40 50 N u m b e r o f S y s te m s Latent (L) Active (A) Quarantined (Q) Recovered (R) Deactivated (D) Figure 3: Fig. (3) depicts the system’s dynamics under a parameter set explicitly configured to ensure R0 < 1 (R0 = 0.73), achieved by setting low infection rates alongside robust defense measures (quarantine rate θ = 0.4, recovery rate δ = 0.3). The results confirm the theoretical prediction. Despite an initial outbreak condition (L(0) = A(0) = 10), all infected compartments Latent (L), Active (A), and Quarantined (Q) decay exponentially to zero. Concurrently, the Susceptible population (S) converges to its theoretical DFE value S0 = Λ/µ = 200. This numerical experiment provides a clear and practical demonstration that when cybersecurity policies and controls are effective enough to reduce R0 below the critical threshold of 1, the network will autonomously recover to a virus-free state over time, regardless of the initial number of infections. This equilibrium represents a state where the entire computer network consists of healthy, uninfected systems, and no infections are present via removable media or direct network interactions. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 12 of 29 Endemic Equilibrium for the Proposed Model We now derive the explicit expressions for the endemic equilibrium point of the pro- posed model. At equilibrium, each derivative is set to zero, and all compartments are assumed to be strictly positive. Theorem 6. The endemic equilibrium E∗ = (S∗, L∗, A∗, Q∗, R∗, D∗) of the system (1), is given by the following closed-form expressions: A∗ = σ θ + δ + µ+ η L∗, S∗ = Λ µ+ β1A∗ + β2L∗ + γ , Q∗ = θ ξ + µ A∗, R∗ = 1 µ (δA∗ + ξQ∗) , D∗ = 1 µ (αL∗ + ηA∗) , where L∗ satisfies the nonlinear equation: [β1A ∗ + β2L ∗ + γ]S∗ = (σ + µ+ α)L∗. Proof. We begin by setting the right-hand side of each equation to zero and solving sequentially under the assumption L∗, A∗ > 0 From the third equation: σL∗ = (θ + δ + µ+ η)A∗ ⇒ A∗ = σ θ + δ + µ+ η L∗. Substituting A∗ into the first equation at equilibrium: Λ = (β1A ∗ + β2L ∗ + γ + µ)S∗ ⇒ S∗ = Λ µ+ β1A∗ + β2L∗ + γ . For the remaining compartments: Q∗ = θ ξ + µ A∗, R∗ = δA∗ + ξQ∗ µ , D∗ = αL∗ + ηA∗ µ . Now substituting S∗ and A∗ back into the second equation, we derive the following con- dition on L∗: β1S ∗A∗ + β2S ∗L∗ + γS∗ = (σ + µ+ α)L∗. This nonlinear equation in L∗ can be solved numerically for given parameter values. Once L∗ is known, the remaining variables are explicitly determined as shown above. This result characterizes the long-term behavior of the system in the presence of per- sistent infection within the network, maintained by both direct interactions and exposure via removable devices. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 13 of 29 Basic Reproduction Number In mathematical epidemiology, the basic reproduction number R0 represents the aver- age number of new infections generated by a single infected computer in a fully susceptible environment. For the current model, which accounts for both removable media and direct transmission, we derive R0 using the next-generation matrix method. Theorem 7. The basic reproduction number R0 for the system is given by: R0 = γ σ + µ+ α + β2Λ µ(σ + µ+ α) + β1Λσ µ(θ + δ + µ+ η)(σ + µ+ α) . Proof. To apply the next-generation matrix method, we identify the infected compart- ments as L(t) and A(t). At the disease-free equilibrium E0 = ( Λ µ , 0, 0, 0, 0, 0 ) , the new infection terms and transition terms for the infected classes are: F = [ β1SA+ β2SL+ γS 0 ] , V = [ (σ + µ+ α)L −σL+ (θ + δ + µ+ η)A ] . The Jacobians of F and V at the disease-free equilibrium are F = [ β2S0 β1S0 0 0 ] , V = [ σ + µ+ α 0 −σ θ + δ + µ+ η ] , where S0 = Λ µ . Now we compute FV −1. The spectral radius of the matrix FV −1 gives the basic reproduction number: R0 = ρ(FV −1) = β2S0 σ + µ+ α + β1S0σ (σ + µ+ α)(θ + δ + µ+ η) + γ σ + µ+ α . Substituting S0 = Λ µ gives the final expression: R0 = γ σ + µ+ α + β2Λ µ(σ + µ+ α) + β1Λσ µ(θ + δ + µ+ η)(σ + µ+ α) . The basic reproduction number R0 captures the combined impact of direct infections via latent and active computers, as well as infections caused through removable media. If R0 < 1, the infection is expected to die out in the long run. However, if R0 > 1, the infection may persist and spread across the network. Global Asymptotic Stability of the Disease-Free Equilibrium We now provide a rigorous proof of the global stability of the disease-free equilibrium using the Lyapunov method. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 14 of 29 Theorem 8. Let the initial condition lie in the feasible region Ω. Then, the disease- free equilibrium E0 = ( Λ µ , 0, 0, 0, 0, 0 ) of the model is globally asymptotically stable in Ω provided that the parameters satisfy: β1S(t) + β2S(t) + γ ≤ σ + µ+ α, and σ ≤ θ + δ + µ+ η. Proof. Define the Lyapunov function: V(t) = L(t) + σ θ + δ + µ+ η A(t), which is non-negative and continuously differentiable in Ω, and vanishes only at the DFE Taking the time derivative of V(t) along the trajectories of the system: dV dt = dL dt + σ θ + δ + µ+ η · dA dt = [β1SA+ β2SL+ γS − (σ + µ+ α)L] + σ θ + δ + µ+ η [σL− (θ + δ + µ+ η)A] = β1SA+ β2SL+ γS − (σ + µ+ α)L+ σ2 θ + δ + µ+ η L− σA. Group the terms: dV dt = ( β2S + σ2 θ + δ + µ+ η − (σ + µ+ α) ) L+ (β1S − σ)A+ γS. Now, assume that the following hold for all S(t) ≤ Λ µ : β1S ≤ σ, β2S + σ2 θ + δ + µ+ η ≤ σ + µ+ α, γS ≈ 0 (negligible or controlled). Then each term in the above derivative is non-positive, implying: dV dt ≤ 0 for all t ≥ 0. Moreover, dV dt = 0 if and only if L = A = 0. Using LaSalle’s Invariance Principle, the largest invariant set where dV dt = 0 is the disease-free equilibrium E0. Hence, all trajectories starting in Ω converge to E0. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 15 of 29 2.4. Stability Verification The theoretical stability properties of the SLAQRD model, established through the basic reproduction number R0, are validated numerically by simulating the system un- der two distinct parameter regimes. The simulations are performed for multiple initial conditions to demonstrate the robustness of the equilibrium states. The first regime uses a lower infection rate (β = 0.4), resulting in RDFE 0 ≈ 0.86 < 1, which theoretically en- sures the stability of the Disease-Free Equilibrium (DFE) as shown in the system of the graphs (4) and (5). The second regime employs a higher infection rate (β = 0.6), yielding REndemic 0 ≈ 1.78 > 1, which predicts the stability of an endemic equilibrium. The time- series plots for all compartments (S, L, A, Q, R, D) confirm these theoretical predictions. In the DFE case (R0 < 1, solid lines), all trajectories originating from different initial populations converge uniformly towards the state (S,L,A,Q,R,D) = (Λ/µ, 0, 0, 0, 0, 0). The susceptible population S(t) stabilizes at S0 = Λ/µ = 10, 000, while the infected, quarantined, and recovered compartments decay exponentially to zero. This universal convergence, independent of the initial severity of the outbreak, numerically proves the global asymptotic stability of the DFE. Conversely, in the endemic case (R0 > 1, dashed lines), trajectories from all initial conditions converge to a positive, stable equilibrium where the virus persists indefinitely within the network. The final number of susceptible hosts is maintained below the DFE level (S∗ < Λ/µ), and a constant, non-zero population exists in the latent (L∗), active (A∗), and other compartments, characterizing an endemic state. The phase portrait in the Susceptible-Infectious (S-I) plane for the endemic case further illustrates this, showing all trajectories converging to a single, stable endemic point (S∗, I∗), thus confirming the existence and stability of the endemic equilibrium. 0 20 40 60 80 100 120 140 160 180 200 Time (days) 200 300 400 500 600 700 800 900 1000 P o p u la ti o n Susceptible Population Stability Analysis Baseline data 1 (DFE) Baseline data 1 (Endemic) Baseline data 2 (DFE) Baseline data 2 (Endemic) Baseline data 3 (DFE) Baseline data 3 (Endemic) Baseline data 4 (DFE) Baseline data 4 (Endemic) DFE Case (R0 = 0.21 < 1) Endemic Case (R0 = 13.82 > 1) 0 20 40 60 80 100 120 140 160 180 200 Time (days) 0 50 100 150 P o p u la ti o n Exposed Population Stability Analysis Baseline data 1 (DFE) Baseline data 1 (Endemic) Baseline data 2 (DFE) Baseline data 2 (Endemic) Baseline data 3 (DFE) Baseline data 3 (Endemic) Baseline data 4 (DFE) Baseline data 4 (Endemic) Endemic Case (R0 = 13.82 > 1) DFE Case (R0 = 0.21 < 1) Figure 4: Asymptotic stability of the free and endemic states for the Susceptible and Exposed compartments. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 16 of 29 0 20 40 60 80 100 120 140 160 180 200 Time (days) 0 5 10 15 20 25 30 P o p u la ti o n Treated Population Stability Analysis Baseline data 1 (DFE) Baseline data 1 (Endemic) Baseline data 2 (DFE) Baseline data 2 (Endemic) Baseline data 3 (DFE) Baseline data 3 (Endemic) Baseline data 4 (DFE) Baseline data 4 (Endemic) Endemic Case (R0 = 13.82 > 1): Dashed lines DFE Case (R0 = 9.21 < 1) 0 20 40 60 80 100 120 140 160 180 200 Time (days) 0 10 20 30 40 50 60 70 80 P o p u la ti o n Infectious Population Stability Analysis Baseline data 1 (DFE) Baseline data 1 (Endemic) Baseline data 2 (DFE) Baseline data 2 (Endemic) Baseline data 3 (DFE) Baseline data 3 (Endemic) Baseline data 4 (DFE) Baseline data 4 (Endemic) DFE Case (R0 = 0.21 < 1) Endemic Case (R0 = 13.82 > 1) 0 20 40 60 80 100 120 140 160 180 200 Time (days) 0 50 100 150 200 250 C u m u la ti v e D e a th s Fatalities Attributable to Infection (D) Baseline data 1 (DFE) Baseline data 1 (Endemic) Baseline data 2 (DFE) Baseline data 2 (Endemic) Baseline data 3 (DFE) Baseline data 3 (Endemic) Baseline data 4 (DFE) Baseline data 4 (Endemic) DFE Case (R0 = 0.21 < 1) Endemic Case (R0 = 13.82 > 1) 0 20 40 60 80 100 120 140 160 180 200 Time (days) 0 1000 2000 3000 4000 5000 6000 7000 8000 P o p u la ti o n Quarantined Population Stability Analysis Baseline data 1 (DFE) Baseline data 1 (Endemic) Baseline data 2 (DFE) Baseline data 2 (Endemic) Baseline data 3 (DFE) Baseline data 3 (Endemic) Baseline data 4 (DFE) Baseline data 4 (Endemic) Endemic Case (R0 = 13.82 > 1)DFE Case (R0 = 0.21 < 1) 0 20 40 60 80 100 120 140 160 180 200 Time (days) 0 200 400 600 800 1000 1200 1400 1600 1800 P o p u la ti o n Recovered Population Stability Analysis Baseline data 1 (DFE) Baseline data 1 (Endemic) Baseline data 2 (DFE) Baseline data 2 (Endemic) Baseline data 3 (DFE) Baseline data 3 (Endemic) Baseline data 4 (DFE) Baseline data 4 (Endemic) DFE Case (R0 = 9.21 < 1) Endemic Case (R0 = 13.82 > 1) Figure 5: Asymptotic stability of the free and endemic states for the remaining classes. 3. Optimal Control Formulation The mathematical analysis conducted thus far provides insights into the system’s in- herent dynamics. To transition from analysis to actionable cybersecurity policy, we now introduce an optimal control framework. This allows us to determine the most effective and cost-efficient strategies to contain a virus outbreak by dynamically allocating limited G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 17 of 29 resources. 3.1. Model with Time-Dependent Controls We introduce three time-dependent control variables into the SLAQRD model. The control variable u1(t) is the effort invested in preventing infection via removable media. This encompasses strategies like disabling auto-run features, enforcing device encryption, and user awareness training. We assume this control reduces the transmission rate γ by a factor of (1 − u1(t)), where 0 ≤ u1(t) ≤ umax 1 ≤ 1, u2(t) is the effort invested in en- hancing quarantine measures for actively infected computers. This represents improving intrusion detection systems (IDS) and automated isolation protocols. We assume this con- trol increases the quarantine rate θ to θ + u2(t), where 0 ≤ u2(t) ≤ umax 2 . and last u3(t) is the effort invested in enhancing recovery efforts for both active and quarantined ma- chines. This includes deploying more effective antivirus patches and streamlining cleanup processes. We assume this control increases the recovery rates δ and ξ to δ + u3(t) and ξ + u3(t) respectively, where 0 ≤ u3(t) ≤ umax 3 . These controls correspond to tangible cy- bersecurity measures. The media protection control u1(t) reflects policies such as enforcing USB encryption (e.g., using BitLocker), implementing device allow-listing software, and conducting mandatory user awareness training to prevent unauthorized removable media use. The quarantine control u2(t) represents the deployment of advanced Intrusion Detec- tion Systems (IDS) capable of automated threat isolation, such as network segmentation or endpoint isolation protocols. The recovery control u3(t) aligns with the rapid deploy- ment of antivirus patches, the use of automated remediation tools, and dedicated incident response teams for system cleaning and restoration. By framing the controls this way, the model offers a direct link between theoretical optimization and practical cybersecurity re- source allocation. Applying these controls, the controlled system of differential equations becomes: dS dt = Λ− β1SA− β2SL− (1− u1(t))γS − µS, dL dt = β1SA+ β2SL+ (1− u1(t))γS − σL− µL− αL, dA dt = σL− (θ + u2(t))A− (δ + u3(t))A− µA− ηA, dQ dt = (θ + u2(t))A− (ξ + u3(t))Q− µQ, dR dt = (δ + u3(t))A+ (ξ + u3(t))Q− µR, dD dt = αL+ ηA− µD. The objective is to minimize the total cost over a fixed time period [0, T ], which includes the cost of the infected populations (damage to the network) and the cost of implementing G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 18 of 29 the controls (resource expenditure). We define the objective functional J as: J(u1, u2, u3) = ∫ T 0 [ A1L(t) +A2A(t) +A3Q(t) + 1 2 ( B1u 2 1(t) +B2u 2 2(t) +B3u 2 3(t) )] dt. Here, A1, A2, A3 > 0 are weight constants balancing the cost associated with the latent, active, and quarantined populations, respectively. The quadratic terms 1 2Biu 2 i (t), with Bi > 0, reflect the nonlinear costs of implementing the controls, indicating that increasing effort becomes progressively more expensive. Our goal is to find an optimal control triple u∗ = (u∗1, u ∗ 2, u ∗ 3) such that: J(u∗) = min u∈U J(u1, u2, u3), where the control set U is defined as: U = {(u1, u2, u3) | ui(t) is Lebesgue measurable on [0, T ], 0 ≤ ui(t) ≤ umax i , i = 1, 2, 3} . 3.2. Existence of an Optimal Control The first step is to establish that such an optimal solution actually exists for our system. Theorem 9. There exists an optimal control triple u∗ = (u∗1, u ∗ 2, u ∗ 3) ∈ U and corresponding state solutions S∗, L∗, A∗, Q∗, R∗, D∗ that minimizes the objective functional J(u1, u2, u3) subject to the controlled system. Proof. The existence of an optimal control is established using standard results from Fleming and Rishel (1975) [29]. The following conditions must be verified: (i) The set of controls and corresponding state variables is non-empty. (ii) The control set U is closed and convex. (iii) The right-hand side of the state system is bounded by a linear function in the state and control variables. (iv) The integrand L of the objective functional, L(t, x, u) = A1L+A2A+A3Q+ 1 2(B1u 2 1+ B2u 2 2 +B3u 2 3), is convex on U . (v) There exist constants c1, c2 > 0 and ρ > 1 such that L(t, x, u) ≥ c1(|u1|2 + |u2|2 + |u3|2)ρ/2 − c2. Condition (1) is satisfied due to the boundedness and Lipschitz continuity of the state system’s right-hand side. Condition (2) is satisfied by the definition of U . The bound- edness of solutions (Theorem 2) ensures condition (3) holds. The integrand L is clearly convex in u as the Hessian matrix with respect to (u1, u2, u3) is a positive definite di- agonal matrix (B1, B2, B3 > 0). Finally, condition (5) is satisfied by choosing ρ = 2, c1 = 1 2 min{B1, B2, B3}, and a sufficiently large c2. Thus, an optimal control exists. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 19 of 29 3.3. Characterization of the Optimal Control Having established existence, we derive the necessary conditions that an optimal con- trol must satisfy using Pontryagin’s Maximum Principle. Theorem 10. Given an optimal control triple u∗ and corresponding state solutions x∗ = (S∗, L∗, A∗, Q∗, R∗, D∗), there exist adjoint variables λi(t), i = 1, . . . , 6 satisfying the fol- lowing adjoint system: dλ1 dt = −λ1(β1A ∗ + β2L ∗ + (1− u∗1)γ + µ) + λ2(β1A ∗ + β2L ∗ + (1− u∗1)γ), dλ2 dt = A1 − λ1β2S ∗ + λ2(β2S ∗ − (σ + µ+ α)) + λ3σ + λ6α, dλ3 dt = A2 − λ1β1S ∗ + λ2β1S ∗ −G+ λ4(θ + u∗2) + λ5(δ + u∗3) + λ6η, dλ4 dt = A3 − λ4((ξ + u∗3) + µ) + λ5(ξ + u∗3), dλ5 dt = −λ5µ, dλ6 dt = −λ6µ. with the transversality conditions (terminal conditions) at final time T , where G = λ3((θ+ u∗2) + (δ + u∗3) + µ+ η): λi(T ) = 0, for i = 1, 2, . . . , 6. Furthermore, the optimal controls u∗1, u∗2, and u∗3 are characterized by: u∗1(t) = min { umax 1 , max { 0, γS∗(λ2 − λ1) B1 }} , u∗2(t) = min { umax 2 , max { 0, A∗(λ3 − λ4) B2 }} , u∗3(t) = min { umax 3 , max { 0, A∗(λ3 − λ5) +Q∗(λ4 − λ5) B3 }} . Proof. Pontryagin’s Maximum Principle is applied by defining the Hamiltonian H: H = A1L+A2A+A3Q+ 1 2 (B1u 2 1 +B2u 2 2 +B3u 2 3) + 6∑ i=1 λifi, where fi are the right-hand sides of the controlled state equations. The adjoint system is derived from the conditions dλi dt = − ∂H ∂xi , where xi represents each state variable. The transversality conditions arise from the free terminal state. The characterization of the optimal controls is obtained by solving the optimality condition ∂H ∂uj = 0 for j = 1, 2, 3 at u∗ subject to the control constraints. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 20 of 29 3.4. Optimal Control Strategy and Comparative Effectiveness The efficacy of the proposed optimal control measures u1(t) (removable media protec- tion), u2(t) (enhanced quarantine), and u3(t) (enhanced recovery)—is evaluated through numerical simulation over a 100-day period. The system, initialized with S(0) = 800, L(0) = 50, A(0) = 30, Q(0) = 10, R(0) = 5, and D(0) = 2 computers, is simulated under two scenarios: a baseline case with no control and an intervention case with the optimal control strategy applied. The results, summarized in Table 2, demonstrate the profound impact of the intervention. Table 2: Comparative summary of key outcomes with and without the application of optimal control strategies. Metric Without Control With Optimal Control Peak Total Infections (L+A) 129.4 83.2 Final Recovered Systems (R) 192.5 312.8 Final Deactivated Systems (D) 58.3 31.6 % Reduction in Peak Infections – 35.7% % Increase in Final Recovered – 62.5% % Reduction in Final Deactivated – 45.8% Fig. (6) presents a comparative time-series analysis of the primary model compart- ments. The most critical observation is the stark contrast in the Total Infected population (L + A). The uncontrolled scenario (solid blue line) exhibits a sharp outbreak, reaching a peak of approximately 129 infected systems around day 20. In contrast, the controlled scenario (dashed red line) demonstrates a significantly mitigated outbreak, suppressing the peak infection by 35.7% to just 83 systems. This suppression is directly attributable to the application of controls u2(t) and u3(t), which increase the quarantine and recov- ery rates, respectively, effectively removing infectious individuals from the transmission chain more rapidly. Consequently, the controlled system maintains a higher population of Susceptible computers and achieves a substantially larger final population of Recovered systems (312.8 vs. 192.5), while minimizing the number of permanently Deactivated sys- tems, as quantified in Table 2. To further demonstrate the superiority of dynamic resource allocation, the optimal control strategy was compared against a static policy where control efforts were held constant at their average values over the simulation period (u1 = 0.4, u2 = 0.35, u3 = 0.5). While the static policy still outperformed the no-control scenario, the dynamic optimal strategy achieved a 18.3% further reduction in peak infections and a 22.1% higher final recovery rate. This comparison underscores that the time-varying, adaptive allocation of resources is crucial for maximizing containment efficiency and min- imizing operational damage, rather than merely applying uniform effort throughout an outbreak. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 21 of 29 0 10 20 30 40 50 60 70 80 90 100 Time (days) 200 300 400 500 600 700 800 N u m b e r o f C o m p u te rs Comparison of Susceptible Computers Without Control With Optimal Control 0 10 20 30 40 50 60 70 80 90 100 Time (days) 0 50 100 150 200 250 N u m b e r o f C o m p u te rs Comparison of Total Infected Computers (L + A) Without Control With Optimal Control 0 10 20 30 40 50 60 70 80 90 100 Time (days) 0 5 10 15 20 25 30 35 40 45 N u m b e r o f C o m p u te rs Comparison of Actively Infected Computers Without Control With Optimal Control 0 10 20 30 40 50 60 70 80 90 100 Time (days) 0 5 10 15 20 25 30 35 N u m b e r o f C o m p u te rs Comparison of Quarantined Computers Without Control With Optimal Control 0 10 20 30 40 50 60 70 80 90 100 Time (days) 0 50 100 150 200 250 300 350 400 N u m b e r o f C o m p u te rs Comparison of Recovered Computers Without Control With Optimal Control 0 10 20 30 40 50 60 70 80 90 100 Time (days) 0 50 100 150 200 250 300 N u m b e r o f C o m p u te rs Comparison of Deactivated Computers Without Control With Optimal Control Figure 6: Compartmental dynamics and control trajectories. Fig. (7) illustrates the temporal profiles of the three optimal control efforts. The strategy employs a high-intensity initial response: u1(t) (removable media protection) and u2(t) (enhanced quarantine) are initially applied at near-maximum capacity (approx. 0.81 and 0.72, respectively) to aggressively curb the initial exponential growth phase of the outbreak. After suppressing the peak (around day 30), the strategy adapts; u1(t) and u2(t) are gradually reduced to conserve resources, while u3(t) (enhanced recovery, G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 22 of 29 initially at 0.43) is ramped up to its maximum value of 0.85 to efficiently clean and return quarantined and infected systems to the recovered state. This time-dependent, adaptive allocation of resources exemplifies a cost-effective strategy that aligns with the theoretical solution from Pontryagin’s Maximum Principle. 0 10 20 30 40 50 60 70 80 90 100 Time (days) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 C o n tr o l E ff o rt Optimal Control Strategies u1: Removable Media Protection u2: Enhanced Quarantine u3: Enhanced Recovery Figure 7: Optimal Control Profiles Fig. (8) provides a bar chart analysis quantifying the percentage improvement achieved by the optimal control strategy. The results, derived from the data in Table 2, offer a clear and concise summary of the intervention’s success. The strategy yields a 35.7% reduction in the maximum number of simultaneously infected computers, directly translating to reduced operational disruption. Furthermore, it facilitates a 62.5% increase in the number of systems successfully recovered and returned to service, enhancing network resilience. Most importantly, it achieves a 45.8% reduction in the number of systems permanently deactivated due to severe infection, representing significant cost savings by avoiding the total loss of critical assets. This tri-fold improvement underscores the practical value of the proposed optimal control framework for cybersecurity resource allocation. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 23 of 29 Effectiveness of Optimal Control Strategies Peak Infections Final Recovered Final Deactivated 0 50 100 150 200 250 300 N u m b e r o f C o m p u te rs Without Control With Optimal Control Figure 8: Quantitative effectiveness assessment Fig. (9) synthesizes the key results of the optimal control simulation, quantifying the percentage improvement achieved by the strategy u∗ = (u∗1, u ∗ 2, u ∗ 3) over the uncontrolled baseline scenario. The bars demonstrate a significant reduction in the peak infection burden (∼34.5%) and final system deactivation (∼27.9%), alongside a substantial increase (∼5.2%) in successfully recovered systems. This tri-fold improvement underscores the practical value of the proposed dynamic resource allocation framework for minimizing the operational and financial damage of a cyber outbreak. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 24 of 29 Percentage Improvement with Optimal Control34.5% 5.2% 27.9% Reduction in Peak Infections Increase in Recovery Reduction in Deactivation 0 5 10 15 20 25 30 35 P e rc e n ta g e I m p ro v e m e n t (% ) Figure 9: Percentage Improvement with Optimal Control 4. Discussion This study has developed and analyzed a novel six-compartment SLAQRD model to investigate the propagation dynamics of computer viruses, with a specific emphasis on the critical and often overlooked role of removable media as a transmission vector. The model integrates realistic cybersecurity features, including latent infections, automated quarantine, system recovery, and the permanent deactivation of critically compromised machines. The mathematical analysis established the model’s well-posedness, derived the basic reproduction number R0 as an outbreak threshold, and characterized equilibrium states. Furthermore, an optimal control framework was formulated to identify dynamic, cost-effective intervention strategies. The numerical simulations yielded several pivotal insights. First, the phase portrait anal- ysis revealed complex system behaviors, including stable foci and saddle points, under- scoring the non-linear and sensitive nature of virus spread. This sensitivity highlights that minor changes in initial conditions or parameters can lead to drastically different outcomes, G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 25 of 29 reinforcing the need for proactive monitoring. Second, the simulations robustly validated the theoretical stability analysis; the system consistently converged to the disease-free equilibrium when R0 < 1 and to an endemic state when R0 > 1. The parameter sen- sitivity analysis, inherent in these simulations, clearly identified the quarantine rate (θ) and the infection rates (β1, β2) as the most influential parameters on R0. This provides a quantitative basis for cybersecurity policy, suggesting that investments in intrusion de- tection systems (to increase θ) and user training to reduce infection risk yield the highest return in containing outbreaks. The most significant practical contribution of this work is the formulation and evalua- tion of an optimal control strategy. The results demonstrate that a dynamic approach aggressively deploying quarantine (u2(t)) and media protection (u1(t)) efforts early in an outbreak, followed by a strategic shift towards enhanced recovery (u3(t)) is vastly superior to a static defense posture. This strategy achieved a 35.7% reduction in peak infections, a 62.5% increase in recovered systems, and a 45.8% reduction in irrecoverably deactivated assets. This tri-fold improvement translates directly into minimized operational disrup- tion, preserved network capacity, and significant cost avoidance, offering a clear blueprint for efficient resource allocation during a cyber incident. The theoretical findings of this study have direct implications for real-world cybersecurity policy and resource allocation. The model parameters, such as infection rates (β1, β2), quarantine rate (θ), and recovery rates (δ, ξ), can be estimated from organizational network data, including intrusion detec- tion system (IDS) logs, incident response reports, and malware analysis. For instance, the propagation dynamics observed during incidents involving worm-like malware (e.g., Con- ficker) or ransomware attacks leveraging removable media (e.g., past USB-based attacks) provide a contextual basis for parameter ranges. The optimal control strategy proposed— prioritizing early quarantine and media protection followed by recovery efforts aligns with established cybersecurity best practices, such as the NIST Cybersecurity Framework’s ”Respond” and “Recover” functions. This alignment demonstrates the model’s utility as a quantitative decision-support tool for security administrators to optimize investment in defensive measures under budget constraints. Despite its contributions, this work has limitations that present avenues for future research. The model assumes a homogeneous network, whereas real-world organizational networks have complex, hierarchical structures. This homogeneity assumption may lead to an underestimation of the potential outbreak scale and speed, as it fails to capture the accelerated propagation dynamics that occur through highly connected nodes (hubs) in scale-free network topologies common in real IT infrastructures. In such networks, a malware incursion targeting a central server or a highly connected user device could lead to a far more severe and rapid epidemic than predicted by our current model, potentially biasing the estimated effectiveness of control measures if not properly accounted for. Future work could extend this model onto complex network topologies (e.g., scale-free, small-world) to explore the impact of node connectiv- ity and critical hubs on propagation dynamics. Furthermore, the model parameters were held constant; integrating stochasticity to account for the unpredictable nature of user behavior and threat evolution would enhance its realism. Finally, a promising direction is the integration of machine learning for real-time parameter estimation and predictive G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 26 of 29 control, moving from a theoretical optimal control to an adaptive, learning-based defense system that can respond to evolving threats. 4.1. Implementation Feasibility and Challenges While the optimal control strategy demonstrates significant theoretical benefits, its real-world implementation faces several challenges. Scaling dynamic quarantine measures (u2(t)) requires robust intrusion detection systems (IDS) capable of real-time threat anal- ysis and automated isolation without causing operational disruption. Organizations may struggle with the initial cost of deploying such systems and the need for continuous mon- itoring. Similarly, scaling recovery efforts (u3(t)) demands efficient patch management tools and dedicated incident response teams, which may be resource-intensive for smaller organizations. Additionally, the effectiveness of media protection controls (u1(t)) hinges on employee compliance with security policies, such as restricting unauthorized USB device usage, which can be difficult to enforce consistently. These challenges highlight the need for balanced investment in both technology and human factors to successfully translate the proposed strategy into practice. 5. Conclusion This research has presented a comprehensive mathematical framework for understand- ing and combating the propagation of computer viruses that leverage both network-based and removable media vectors. The novel SLAQRD compartmental model developed in this work addresses a critical gap in cybersecurity literature by incorporating realistic states such as latent infections, automated quarantine, system recovery, and deactivation, thereby providing a more nuanced depiction of modern cyber threats. Theoretical analysis established the model’s mathematical robustness, defining a feasible region and proving the existence and uniqueness of solutions. The derivation of the basic reproduction num- ber R0 provided a crucial epidemiological threshold, determining the conditions under which an outbreak will fail or persist. Stability analysis, supported by numerical simula- tions, confirmed that the system converges predictably to a disease-free equilibrium when R0 < 1 or to an endemic state when R0 > 1, validating the model’s predictive power. The most significant contribution of this study is the formulation and application of an optimal control framework. By dynamically allocating resources to three key defense strategies removable media protection (u1), enhanced quarantine (u2), and accelerated recovery (u3) the model demonstrates a path to significantly mitigating cyber outbreaks. The results prove that a time varying, adaptive strategy is vastly superior to static de- fenses, achieving a substantial reduction in peak infections (35.7%), a major increase in recovered systems (62.5%), and a sharp decrease in irrecoverable system loss (45.8%). This translates directly into preserved operational continuity, maintained asset availability, and substantial cost avoidance for organizations. In conclusion, this work transitions cyberse- curity policy from a reactive to a predictive and optimized paradigm. It provides network administrators and security policymakers with a quantitative, decision-support tool for G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 27 of 29 resource allocation, emphasizing that investments in proactive detection, rapid isolation, and efficient recovery yield the highest returns in resilience. By translating the principles of mathematical epidemiology into actionable cybersecurity insights, this research offers a powerful foundation for building more defensible and resilient digital infrastructures in an increasingly threatened landscape. Acknowledgements The authors would like to acknowledge Prince Sultan University and EIAS Lab for their valuable support. Further, the authors would like to acknowledge Prince Sultan University for paying the Article Processing Charges (APC) of this publication. This paper is derived from a research grant “Cybersecurity Research and Innovation Pioneers Grants Initiative” funded by The National Program for RDI in Cybersecurity (National Cybersecurity Authority) Kingdom of Saudi Arabia - with grant number (CRPG-25-3168). The authors would like to acknowledge Prince Sultan University and EIAS Lab for their valuable support. Conflict of interest We have no conflict of interest regarding this work. Data Availability Statement All the data used in this article are included in the manuscript. No additional datasets were created or analyzed. AI Involvement Declaration The authors affirm that this research was conducted entirely through human intellec- tual effort, without the use of artificial intelligence tools during conceptualization, analysis, writing, or editing. References [1] P. Wang. Analysis of computer virus defense strategy based on network security. Academic Journal of Computing & Information Science, 5(14):33–39, 2022. [2] R. Anderson. Security Engineering: A Guide to Building Dependable Distributed Systems. Wiley Publishing, 2 edition, 2008. [3] P. Rai and V. Nain. Stuxnet unveiled: The blueprint for modern cyber conflict. In 2024 3rd Edition of IEEE Delhi Section Flagship Conference (DELCON), pages 1–4, 2024. [4] S. Mohurle and M. Patil. A brief study of wannacry threat: Ransomware attack 2017. International Journal, 8(5), 2017. G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 28 of 29 [5] A. Din and Y. Li. Optimizing hiv/aids dynamics: stochastic control strategies with education and treatment. The European Physical Journal Plus, 139(9):812, 2024. [6] J. O. Kephart and S. R. White. Measuring and modeling computer virus prevalence. In Proceedings of the 1993 IEEE Computer Society Symposium on Research in Security and Privacy, pages 2–15, 1993. [7] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani. Epidemic processes in complex networks. Reviews of Modern Physics, 87(3):925, 2015. [8] A. A. Ghorbani, W. Lu, and M. Tavallaee. Network Intrusion Detection and Pre- vention: Concepts and Techniques, volume 47 of Advances in Information Security. Springer Science & Business Media, 2009. [9] Z. Chen, L. Gao, and K. Kwiat. Modeling the spread of active worms. In Proceedings IEEE INFOCOM 2003, volume 3, pages 1890–1900, 2003. [10] H. J. Hadi, A. Khalid, F. B. Hussain, N. Ahmad, and M. A. Alshara. Flsh: A framework leveraging similarity hashing for android malware and variant detection. IEEE Access, 2025. [11] N. Ahmad, A. S. Rana, H. J. Hadi, F. B. Hussain, P. Chakrabarti, M. A. Alshara, and T. Chakrabarti. Geaad: generating evasive adversarial attacks against android malware defense. Scientific Reports, 15(1):11867, 2025. [12] A. Hawana, E. S. Hassan, W. El-Shafai, and S. A. El-Dolil. Enhancing malware detection with deep learning convolutional neural networks: Investigating the impact of image size variations. Security and Privacy, 8(2):e70000, 2025. [13] S. Staniford, V. Paxson, and N. Weaver. How to own the internet in your spare time. In Proceedings of the 11th USENIX Security Symposium, pages 149–167, 2002. [14] P. J. Taylor, T. Dargahi, A. Dehghantanha, R. M. Parizi, and K. K. R. Choo. A systematic literature review of blockchain cybersecurity. Digital Communications and Networks, 6(2):147–156, 2019. [15] M. Guri. Usbee: Air-gap covert-channel via electromagnetic emission from usb. In 2016 IEEE 35th Symposium on Reliable Distributed Systems (SRDS), pages 1–10, 2016. [16] L. X. Yang and X. Yang. The spread of computer viruses under the influence of removable storage devices. Applied Mathematics and Computation, 219(8):3914–3922, 2012. [17] M. T. Jafar, L. X. Yang, G. Li, and X. Yang. Optimal control of malware propagation in iot networks, 2024. arXiv:2401.11076. [18] S. Muthukumar, A. Balakumar, and V. Chinnadurai. The dynamics of the frac- tional seiqr malware spread model on wireless sensor networks. Journal of Analysis, 32(4):2349–2370, 2024. [19] S. Hosseini. Defense against malware propagation in complex heterogeneous networks. Cluster Computing, 24(2):1199–1215, 2021. [20] B. K. Mishra and S. K. Pandey. Dynamic model of worms with vertical transmission in computer network. Applied Mathematics and Computation, 217(21):8438–8446, 2018. [21] H. Yuan and G. Chen. Network virus-epidemic model with the point-to-group infor- G. Ali et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6967 29 of 29 mation propagation. Applied Mathematics and Computation, 219(2):515–522, 2011. [22] I. Shah, I. Alrabaiah, H. Alrabaiah, and B. Ozdemir. Using advanced analysis together with fractional order derivative to investigate a smoking tobacco cancer model. Results in Physics, 1(1):106700, 2023. [23] I. Shah, I. Ali, A. Ali, I. Ahmad, S. Islam, G. Rasool, S. Formanova, and M. Kallel. Optimal control and sensitivity analysis of a mathematical model for mdr-tb trans- mission with advanced treatment strategies. The European Physical Journal Plus, 140(6):1–15, 2025. [24] S. Ahmad, N. Ahmad, and I. Shah. Stability and sensitivity analysis of cyberattack propagation models in computer networks. European Journal of Pure and Applied Mathematics, 18(3):6336, 2025. [25] S. Ahmad, M. A. Elaffendi, N. Ahmad, and I. Shah. Machine learning-enhanced simulation of multi-vector email malware spread in organizational networks. European Journal of Pure and Applied Mathematics, 18(3):6542, 2025. [26] S. R. Chawla, S. Ahmad, W. Albalawi, A. Khan, I. Shah, and M. R. Eid. Stability analysis of a modified general seir model with harmonic mean type of incidence rate. Alexandria Engineering Journal, 2025. [27] F. A. Aboaoja, A. Zainal, et al. Malware detection issues, challenges, and future directions: A survey. Applied Sciences, 12(17):8482, 2022. [28] A. Sturaro, S. Silvestri, M. Conti, and S. K. Das. A realistic model for failure prop- agation in interdependent cyber-physical systems. IEEE Transactions on Network Science and Engineering, 7(2):817–831, 2018. [29] W. H. Fleming and R. W. Rishel. Deterministic and Stochastic Optimal Control, volume 1. Springer Science & Business Media, 2012.