EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6347 ISSN 1307-5543 – ejpam.com Published by New York Business Global Mathematical Modeling of SARS-CoV-2 Epidemics Using Fractional Calculus and Optimal Interventions Nadeem Abbas1, Wasfi Shatanawi1,2,∗, Syeda Alishwa Zanib3 1 Department of Mathematics and Sciences, College of Humanities and Sciences, Prince Sultan University, Riyadh, 11586, Saudi Arabia 2 Department of Mathematics, Faculty of Science, The Hashemite University, P.O Box 330127, Zarqa 13133, Jordan 3Department of Mathematics, Riphah International University, Main Satyana Road, Faisalabad 44000, Pakistan Abstract. In December 2019, the SARS-CoV-2 (COVID-19) virus was identified and quickly spread worldwide, causing a major global health crisis. To investigate its transmission dynam- ics, we developed a ten-compartment mathematical model, named CoVCom10, which includes key stages such as asymptomatic (F ), pre-symptomatic (E), and vaccinated (V ) individuals. The basic reproduction number (R0) has been calculated to evaluate how easily the virus can spread. We analyzed the local and global stability of the disease-free equilibrium and prove that the disease under control after vaccination when R0 < 1. A sensitivity analysis was conducted to assess the impact of key parameters, including the vaccination rate from susceptible individuals (β), trans- mission from susceptible to pre-symptomatic individuals (ϕ), and the rate of vaccination from pre-symptomatic individuals (γ). To evaluate intervention strategies, we extended the model by incorporating time-dependent control variables representing vaccination (a1), hospitalization (a2), and isolation of asymptomatic individuals (a3). The Pontryagin Maximum Principle was applied to identify optimal control strategies. Numerical simulations reveal that these interventions signif- icantly reduce virus transmission, particularly as the fractional-order parameter (ς) approaches 1, which aligns with observed real-world disease dynamics. The study emphasizes the effectiveness of integrated vaccination and treatment strategies in controlling the spread of COVID-19. 2020 Mathematics Subject Classifications: 26A33, 34A08, 03C65 Key Words and Phrases: SARS-CoV-2, Fractional-Order Model, Compartmental Model, Basic Reproduction Number, Stability Analysis, Pontryagin’s Maximum Principle ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6347 Email addresses: nabbas@psu.edu.sa (N. Abbas), wshatanawi@psu.edu.sa (W. Shatanawi), 19907@riphahfsd.edu.pk (S. A. Zanib) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 2 of 34 1. Introduction In the middle of the 1960s, human coronaviruses were first discovered. The well-known coronaviruses that may infect individuals include SARS-CoV (severe acute respiratory syndrome), MERS-CoV (Middle East Respiratory Syndrome), and SARS-CoV-2 (the new coronavirus that causes coronavirus disease 2019, or COVID-19), which is the subject of this research. Corona, virus, and disease are represented, respectively, by the letters CO, VI, and D in COVID-19. In December 2019, a novel virus was discovered during an outbreak in Wuhan, China [1–4] shown in Figure 1. There were attempts to control it, but they did not succeed, allowing the virus to spread to other regions of China and, subsequently the propagation of the world. According to a survey of individuals who passed away, the majority of them were elderly or had serious diseases like parkinson’s, lung, diabetes, chronic heart, or kidney disease. Flu and other viruses that can propagate contact with the mouth and touch with the nose can spread quickly. Coronaviruses are very dangerous and spread readily from person to person [5]. Figure 1: SARS-CoV-2 Structure Many mathematicians are constantly working to construct new, more effective models that may be used for modeling to evaluate the relationship between death and infection, fluid dynamic and predict how it will spread in the future [6–9]. Fractional calculus plays a vital role in biological modeling, offering a more accurate and flexible framework for captur- ing the memory and hereditary characteristics inherent in biological systems. Numerous researchers have developed mathematical models using fractional calculus to enhance the precision of numerical simulations and better reflect real-world disease dynamics [10]. In a short period, various types of research on COVID-19 model, pandemic have been con- ducted in the literature. Diagne et al., (2021) [11] studied formulation of COVID-19 model with vaccination. In their model, Various epidemiological stages were created for the entire population N based on each person’s health at any given moment t. They was shown how to regulate the model to minimize the spread of COVID-19 by employing Pontryagin’s maximal principle. Acheampong et al., (2022) [12] developed a model where N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 3 of 34 Pre-Symptomatic and infectious epidemiological classes generate infectious disease mod- els, leaving them as abstract ideas. Because of the prevalence of asymptomatic carriers, it was difficult to identify those who have been Pre-Symptomatic (represented by E) to or infected (represented by I) with SARS-CoV-2. They had produced two epidemiological classes: (1) a known group of Pre-Symptomatic individuals thought to be SARS-CoV-2 (represented by Q), and (2) those whose SARS-CoV-2 status has been clinically verified (represented by P ). Those who have been recognized as vulnerable were marked with the letter Q as they were required to quarantine under Ghana’s COVID-19 principles. The same was applied to confirm positive (P ) where clinical tests have shown that they had SARS-CoV-2. They named their model CoVCom-9. Butt et al. (2023) [13] developed a nonlinear SEIQHR fractional model for COVID-19 using the Atangana–Baleanu (ABC) derivative to capture the complex dynamics of disease transmission. They analyzed equi- librium points, performed sensitivity and bifurcation analyses, and applied an optimal con- trol framework with the Toufik–Atangana numerical method to assess control strategies. Zanib et al., (2024) [14] developed a modified compartmental COVID-19 model incorpo- rating vaccination strategies using conformable fractional derivatives to capture complex transmission dynamics. The basic reproduction number R0, its sensitivity indices, and the model’s stability were analyzed, with a finite difference method providing accurate numer- ical solutions and highlighting the role of vaccination in disease control. Butt et al., (2024) [15] developed a nonlinear fractional bi-susceptible S1S2V1V2IHR model using the Atan- gana–Baleanu Caputo derivative to investigate the dynamics and control of COVID-19. They analyzed the model’s stability, validated results using the Toufik–Atangana numer- ical method, and demonstrated the effectiveness of optimized control strategies through simulations. Abboubakar and Racke (2025) [16] developed a COVID-19 model using inte- ger and Caputo fractional-order derivatives, incorporating vaccination, confinement, and treatment with limited resources. Using German data, it was shown that the fractional model provided more accurate long-term forecasts, with a reproduction number around 1.90, indicating endemic persistence. Kumar et al., (2025)[17] developed an age-structured SEIR model to analyze the spread of COVID-19 and estimate key parameters such as the basic reproduction number R0 and case fatality ratio (CFR). The model, validated using epidemiological data and uncertainty analysis, provided insights for public health inter- ventions and was suggested to be extended using agent-based modeling. While reviewing the literature, we observed that many existing COVID-19 models skip critical compo- nents, particularly the inclusion of a vaccinated compartment. This limitation reduces their capacity to realistically capture disease dynamics in the post-vaccine era. To ad- dress this gap, we developed a comprehensive model that consist for all possible stages of COVID-19 progression. Specifically, we propose a ten-compartment mathematical model, termed CoVCom10, which incorporates key categories such as asymptomatic (F ), pre- symptomatic (E), and vaccinated (V ) individuals. CoVCom10 builds upon the previous CoVCom9 framework by incorporating the vaccination compartment, enabling a more accurate simulation of immunization effects on SARS-CoV-2 transmission. The addition of asymptomatic and pre-symptomatic compartments is essential for capturing hidden transmission routes and reflecting the full spectrum of COVID-19 progression. These N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 4 of 34 improvements enhance the epidemiological relevance of the model and support the evalua- tion of various intervention strategies. To further improve the model’s realism, we employ the conformable fractional derivative, which accounts for memory and hereditary effects. This approach enables the model to incorporate the influence of past states on current dynamics, thereby capturing complex behaviors such as delayed responses and long-term persistence of infections. The conformable fractional framework preserves key properties of classical calculus while offering improved computational tractability compared to other fractional definitions. Notably, it reduces to the classical model when the fractional or- der approaches one, ensuring consistency with traditional differential models. To validate the model’s applicability, we compared its simulation results with real-world COVID-19 data [18]. This comparison confirms the model’s capability to realistically represent the pandemic’s progression and evaluate the effectiveness of various public health measures. Despite the computational challenges associated with solving fractional differential equa- tions, the proposed model offers a valuable and flexible framework for analyzing the roles of vaccination, asymptomatic carriers, and control strategies in managing COVID-19 out- breaks. N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 5 of 34 2. Model Description and Formulation Table 1: List of Symbols, Parameters, and Abbreviations State Variables S(t) Number of susceptible individuals at time t E(t) Number of pre-symptomatic individuals at time t U(t) Number of infected individuals at time t Q(t) Number of quarantined individuals at time t P (t) Number of confirmed positive individuals at time t H(t) Number of hospitalized individuals in ordinary wards at time t C(t) Number of individuals in intensive care unit at time t F (t) Number of asymptomatic individuals at time t V (t) Number of vaccinated individuals at time t R(t) Number of recovered individuals at time t Model Parameters ϕ Transmission rate from S to E λ1 Transition rate from E to U λ2 Transition rate from E to Q γ Transition rate from E to V α1 Recovery rate of U α2 Transition rate from U to P b1 Transition rate from Q to S b2 Transition rate from Q to V b3 Transition rate from Q to P φ1 Transition rate from P to H φ2 Transition rate from P to C φ3 Transition rate from P to F m1 Recovery rate from H m2 Transition rate from H to C m3 Transition rate from H to F σ1 Self-isolation rate of F σ2 Transition rate from F to H η Recovery rate from C to H β Vaccination rate of S Λ Recruitment/birth rate into S µ Natural death rate d1–d7 Disease-induced death rates in E, U , Q, P , H, C, and F τ Loss of immunity rate from R to S Abbreviations DFEP Disease-Free Equilibrium Point SARS-CoV-2 Severe Acute Respiratory Syndrome Coronavirus 2 COVID-19 Coronavirus Disease 2019 CoVCom10 Coronavirus Compartment Model with 10 compartments CFD Conformable Fractional Derivative N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 6 of 34 The total population N(t) is stratified into ten distinct epidemiological compartments to comprehensively capture the transmission dynamics of COVID-19, as illustrated in the flowchart in Figure 2. The transitions between these compartments occur in continuous time and are governed by a system of nonlinear ordinary differential equations. Figure 2: Schematic diagram of the CoVCom10 compartmental model. The governing system for the model is system of differential equations shown below: dS dt = Λ+ τR+ b2Q− (ϕE + β + µ)S, dE dt = ϕES + γ V E − (λ1U + λ2Q+ µ+ d1)E, dU dt = λ1EU − (α1 + α2 + µ+ d2)U, dQ dt = λ2EQ− (b1 + b2 + b3 + µ+ d3)Q, dP dt = α2U + b1Q− (φ1 + φ2 + φ3 + µ+ d4)P, dH dt = φ1P + η C + σ2F − (m1 +m2 +m3 + µ+ d5)H, dC dt = φ2P +m2H − (η + µ+ d6)C, dF dt = φ3P +m3H − (σ1 + σ2 + µ+ d7)F, dV dt = β S + b3Q− (γ E + µ)V, dR dt = σ1F + α1U +m1H − (τ + µ)R. (2.1) N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 7 of 34 Therefore, N(t) = S(t) + E(t) + U(t) +Q(t) + P (t) +H(t) + C(t) + F (t) + V (t) +R(t). Each compartment, along with its associated abbreviation and biological interpretation, is detailed in Table 1, which also includes a complete list of model parameters and their descriptions. Fractional-Order Model Formulation Khalil et al. [19] introduced the conformable fractional derivative (CFD), a mathemati- cal operator that generalizes classical calculus while maintaining key differential properties. The CFD of order ζ ∈ (0, 1] is defined as: Dζ tM(t) = lim ϵ→0 M(t+ ϵt1−ζ)−M(t) ϵ , (2.2) which simplifies to the classical derivative when ζ = 1. This operator satisfies the compo- sition rule: Dζ tM(t) = t1−ζ dM dt . (2.3) To capture memory effects in COVID-19 transmission, we reformulate the CoVCom10 model (2.1) using the CFD framework [19]: Dζ t S = Λ+ τR+ b2Q− (ϕE + β + µ)S, Dζ tE = ϕES + γV E − (λ1U + λ2Q+ µ+ d1)E, Dζ tU = λ1EU − (α1 + α2 + µ+ d2)U, Dζ tQ = λ2EQ− (b1 + b2 + b3 + µ+ d3)Q, Dζ tP = α2U + b1Q− (φ1 + φ2 + φ3 + µ+ d4)P, Dζ tH = φ1P + ηC + σ2F − (m1 +m2 +m3 + µ+ d5)H, Dζ tC = φ2P +m2H − (η + µ+ d6)C, Dζ tF = φ3P +m3H − (σ1 + σ2 + µ+ d7)F, Dζ t V = βS + b3Q− (γE + µ)V, Dζ tR = σ1F + α1U +m1H − (τ + µ)R. (2.4) The system is solved under non-negative initial conditions: S(0) = S0 ≥ 0, E(0) = E0 ≥ 0, U(0) = U0 ≥ 0, Q(0) = Q0 ≥ 0, P (0) = P0 ≥ 0, H(0) = H0 ≥ 0, C(0) = C0 ≥ 0, F (0) = F0 ≥ 0, V (0) = V0 ≥ 0, R(0) = R0 ≥ 0. (2.5) N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 8 of 34 Mathematical Analysis The biologically feasible region for the fractional-order model (2.4) is defined as: D = { (S(t), E(t), U(t), Q(t), P (t), H(t), C(t), F (t), V (t), R(t)) ∈ R10 + : S + E + U +Q+ P +H + C + F + V +R ≤ Λ µ } , (2.6) where D forms a positive invariant set that captures all epidemiologically meaningful states of the system. The mathematical analysis of the model demonstrate the the feasible region D in terms of bounded and non-negative, stability results, and epidemiological thresholds. Solutions of equations (2.4) with initial conditions (2.5) keep in D discuss in the below theorem. Theorem 1. (Positive Invariance of Feasible Region) For the fractional-order sys- tem (2.4) with non-negative initial conditions (2.5), the closed set D defined in (2.6) forms a positively invariant region under conformable fractional dynamics when all parameters satisfy {ϕ, β, µ, . . .} ∈ R+. Proof. Let N(t) = ∑10 i=1Xi(t) represent the total population, where Xi denotes each compartment. Applying the conformable fractional derivative operator: Dζ tN(t) = Λ− µN(t)− 7∑ i=1 diXi(t). (2.7) Using the conformable derivative property Dζ tN = t1−ζ dN dt , we rewrite: t1−ζ dN dt = Λ− µN− 7∑ i=1 diXi(t)︸ ︷︷ ︸ ≥0 . (2.8) This establishes the inequality: dN dt ≤ tζ−1(Λ− µN). (2.9) Solving this fractional differential inequality through separation of variables:∫ N(t) N(0) dN Λ− µN ≤ ∫ t 0 τ ζ−1dτ, (2.10) taking integrated factor and after simplify: N(t) ≤ Λ µ − ( Λ µ − N(0) ) exp ( −µ ζ tζ ) . (2.11) For N(0) ≤ Λ µ , the exponential term ensures N(t) ≤ Λ µ ∀t ≥ 0. Thus, all solutions remain bounded within D, making it positively invariant under the fractional-order dynamics. N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 9 of 34 Positivity and Boundedness of Solutions Theorem 2. For the fractional-order system (2.4) with non-negative initial conditions (2.5) and non-negative parameters, the following holds: • All solutions {S(t), E(t), U(t), Q(t), P (t), H(t), C(t), F (t), V (t), R(t)} remain non- negative for all t ≥ 0. • The total population satisfies lim sup t→∞ N(t) ≤ Λ µ . Proof. Positivity Consider the conformable fractional derivative formulation Dζ tX = t1−ζX ′. For each compartment X ∈ {S,E,U,Q, P,H,C, F, V,R}: Lemma 1. Let Dζ tX ≥ −ψX with ψ ≥ 0 and X(0) ≥ 0. Then X(t) ≥ 0 ∀t ≥ 0. For the susceptible population: Dζ t S = Λ+ τR+ b2Q− (ϕE + β + µ)S ≥ −(ϕE + β + µ)S, (2.12) Applying the comparison principle for fractional differential equations: S(t) ≥ S(0) exp ( −1 ζ (ϕE + β + µ)tζ ) ≥ 0. (2.13) Similar analysis for other compartments yields: Dζ tE ≥ −(λ1U + λ2Q+ µ+ d1)E, Dζ tU ≥ −(α1 + α2 + µ+ d2)U, ... (2.14) By sequential application of the comparison lemma, all compartments maintain non- negativity. Boundedness From Theorem 1 (Positive Invariance of Feasible Region), the total population dynamics satisfy: Dζ tN = Λ− µN− 7∑ i=1 diXi ≤ Λ− µN. (2.15) Solving the fractional inequality: N(t) ≤ Λ µ − ( Λ µ − N(0) ) exp ( −µ ζ tζ ) . (2.16) As t→ ∞, the exponential term vanishes, yielding: lim sup t→∞ N(t) ≤ Λ µ . (2.17) N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 10 of 34 2.1. Disease-free Equilibrium Point The disease-free equilibrium point (DFE) [20] of the model is determined by setting all compartments associated with infection to zero, reflecting the absence of disease in the population. Specifically, we set E = U = Q = P = H = C = F = R = 0. At equilibrium, the rates of change for the susceptible (S) and vaccinated (V ) compart- ments are given by the following steady-state equations: dS dt = Λ− (β + µ)S = 0, dV dt = βS − µV = 0. Solving these equations yields the equilibrium values: S∗ = Λ β + µ , V ∗ = Λβ (β + µ)µ . (2.18) Therefore, the disease-free equilibrium point is E0 = (S,E,U,Q, P,H,C, F, V,R) = ( Λ β + µ , 0, 0, 0, 0, 0, 0, 0, Λβ (β + µ)µ , 0 ) . 3. Basic Reproduction Number To analyze the transmission potential of the CoVCom10 model (2.4), we employ the next-generation matrix approach developed by Van den Driessche and Watmough [21]. This method involves decomposing the system into two main components: the transmis- sion matrix, which describes the generation of new infections, and the transition matrix, which captures the movement of individuals between different epidemiological states. The transmission matrix A(x∗) and the transition matrix B(x∗) are constructed as follows: A(x∗) =  γ V E + ϕES λ1EU λ2EQ 0 0 0 0  , (3.19) N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 11 of 34 B(x∗) =  −Π1E −Π2U −Π3Q α2U + b1Q−Π4P η C + σ2F −Π5H + φ1P −Π6C +m2H + φ2P −Π7F +m3H + φ3P  . (3.20) The next-generation matrix AB−1 evaluated at the disease-free equilibrium yields, AB−1 =  1 Π1 ( ϕΛ β+µ + β Λ γ (β+µ)µ ) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  . (3.21) Here, the composite parameters Πi are defined to combine various transition rates: Π1 = λ1U + λ2Q+ µ+ d1, Π2 = α1 + α2 + µ+ d2, Π3 = b1 + b2 + b3 + µ+ d3, Π4 = φ1 + φ2 + φ3 + µ+ d4, Π5 = m1 +m2 +m3 + µ+ d5, Π6 = η + µ+ d6, Π7 = σ1 + σ2 + µ+ d7. (3.22) By evaluating the Jacobians of A(x∗) and B(x∗) at the disease-free equilibrium (DFE) and computing the spectral radius of AB−1, we obtain the basic reproduction number: R0 = Λ (γ β + µϕ) µ (β + µ) (µ+ d1) . (3.23) 3.1. Sensitivity Analysis Sensitivity analysis quantifies how variations in model parameters affect the basic reproduction number R0 [22]. The normalized sensitivity index, defined as ΥR0 ξ , represents N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 12 of 34 the relative change in R0 resulting from a relative change in a parameter ξ, mathematically given by: ΥR0 ξ = ∂R0 ∂ξ × ξ R0 . (3.24) Applying this definition to model reproduction number: R0 = Λ (γ β + µϕ) µ (β + µ) (µ+ d1) , (3.25) we compute sensitivity indices for key parameters: • Sensitivity with respect to β: ΥR0 β = γ β γ β + µϕ − β β + µ . (3.26) An increase in β will increase R0, indicating a higher transmission potential. • Sensitivity with respect to γ: ΥR0 γ = γ β γ β + µϕ . (3.27) A higher γ elevates R0, highlighting increased transmission from asymptomatic or pre-symptomatic individuals. • Sensitivity with respect to ϕ: ΥR0 ϕ = µϕ γ β + µϕ . (3.28) An increase in exposure rate ϕ raises R0, demonstrating greater susceptibility within the population. • Sensitivity with respect to µ: ΥR0 µ = µϕ γ β + µϕ − µ β + µ − µ µ+ d1 − 1. (3.29) Increasing the natural death rate µ generally reduces R0, as fewer individuals remain susceptible to infection. • Sensitivity with respect to d1: ΥR0 d1 = − d1 µ+ d1 . (3.30) Increasing disease-induced death rate d1 lowers R0, due to a reduction in the number of infectious contacts. N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 13 of 34 • Sensitivity with respect to Λ: ΥR0 Λ = 1. (3.31) An increase in Λ which is birth rate will increase R0, indicating a higher transmission potential. This sensitivity analysis identifies critical parameters influencing disease dynamics, guiding effective strategies for controlling the epidemic shown in Figure 3. d1 Sensitive Parameters 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 Se ns iti vi ty In di ce s 0.013 0.996 0.004 -1.036 -0.957 1.000 Figure 3: Sensitivity analysis 3.2. Stability Analysis Local Stability of Disease-Free Equilibrium Point To assess the behavior of the CoVCom10 model (2.4) near the disease-free equilibrium, it is essential to analyze its local stability. This analysis determines whether small distur- bances from the disease-free state will decay or grow, which is crucial for understanding epidemic control strategies. The following theorem and proof outline the conditions under which the disease-free equilibrium point (DFEP) is locally asymptotically stable. Theorem 3. The disease-free equilibrium point E0 of the CoVCom10 model (2.4) is locally asymptotically stable if the basic reproduction number R0 < 1 and unstable if R0 > 1 [23]. Proof. To establish the local stability of the DFEP, we compute the Jacobian matrix of the system at E0. The Jacobian matrix J0 is given by: N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 14 of 34 J0 =  −β − µ − ϕΛ β+µ 0 b2 0 0 0 0 0 τ 0 γβΛ (β+µ)µ + ϕΛ β+µ −Π1 0 0 0 0 0 0 0 0 0 0 −Π2 0 0 0 0 0 0 0 0 0 0 −Π3 0 0 0 0 0 0 0 0 α2 b1 −Π4 0 0 0 0 0 0 0 0 0 φ1 −Π5 η σ2 0 0 0 0 0 0 φ2 m2 −Π6 0 0 0 0 0 0 0 φ3 m3 0 −Π7 0 0 β − γβΛ (β+µ)µ 0 b3 0 0 0 0 −µ 0 0 0 α1 0 0 m1 0 σ1 0 −τ − µ  (3.32) The characteristic equation of J0 is obtained from det(J0−∆I) = 0, where I is the identity matrix. The coefficients of the characteristic polynomial are: c1 = 1, c2 = Π5 +Π6 +Π7, c3 = (Π5 +Π7)Π6 − ηm2 −m3σ2 +Π5Π7, c4 = (Π5Π7 −m3σ2)Π6 − ηΠ7m2. (3.33) The characteristic polynomial can be written as: 1 µ(β + µ) [ (β + µ+∆)(∆µβ +∆µ2 − Λβγ − Λµϕ+ βµΠ1 + µ2Π1) × (Π2 +∆)(Π3 +∆)2(Π4 +∆)(µ+∆)(τ + µ+∆)2 ×(β + µ+∆)(c1∆ 3 + c2∆ 2 + c3∆+ c4) ] = 0 (3.34) The first seven roots of the characteristic equation are: ∆1 = −(β + µ), ∆2 = (R0 − 1) ( 1 µ(β + µ) ) , ∆3 = −Π4, ∆4 = −µ, ∆5 = −(τ + µ), ∆6 = −Π3, ∆7 = −(τ + µ). (3.35) All explicit eigenvalues have negative real parts when R0 < 1. For the cubic polynomial c1∆ 3 + c2∆ 2 + c3∆+ c4 = 0, the Routh-Hurwitz stability conditions are: c1 > 0, c2 > 0, c1c2c3 > c23 + c21c4. (3.36) N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 15 of 34 These conditions are satisfied when R0 < 1, as confirmed by the structure of the composite parameters Πi. Therefore, all eigenvalues of the Jacobian matrix J0 have negative real parts if and only if R0 < 1, ensuring local asymptotic stability of the DFEP. Conversely, if R0 > 1, the DFEP becomes unstable, indicating the potential for an epidemic outbreak. Global Stability of Disease-Free Equilibrium Point Lemma 2. (Castillo-Chavez Method [24]) The DFEP E0 = (X0,0) of system (2.4) is globally asymptotically stable if: (Z1) For dX dt = F (X, 0), X0 is globally asymptotically stable (Z2) H(X,Y) = PY− Ĥ(X,Y) satisfies Ĥ(X,Y) ≥ 0 in Ω, where: – 1. P = DYH(X0, 0) is Metzler (non-negative off-diagonal elements) – 2. Ω is the biologically feasible region Theorem 4. The CoVCom10 model (2.4) is globally asymptotically stable at DFEP E0 when R0 < 1, satisfying both Castillo-Chavez conditions [24]. Proof. Firstly, to satisfy condition (Z1), the model (2.4) are rewrite by setting, TH = (S, V ) and, GH = (E,U,Q, P,H,C, F,R). Then, disease-free equilibrium point is given by the fixed point, E0 = ( X0, 0 ) = ( Λ β + µ , β Λ µ (β + µ) ) , the system dTH dt = F (TH , 0) becomes, dS∗ dt = Λ− (β + µ)S, dV ∗ dt = β S − (µ)V. (3.37) By solving Eq. (3.37), the equation has a unique equilibrium point, (S∗, V ∗) = ( Λ β + µ , β Λ µ (β + µ) ) , (3.38) hence X0 is globally asymptotically stable. So we can say the condition (Z1) is fulfilled. Now, to satisfy the second condition (Z2), H(TH , GH) = PHGN − Ĥ(TH , GH), and N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 16 of 34 Ĥ(TH , GH) ≥ 0, For that, system of equations (2.4). We have, H(TH , GH) =  EϕS + γ V E − (λ1U + λ2Q+ µ+ d1)E λ1EU − (α1 + α2 + µ+ d2)U λ2EQ− (b1 + b2 + b3 + µ+ d3)Q α2U + b1Q− (φ1 + φ2 + φ3 + µ+ d4)P η C + σ2F + φ1P − (m1 +m2 +m3 + µ+ d5)H m2H + φ2P − (η + µ+ d6)C φ3F +m3H − (σ1 + σ2 + µ+ d7)F σ1F +m1H + α1U − (τ + µ)R  , (3.39) Ĥ(TH , GH) = PHGN −H(TH , GH) =  ((ϕ(S∗ − S) + γ(V ∗ − V )− λ1U − λ2Q)E λ1EU λ2EQ 0 0 0 0 0  , (3.40) this shows that, Ĥ(TH , GH) ≥ 0, where GN represent an M matrix, it contains a non- negative off-diagonal element. Therefore, Both conditions (Z1) and (Z2) are satisfied when R0 < 1, proving global asymptotic stability of DFEP E0 by Lemma 2. 3.3. Existence and Uniqueness of Solution The mathematical results presented in this section that model (2.4) and (2.5) (CoV- Com10) have a unique solution under certain reasonable assumptions. From an epidemio- logical viewpoint, the existence of solutions implies that the model reliably predicts disease dynamics, confirming that realistic initial conditions and parameters will always yield a meaningful trajectory of the disease. Uniqueness assures that the models predictions are consistent and reproducible, crucial for decision-making in public health. The CoVCom10 model’s solutions, which are provided in system of equation (2.4) and (2.5) are described in this section by their qualitative qualities. The following Volterra-type integral equation results from first taking the both sides integral, where ∫ ζ t is the integration function having N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 17 of 34 the order ζ with respect to t Now by using the definition of Khalilzadeh [19, 20], we get, S(t)− S(0) = ∫ t 0 ρζ−1 [Λ + τR(ρ) + b2Q(ρ)− (ϕE(ρ) + β + µ)S(ρ)] dρ, E(t)− E(0) = ∫ t 0 ρζ−1 [ϕE(ρ)S(ρ) + γ V (ρ)E(ρ)− (λ1U(ρ) + λ2Q(ρ) + µ+ d1)E(ρ)] dρ, U(t)− U(0) = ∫ t 0 ρζ−1 [λ1E(ρ)U(ρ)− (α1 + α2 + µ+ d2)U(ρ)] dρ, Q(t)−Q(0) = ∫ t 0 ρζ−1 [λ2E(ρ)Q(ρ)− (b1 + b2 + b3 + µ+ d3)Q(ρ)] dρ, P (t)− P (0) = ∫ t 0 ρζ−1 [α2I(ρ) + b1Q(ρ)− (φ1 + φ2 + φ3 + µ+ d4)P (ρ)] dρ, H(t)−H(0) = ∫ t 0 ρζ−1 [φ1P (ρ) + η C(ρ) + σ2F (ρ)− (m1 +m2 +m3 + µ+ d5)H(ρ)] dρ, C(t)− C(0) = ∫ t 0 ρζ−1 [φ2P (ρ) +m2H(ρ)− (η + µ+ d6)C(ρ)] dρ, F (t)− F (0) = ∫ t 0 ρζ−1 [φ3P (ρ) +m3H(ρ)− (σ1 + σ2 + µ+ d7)F (ρ)] dρ, V (t)− V (0) = ∫ t 0 ρζ−1 [β S(ρ) + b3Q(ρ)− (γ E(ρ) + µ)V (ρ)] dρ, R(t)−R(0) = ∫ t 0 ρζ−1 [σ1F (ρ) + α1I(ρ) +m1H(ρ)− (τ + µ)R(ρ)] dρ, (3.41) define the kernels in following, Φ1(t, S) = Λ + τR(t) + b2Q(t)− (ϕE(t) + β + µ)S(t), Φ2(t, E) = ϕE(t)S(t) + γ V (t)E(t)− (λ1U(t) + λ2Q(t) + µ+ d1)E(t), Φ3(t, U) = λ1E(t)U(t)− (α1 + α2 + µ+ d2)U(t), Φ4(t, Q) = λ2E(t)Q(t)− (b1 + b2 + b3 + µ+ d3)Q(t), Φ5(t, P ) = α2U(t) + b1Q(t)− (φ1 + φ2 + φ3 + µ+ d4)P (t), Φ6(t,H) = φ1P (t) + η C(t) + σ2F (t)− (m1 +m2 +m3 + µ+ d5)H(t), Φ7(t, C) = φ2P (t) +m2H(t)− (η + µ+ d6)C(t), Φ8(t, F ) = φ3P (t) +m3H(t)− (σ1 + σ2 + µ+ d7)F (t), Φ9(t, V ) = β S(t) + b3Q(t)− (γ E(t) + µ)V (t), Φ10(t, R) = σ1F (t) + α1U(t) +m1H(t)− (τ + µ)R(t).,m (3.42) Theorem 5. (Lipschitz Continuity and Contraction Mapping of Kernel Operators) Let Φi : R+ × X → X (i = 1, . . . , 10) be kernel operators defined on a Banach space X with norm ∥ · ∥. Assume: N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 18 of 34 • Population compartments are bounded: ∥ S ∥≤ r1, ∥ E ∥≤ r2, ∥ U ∥≤ r3, ∥ Q ∥≤ r4, ∥ P ∥≤ r5, ∥ H ∥≤ r6, ∥ C ∥≤ r7, ∥ F ∥≤ r8, ∥ V ∥≤ r9, ∥ R ∥≤ r10 • Lipschitz coefficients satisfy: s∗1 = ϕr2 + β + µ s∗2 = ϕr1 + γr9 +Π1 s∗3 = λ1r2 +Π2 s∗4 = λ2r2 +Π3 s∗5 = Π4, s ∗ 6 = Π5, s ∗ 7 = Π6, s ∗ 8 = Π7 s∗9 = γr2 + µ s∗10 = τ + µ • Parameter constraint: 0 ≤ s∗i < 1 ∀i ∈ {1, . . . , 10} Then each Φi satisfies the Lipschitz condition and constitutes a contraction mapping. Proof. We demonstrate the result for Φ1; analogous arguments apply to Φ2, . . . ,Φ10. Let S1, S2 ∈ X be arbitrary functions. Then: ∥ Φ1(t, S1)− Φ1(t, S2) ∥ =∥ [ Λ + τR(t) + b2Q(t)− (ϕE(t) + β + µ)S1(t) ] − [ Λ + τR(t) + b2Q(t)− (ϕE(t) + β + µ)S2(t) ] ∥ =∥ −(ϕE(t) + β + µ)(S1(t)− S2(t)) ∥ ≤ (ϕ ∥ E ∥ +β + µ) ∥ S1 − S2 ∥ (by triangle inequality) ≤ s∗1 ∥ S1 − S2 ∥, (3.43) where s∗1 = ϕr2 + β + µ by the boundedness assumption ∥ E ∥≤ r2. The contraction property follows from 0 ≤ s∗1 < 1. Similar calculations for Φ2, . . . ,Φ10 yield corresponding Lipschitz constants s∗2, . . . , s ∗ 10 with contraction properties under the stated parameter con- straints. Therefore, all kernel operators satisfy both Lipschitz continuity and contraction mapping requirements. By considering the kernels Φi, i = 1, 2, 3, . . . , 10. Now the system of equation have (3.41) then, Recursive formula can be proceed in the following, 01v = S(t)− S(0) = ∫ t 0 ρζ−1 (Φ1(ρ, Sv−1)− Φ1(ρ, Sv−2)) dρ, (3.44) triangle inequality law will be apply on Eq. (3.44), ∥ 01v ∥=∥ Sv(t)− Sv−1(t) ∥≤ s∗1 ∫ t 0 ρζ−1 ∥ (Sv−1 − Sv−2)) ∥ dρ, (3.45) N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 19 of 34 by applying Lipschitz conditions (5), ∥ 01v ∥≤ s∗1 ∫ t 0 ∥ 01v−1 ∥ dρ, (3.46) it can write, Sv(t) = v∑ i=1 01v(t). (3.47) Similarly for others. Thus the following theorem may be derived from these findings. Theorem 6. (Existence of CoVCom10 Model Solutions) Let X = C([0, tmax],R+) be the Banach space of continuous functions with norm ∥ · ∥. The CoVCom10 model admits a unique solution (S,E,U,Q, P,H,C, F, V,R) ∈ X 10 if: • Lipschitz constants s∗i from Theorem 1 satisfy: sitmax ≤ 1. ∀i ∈ {1, . . . , 10} where tmax > 0 is the maximal existence time • Initial conditions (S0, . . . , R0) are bounded in X Proof. Using the Picard iteration scheme for the integral equations, define successive approximations: Sv+1(t) = S0 + ∫ t 0 Φ1(ρ, Sv(ρ))dρ, with analogous definitions for other compartments. From Theorem 1’s Lipschitz condi- tions: ∥ 01v ∥ ≤∥ S0 ∥ (s∗1tmax) v, ≤∥ S0 ∥ .(s∗1tmax) v. (by induction hypothesis) (3.48) For residual terms ℶ1v(t) := S(t)− Sv(t): ∥ ℶ1v(t) ∥ ≤ ∫ t 0 ∥ Φ1(ρ, S)− Φ1(ρ, Sv−1) ∥ dρ, ≤ s∗1 ∫ t 0 ∥ S − Sv−1 ∥ dρ, ≤ (s∗1tmax) v ∥ S0 ∥ . (via recursive estimation) Under condition s∗i tmax < 1, we get: lim v→∞ ∥ ℶ1v(t) ∥≤ lim v→∞ (s∗1tmax) v ∥ S0 ∥= 0. Similar convergence holds for other compartments by identical reasoning. By Banach fixed-point theorem, this establishes existence and uniqueness of solutions in X 10. N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 20 of 34 Theorem 7. (Uniqueness of Solutions for Fractional CoVCom10 Model) Let ζ ∈ (0, 1] be the fractional order and t ∈ [0, T ] with T < min{1/s∗i }10i=1. If the kernels Φi satisfy: • Lipschitz continuity: ∃s∗i > 0 such that ∥ Φi(t, x)− Φi(t, y) ∥≤ s∗i ∥ x− y ∥ . ∀x, y ∈ L1([0, T ]) • Time-domain constraint: 1− s∗i t ζ ≥ 0. ∀i ∈ {1, . . . , 10}, t ∈ [0, T ] then the fractional CoVCom10 system admits a unique solution (S,E,U,Q, P,H,C, F, V,R) ∈ (C[0, T ])10. Proof. Assume two distinct solutions X = (S, . . . , R) and X∗ = (S∗, . . . , R∗) exist. For the S-compartment: ∥ S(t)− S∗(t) ∥ ≤ ∫ t 0 ρζ−1 ∥ Φ1(ρ, S)− Φ1(ρ, S ∗) ∥ dρ, ≤ s∗1 ∫ t 0 ρζ−1 ∥ S(ρ)− S∗(ρ) ∥ dρ. (by Lipschitz condition) (3.49) Applying the generalized Gronwall inequality for fractional integrals: ∥ S(t)− S∗(t) ∥≤∥ S0 − S∗ 0 ∥ Eζ(s ∗ 1Γ(ζ)t ζ), where Eζ is the Mittag-Leffler function. Given identical initial conditions S0 = S∗ 0 and T < (s∗1) −1/ζ , the growth estimate implies: ∥ S(t)− S∗(t) ∥≤ 0 · Eζ(· · · ) = 0. ∀t ∈ [0, T ] Thus S(t) ≡ S∗(t). Repeating this steps for E(t), . . . , R(t) using their respective Lipschitz constants s∗2, . . . , s ∗ 10 completes the proof. 3.4. Optimal Control To stop COVID-19 from spreading, the effects will be examined by using medicinal treatments. To do this, a set of time-dependent control variables, a1, a2, and a3 have been introduced, • The implementation of continuous vaccination is represented by a1, • The social distancing and lockdown measures is represented by a2, • Testing and quarantine strategy is represented by a3. N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 21 of 34 The COVID-19 model with proposed optimal control a1, a2, and a3 are part of the following nonautonomous system of nonlinear ordinary differential equations. S′ = Λ+ τR+ b2Q− (ϕE + a1 + µ)S, (3.50) E′ = ϕES + γ V E − (λ1U + λ2Q+ µ+ d1)E, (3.51) U ′ = λ1EU − (α1 + α2 + µ+ d2)U, (3.52) Q′ = λ2EQ− (b1 + b2 + b3 + µ+ d3)Q, (3.53) P ′ = α2U + b1Q− (φ1 + φ2 + φ3 + µ+ d4)P, (3.54) H ′ = φ1P + η C + σ2F − (a2 +m2 +m3 + µ+ d5)H, (3.55) C ′ = φ2P +m2H − (η + µ+ d6)C, (3.56) F ′ = φ3P +m3H − (a3 + σ2 + µ+ d7)F, (3.57) V ′ = a1 S + b3Q− (γ E + µ+ a)V, (3.58) R′ = a3F + α1U + a2H − (τ + µ)R. (3.59) The discussion of optimal control in the model is done to order to determine the optimum values of a1,a2, and a3 that minimise the objective function J(a1(t), a2(t), a3(t)) affected by the differential equations (3.50-3.59). The provided objective function is, J(a1, a2, a3) = min a1,a2,a3 ∫ T 0 [X1E +X2U +X3Q+X4P +X5F +X6C +X7H +Y1 (a1 (t)) 2 + Y2 (a2 (t)) 2 + Y3 (a3 (t)) 2]dt. (3.60) Where T is final time,X1, X2, X3, X4, X5, X6, andX7 are the weight cost of Pre-Symptomatic, Infected, Quarantine, Confirm-positive, Asymptomatic, Hospitalized at intensive care unit and Hospitalized at ordinary ward individuals respectively. Where Y1, Y2 and Y3 are weight costs for each control measurable individuals. In this paper, the weight cost of the optimal control as applied by [25, 26] is measured using a quadratic function that fulfills the optimality criteria. The objective is to identify the optimal control (a∗1, a ∗ 2, a ∗ 3). J(a∗1, a ∗ 2, a ∗ 3) = min(a1, a2, a3). Where, (a1, a2, a3) ∈ U, ai(t) is measurable lebesgue on [0, T ], 0 ≤ ai(t) ≤ 1, i = 1, 2, 3. Hamiltonian and Optimality Equation Pontryangin’s Maximum Principle [27] is used to determine the requirements that an optimal control must satisfy. Equations (3.50) and (3.60) are transformed into a problem N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 22 of 34 of minimising the point-wise Hamiltonian (H) with respect to a1(t), a2(t), and a3(t) as a result of his principle. H = X1E +X2U +X3Q+X4P +X5H +X6C +X7F + Y1 (a1 (t)) 2 + Y2 (a2 (t)) 2 + Y3 (a3 (t)) 2 + ξ1[Λ + τR+ b2Q− (ϕE + a1 + µ)S] + ξ2[ϕES + γ V E − (λ1U + λ2Q+ µ+ d1)E] + ξ3[λ1EU − (α1 + α2 + µ+ d2)U ] + ξ4[λ2EQ− (b1 + b2 + b3 + µ+ d3)Q] + ξ5[α2U + b1Q− (φ1 + φ2 + φ3 + µ+ d4)P ] + ξ6[φ1P + η C + σ2F − (a2 +m2 +m3 + µ+ d5)H] + ξ7[φ2P +m2H − (η + µ+ d6)C] + ξ8[φ3P +m3H − (a3 + σ2 + µ+ d7)F ] + ξ9[a1 S + b3Q− (γ E + µ+ a)V ] + ξ10[a3F + α1U + a2H − (τ + µ)R], (3.61) where (ξi), i = 1, 2, ...., 10 are adjoint variables associate with S,E,U,Q, P, H,C, F, V and, R. Theorem 8. For the optimal control (a∗2, a ∗ 2, a ∗ 3) and corresponding state solution (S,E,U,Q, P,H,C, F, V,R) that minimize J over U of the corresponding system of equation (2.4) having the adjoint variable ξ1, ...., ξ10 such that, dξ1 dt = (ξ1 − ξ2)ϕE + (ξ1 − ξ9)a1 − aξ1, (3.62) dξ2 dt = (ξ1 − ξ2)ϕS + (ξ9 − ξ2)γV + (ξ2 − ξ3)λ1U − (ξ2 − ξ4)λ2Q+ ξ2(µ+ d1)−X1, (3.63) dξ3 dt = (ξ2 − ξ3)λ1E + (ξ3 − ξ10)α1 + (ξ3 − ξ5)α2 + ξ3(µ+ d2)−X2, (3.64) dξ4 dt = (ξ2 − ξ4)λ2E + (ξ4 − ξ5)b1 + (ξ4 − ξ1)b2 + (ξ4 − ξ9)b3 + ξ4(µ+ d3)−X3, (3.65) dξ5 dt = (ξ5 − ξ6)φ1 + (ξ5 − ξ7)φ2 + (ξ5 − ξ8)φ3 + ξ5(µ+ d4)−X4, (3.66) dξ6 dt = (ξ6 − ξ10)a2 + (ξ6 − ξ7)m2 + (ξ6 − ξ8)m3 + ξ6(µ+ d5)−X5, (3.67) dξ7 dt = (ξ7 − ξ6)η + ξ7(µ+ d6)−X6, (3.68) dξ8 dt = (ξ8 − ξ10)a3 + (ξ8 − ξ6)σ2 + ξ8(µ+ d7)−X7, (3.69) dξ9 dt = (ξ9 − ξ2)γE + ξ9(µ)− ξ9a, (3.70) dξ10 dt = ξ10 (µ+ τ)− ξ1τ. (3.71) ξi(T ) = 0, for i = 1, 2, 3....10 having conditions, a∗1 = max{0,min(1,(ξ1 − ξ9)S 2Y1 )}, a∗2 = max{0,min(1, (ξ6 − ξ10)H 2Y2 )}, a∗3 = max{0,min(1, (ξ8 − ξ10)F 2Y3 )}, (3.72) N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 23 of 34 Proof. Using Pontryagin’s maximal principle [28], adjoint system is written as a result of our differentiation of the Hamiltonian with repect to the various states S,E,U,Q, P,H,C, F, V and R. dξ1 dt = −dH dS , dξ2 dt = −dH dE , dξ3 dt = −dH dI , dξ4 dt = −dH dQ , dξ5 dt = −dH dP , dξ6 dt = − dH dH , dξ7 dt = −dH dC , dξ8 dt = −dH dF , dξ9 dt = −dH dV , dξ10 dt = −dH dR . (3.73) With conditions, ξi(T ) = 0, for i = 1, 2, 3....10. The optimal controls (a∗1, a ∗ 2, a ∗ 3) are characterized by adopting the strategy used by Pon- tryagin et al., [27], based on the following conditions, ∂H ∂ai , for i=1,2,3....10. for a∗i , the result are, a∗1 = (ξ1 − ξ9)S 2Y1 , a∗2 = (ξ6 − ξ10)H 2Y2 , a∗3 = (ξ8 − ξ10)F 2Y3 , (3.74) By including the described control set, initial, and transversal conditions, the optimal control system is generated from the adjoint variable system and the optimal control system. 4. Numerical Results and Discussion Figure 4 presents a comparison between the model-predicted number of infected in- dividuals (U) and the smoothed daily new COVID-19 cases per 100,000 people reported in the United States from January 2022 onward. The model output is shown in red, while the real-world data is plotted in green. The results demonstrate the alignment of model dynamics with observed trends in case data, highlighting the utility of the proposed compartmental framework for capturing the progression of COVID-19 infections A numerical simulation of the non-linear integer-order system (2.4) has been performed using the classical fourth-order Runge-Kutta (RK4) method. The initial conditions for all compartments are given by: S(0) = 500, E(0) = 20, U(0) = 10, Q(0) = 8, P (0) = 6, H(0) = 4, F (0) = 2, C(0) = 1, V (0) = 0, R(0) = 0. The system of ordinary differential equations governing the CoVCom10 model was solved using the fourth-order Runge-Kutta (RK4) method, which provides improved accuracy over basic Euler methods through its multi-stage slope calculations. For a general ODE system: dy dt = f(t,y), y(t0) = y0 where y = [S,E,U,Q, P,H, F,C, V,R]T , the RK4 method computes each time step as: k1 = f(tn,yn) N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 24 of 34 2022-01 2022-02 2022-03 2022-04 2022-05 2022-06 2022-07 2022-08 2022-09 2022-10 2022-11 2022-12 2023-01 2023-02 2023-03 2023-04 2023-05 2023-06 Date 0 2 4 6 8 10 M od el : I nf ec te d 0 50 100 150 200 Re al C as es (p er 1 00 k) COVID-19 Model vs Real Cases (USA) Model Infected (U) Real Cases (per 100k) Figure 4: Comparison of model-predicted infected population (U) with real-world COVID-19 cases per 100,000 people in the United States. The model output is shown in red, and the real-world data is plotted in green. k2 = f ( tn + h 2 ,yn + h 2 k1 ) k3 = f ( tn + h 2 ,yn + h 2 k2 ) k4 = f (tn + h,yn + hk3) with the state update given by: yn+1 = yn + h 6 (k1 + 2k2 + 2k3 + k4) where h represents the step size controlling temporal resolution. The method’s fourth- order accuracy makes it particularly suitable for capturing non-linear epidemiological dy- namics while maintaining numerical stability. N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 25 of 34 0 20 40 60 80 100 Days (Time) 500 1000 1500 2000 2500 3000 In di vi du al s Susceptible (S) 0 20 40 60 80 100 Days (Time) 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 In di vi du al s Pre-Symptomatic (E) 0 20 40 60 80 100 Days (Time) 4 5 6 7 8 9 10 In di vi du al s Infected (U) 0 20 40 60 80 100 Days (Time) 10 15 20 25 30 In di vi du al s Quarantined (Q) 0 20 40 60 80 100 Days (Time) 2 3 4 5 6 In di vi du al s Clinically Positive (P) 0 20 40 60 80 100 Days (Time) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 In di vi du al s Hospitalized (H) 0 20 40 60 80 100 Days (Time) 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 In di vi du al s Asymptomatic (F) 0 20 40 60 80 100 Days (Time) 0.5 0.6 0.7 0.8 0.9 1.0 In di vi du al s ICU (C) 0 20 40 60 80 100 Days (Time) 0 5 10 15 20 25 30 In di vi du al s Vaccinated (V) 0 20 40 60 80 100 Days (Time) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 In di vi du al s Recovered (R) Figure 5: Time series behavior of COVID-19 model using RK4 method The Figure 5 highlight the dynamic evolution of the different compartments in the CoV- Com10 model. The susceptible population increases over time, possibly due to natural recruitment or effective preventive measures reducing disease transmission. The exposed and pre-symptomatic populations decrease, indicating early identification and isolation of cases, which helps prevent further spread. The infected and clinically positive individuals show a gradual decline, suggesting that the disease is being controlled through diagnosis, treatment, and isolation. The quarantined population initially rises and then falls, reflect- ing efficient identification and management of at-risk individuals. A continuous decrease in hospitalized and ICU cases demonstrates effective healthcare response, while the asymp- tomatic population also declines, possibly due to detection and reclassification into other compartments. The vaccinated population shows a consistent increase, playing a critical role in protecting susceptible individuals and limiting transmission. The recovered popu- lation increases initially and then stabilizes, indicating successful treatment and acquired immunity among the affected individuals. Overall, these results emphasize the impor- N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 26 of 34 tance of timely public health interventions—such as vaccination, quarantine, treatment, and healthcare support—in curbing the spread of infection and improving population-level health outcomes. 0 20 40 60 80 100 Days 500 1000 1500 2000 2500 3000 3500 In di vi du al s Susceptible (S) Small Scale Large Scale 0 20 40 60 80 100 Days 0 10 20 30 40 50 60 70 80 In di vi du al s Pre-Symptomatic (E) Small Scale Large Scale 0 20 40 60 80 100 Days 10 20 30 40 50 60 In di vi du al s Infected (U) Small Scale Large Scale 0 20 40 60 80 100 Days 20 40 60 80 100 120 In di vi du al s Quarantined (Q) Small Scale Large Scale 0 20 40 60 80 100 Days 0 5 10 15 20 25 30 In di vi du al s Clinically Positive (P) Small Scale Large Scale 0 20 40 60 80 100 Days 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 In di vi du al s Hospitalized (H) Small Scale Large Scale 0 20 40 60 80 100 Days 0 2 4 6 8 10 In di vi du al s ICU (C) Small Scale Large Scale 0 20 40 60 80 100 Days 0 2 4 6 8 10 In di vi du al s Asymptomatic (F) Small Scale Large Scale 0 20 40 60 80 100 Days 0 5 10 15 20 25 30 35 In di vi du al s Vaccinated (V) Small Scale Large Scale 0 20 40 60 80 100 Days 0 2 4 6 8 10 12 14 In di vi du al s Recovered (R) Small Scale Large Scale Figure 6: Comparison of model dynamics under small- and large-scale population settings. To evaluate the robustness and reliability of the proposed CoVCom10 model, we con- ducted simulations using two distinct sets of initial population values. In the small-scale scenario, the initial conditions were: S(0) = 500, E(0) = 20, U(0) = 10, Q(0) = 8, P (0) = 6, H(0) = 4, F (0) = 2, C(0) = 1, V (0) = 0, R(0) = 0. In contrast, the large-scale scenario used proportionally higher values: S(0) = 1000, E(0) = 80, U(0) = 60, Q(0) = 50, P (0) = 30, H(0) = 20, F (0) = 10, C(0) = 10, V (0) = 0, R(0) = 0. The comparative outcomes, illustrated in Figure 6, demonstrate that the qualitative behavior of all compartments is consistent across both population scales. While the absolute values increase in line with population size, the temporal progression, peak values, and recovery dynamics remain N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 27 of 34 largely unchanged. This consistency confirms that the model robustly captures the intrin- sic transmission dynamics and is not overly sensitive to variations in initial population sizes. Peak E Final R Time to Peak E 0.2 0.4 0.6 0.8 1.0 0.00.00.0 Effect of Varying = 5.0e-05 = 1.0e-04 = 2.0e-04 Peak E Final R Time to Peak E 0.2 0.4 0.6 0.8 1.0 0.00.00.0 Effect of Varying = 1.0e-04 = 2.0e-04 = 4.0e-04 Peak E Final R Time to Peak E 0.2 0.4 0.6 0.8 1.0 0.00.00.0 Effect of Varying = 5.0e-05 = 1.0e-04 = 2.0e-04 Peak E Final R Time to Peak E 0.2 0.4 0.6 0.8 1.0 0.00.00.0 Effect of Varying = 3.0e-04 = 5.0e-04 = 7.0e-04 Figure 7: Radar plots showing the impact of varying ϕ, β, γ, and τ on peak exposed, final recovered, and time to peak exposure. Figure 7 presents radar plots illustrating the impact of four parameters ϕ, β, γ, and τ on key epidemic outcomes. Each subplot compares peak exposed (E), final recovered (R), and time to peak E across three values of the selected parameter. Higher ϕ and τ intensify the outbreak by increasing peak E and reducing final R. Conversely, increasing β and γ improves recovery and reduces the epidemic size. Time to peak remains relatively stable under most changes. These results emphasize the role of transmission control and immunity in outbreak management. 4.1. Impact of Optimal Control Strategies on Epidemic Disease Dynamics This section presents a detailed analysis of the impact of various optimal control strate- gies on COVID-19 transmission dynamics using the CoVCom10 model. We systematically evaluate the effectiveness of different combinations of interventions, such as vaccination, hospitalization, and management of asymptomatic cases—by applying Pontryagin’s Max- imum Principle to determine the optimal allocation and timing of these controls. Figure 8 illustrates the temporal evolution of all compartments under scenarios with and without optimal interventions. The simulation results clearly demonstrate that the coordinated im- plementation of vaccination, enhanced hospitalization capacity, and targeted management N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 28 of 34 of asymptomatic individuals leads to a substantial reduction in the number of infections and severe cases over time. Notably, the model quantifies the optimal intensity and timing for each control measure, providing actionable guidance for public health authorities to maximize resource efficiency and intervention impact. These findings underscore the criti- cal importance of integrated and sustained public health strategies. By combining multiple interventions, the spread of SARS-CoV-2 can be curbed more effectively and sustainably than by relying on any single measure alone. The numerical results offer practical insights for policymakers, highlighting the value of adaptive, data-driven approaches to epidemic control. 0 20 40 60 80 100 Days 500 1000 1500 2000 2500 3000 3500 Nu m be r o f S us ce pt ib le (S ) Susceptible (S) Dynamics Without Control With Control 0 20 40 60 80 100 Days 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Nu m be r o f P re -S ym pt om at ic (E ) Pre-Symptomatic (E) Dynamics Without Control With Control 0 20 40 60 80 100 Days 2 4 6 8 10 Nu m be r o f I nf ec te d (U ) Infected (U) Dynamics Without Control With Control 0 20 40 60 80 100 Days 10 15 20 25 30 Nu m be r o f Q ua ra nt in ed (Q ) Quarantined (Q) Dynamics Without Control With Control 0 20 40 60 80 100 Days 2 3 4 5 6 Nu m be r o f C lin ica lly P os iti ve (P ) Clinically Positive (P) Dynamics Without Control With Control 0 20 40 60 80 100 Days 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Nu m be r o f H os pi ta liz ed (H ) Hospitalized (H) Dynamics Without Control With Control 0 20 40 60 80 100 Days 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 Nu m be r o f A sy m pt om at ic (F ) Asymptomatic (F) Dynamics Without Control With Control 0 20 40 60 80 100 Days 0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0 Nu m be r o f I CU (C ) ICU (C) Dynamics Without Control With Control 0 20 40 60 80 100 Days 0 100 200 300 400 500 600 Nu m be r o f V ac cin at ed (V ) Vaccinated (V) Dynamics Without Control With Control 0 20 40 60 80 100 Days 0 1 2 3 4 5 6 7 Nu m be r o f R ec ov er ed (R ) Recovered (R) Dynamics Without Control With Control Figure 8: Temporal dynamics of COVID-19 compartments with and without optimal control inter- ventions. 4.2. Convergence Analysis of Conformable Fractional Derivative Model The structural convergence of the conformable fractional system (2.4) was rigorously analyzed through numerical simulations and quantitative error metrics. As demonstrated N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 29 of 34 in Figure 9, solution trajectories exhibit continuous dependence on the fractional order parameter ς, asymptotically approaching the classical integer-order model (2.1) as ς → 1. The convergence was quantified using the L2-norm difference: ∥yς(t)− y1(t)∥2 = √√√√ 10∑ i=1 ∫ tf t0 (yς,i(t)− y1,i(t))2dt where yς(t) = [Sς , Eς , . . . , Rς ] T and y1(t) represent the fractional and integer-order so- lutions respectively. Numerical solutions were computed using an implicit Runge-Kutta method with adaptive step size control h ∈ [0.01, 0.1] to maintain relative error tolerance below 10−5. 0 20 40 60 80 100 Days (Time) 500 1000 1500 2000 2500 Su sc ep tib le (S ) =0.1 =0.3 =0.5 =0.7 =0.9 0 20 40 60 80 100 Days (Time) 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Pr e- Sy m pt om at ic (E ) =0.1 =0.3 =0.5 =0.7 =0.9 0 20 40 60 80 100 Days (Time) 6 7 8 9 10 In fe ct ed (U ) =0.1 =0.3 =0.5 =0.7 =0.9 0 20 40 60 80 100 Days (Time) 10 15 20 25 30 Qu ar an tin ed (Q ) =0.1 =0.3 =0.5 =0.7 =0.9 0 20 40 60 80 100 Days (Time) 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Cl in ica lly P os iti ve (P ) =0.1 =0.3 =0.5 =0.7 =0.9 0 20 40 60 80 100 Days (Time) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Ho sp ita liz ed (H ) =0.1 =0.3 =0.5 =0.7 =0.9 0 20 40 60 80 100 Days (Time) 1.5 1.6 1.7 1.8 1.9 2.0 As ym pt om at ic (F ) =0.1 =0.3 =0.5 =0.7 =0.9 0 20 40 60 80 100 Days (Time) 0.70 0.75 0.80 0.85 0.90 0.95 1.00 IC U (C ) =0.1 =0.3 =0.5 =0.7 =0.9 0 20 40 60 80 100 Days (Time) 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 Va cc in at ed (V ) =0.1 =0.3 =0.5 =0.7 =0.9 0 20 40 60 80 100 Days (Time) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Re co ve re d (R ) =0.1 =0.3 =0.5 =0.7 =0.9 Figure 9: Dynamic convergence of conformable fractional model solutions to classical solutions as ς → 1. Shading represents solution variance across ς ∈ [0.1, 0.9]. N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 30 of 34 The conformable fractional derivative introduces time-dependent memory effects into epidemiological modeling, capturing how historical transmission dynamics, immune re- sponses, and behavioral adaptations continuously influence disease spread. Unlike classical integer-order models that assume instantaneous state transitions, the conformable oper- ator Dς ty = t1−ς dy dt incorporates a temporal scaling factor t1−ς , modulating transmission rates through the fractional order ς. For ς < 1, slower infection progression emerges, mim- icking real-world delays in interventions such as lockdowns or vaccine rollouts. This formu- lation avoids non-physical infinite memory artifacts while preserving short-term memory effects, such as waning immunity. By varying ς ∈ (0, 1], the model spans rapid containment (ς → 1) to prolonged transmission (ς ≪ 1), reflecting heterogeneous public health scenar- ios. This approach bridges idealized compartmental frameworks with empirical outbreak patterns, where cumulative interactions dictate epidemic trajectories. Table 2: Quantitative convergence analysis of fractional-order solutions ς L2-Norm Difference Normalized L2-Norm 0.1 14986.57 1.0000 0.3 14317.17 0.9553 0.5 12602.19 0.8409 0.7 9295.42 0.6202 0.9 3712.78 0.2477 1.0 0.00 0.0000 Table 3: Parameters used for the CoVCom10 Model Parameters Value Source Parameters Value Source Λ 50 [12] λ1 0.008 [12] µ 0.009 [12] λ2 0.0011 estimated γ 0.08 estimated α1 0.01004 estimated b1 0.002 [12] α2 0.0083 estimated b2 0.0013 estimated b3 0.00379 estimated φ1 0.001971 [12] m1 0.008619 estimated φ2 0.0018 estimated m2 0.0012 estimated φ3 0.004711 [12] m3 0.0003 estimated σ1 0.0018 estimated η 0.0009 estimated σ2 0.0021 estimated β 0.5 estimated β 0.5 [11] τ 0.05 [11] ϕ 0.02 estimated d1 0.26 [11] d2 0.25 estimated d5 0.22 estimated d3 0.24 estimated d6 0.21 estimated d4 0.23 estimated d7 0.20 estimated N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 31 of 34 5. Conclusion This study investigated the impact of vaccination on virus transmission and proposed effective control strategies using a ten-compartment mathematical model named CoV- Com10, formulated as a system of nonlinear differential equations. The model explicitly includes vaccination, asymptomatic, and pre-symptomatic carriers to better capture the real-world transmission dynamics. Our results highlight the critical role of vaccination in mitigating SARS-CoV-2 spread. The model alters the classical susceptible population expression from S0 = Λ µ (as seen in SEIR models) to S0 = Λ β+µ in CoVCom10, reflecting immunization’s role in reducing susceptibility. The basic reproduction number R0 was computed to quantify disease potential. Stability analysis confirms that the disease-free equilibrium is locally and globally asymptotically stable when R0 < 1, and the endemic equilibrium is stable when R0 > 1. Time-dependent control variables for vaccination (a1), hospitalization (a2), and asymptomatic isolation (a3) were introduced via the Pontryagin Maximum Principle. Simulations demonstrated that optimizing these controls signifi- cantly reduces transmission. Model parameters such as ϕ, β, and γ reflect real-world biological mechanisms of infection, recovery, and exposure. The success of these strategies depends on coordinated policy implementation and community behavior. Moreover, the modularity of CoVCom10 makes it adaptable to other infectious diseases like influenza or monkeypox, or emerging SARS-CoV-2 variants. Adjusting compartments and parameters allows the model to simulate various epidemiological scenarios. The use of the conformable fractional derivative enhances the model’s ability to capture memory effects and persis- tent behavior. We recommend that vaccination campaigns be complemented by personal preventive practices such as mask-wearing, hand hygiene, and social distancing. These combined interventions are essential for managing current and future infectious disease outbreaks. As a limitation of the model, despite its validation against real epidemiologi- cal data, it assumes homogeneous mixing and does not account for spatial heterogeneity or stochastic variations. Additionally, due to limited data availability, some parameters were estimated, which may affect the model’s precision. Future work could enhance this model by incorporating stochastic processes, spatial heterogeneity, with co-infection or network-based transmission to provide more realistic and robust policy guidance. Acknowledgements Authors (Dr. Nadeem Abbas and Prof. Dr. Wasfi Shatanawi) would like to thank Prince Sultan University for their support through the TAS research lab. Data Availability No data availability Conflict of Interest The authors declare that there is no conflict of interest. N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 32 of 34 References [1] Qun Li, Xuhua Guan, Peng Wu, Xiaoye Wang, Lei Zhou, Yeqing Tong, Ruiqi Ren, Kathy SM Leung, Eric HY Lau, Jessica Y Wong, et al. Early transmission dynamics in wuhan, china, of novel coronavirus–infected pneumonia. New England journal of medicine, 382(13):1199–1207, 2020. [2] Sk Shahid Nadim and Joydev Chattopadhyay. Occurrence of backward bifurcation and prediction of disease transmission with imperfect lockdown: A case study on covid-19. Chaos, Solitons & Fractals, 140:110163, 2020. [3] Joseph T Wu, Kathy Leung, and Gabriel M Leung. Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in wuhan, china: a modelling study. The lancet, 395(10225):689–697, 2020. [4] Sansao A Pedro, Frank T Ndjomatchoua, Peter Jentsch, Jean M Tchuenche, Madhur Anand, and Chris T Bauch. Conditions for a second wave of covid-19 due to interac- tions between disease dynamics and social processes. Frontiers in Physics, 8:574514, 2020. [5] Li-Xiang Feng, Shuang-Lin Jing, Shi-Ke Hu, De-Fen Wang, Hai-Feng Huo, et al. Modelling the effects of media coverage and quarantine on the covid-19 infections in the uk. Math. Biosci. Eng, 17(4):3618–3636, 2020. [6] Muhammad Riaz, Kamal Shah, Thabet Abdeljawad, Inas Amacha, Asma Al-Jaser, and Manar Alqudah. A comprehensive analysis of covid-19 nonlinear mathematical model by incorporating the environment and social distancing. Scientific Reports, 14(1):12238, 2024. [7] Eyup Cetin, Serap Kiremitci, and Baris Kiremitci. Developing optimal policies to fight pandemics and covid-19 combat in the united states. European Journal of Pure and Applied Mathematics, 13(2):369–389, 2020. [8] Yasir Nawaz, Muhammad Shoaib Arif, and Kamaleldin Abodayeh. Predictor– corrector scheme for electrical magnetohydrodynamic (mhd) casson nanofluid flow: a computational study. Applied Sciences, 13(2):1209, 2023. [9] Yasir Nawaz, Muhammad Shoaib Arif, Kamaleldin Abodayeh, Muhammad Usman Ashraf, and Mehvish Naz. A new explicit numerical scheme for enhancement of heat transfer in sakiadis flow of micro polar fluid using electric field. Heliyon, 9(10), 2023. [10] Muhammad Farhan, Zahir Shah, Rashid Jan, and Saeed Islam. A fractional modeling approach of buruli ulcer in possum mammals. Physica Scripta, 98(6):065219, 2023. [11] ML Diagne, H Rwezaura, SY Tchoumi, and JM Tchuenche. A mathematical model of covid-19 with vaccination and treatment. Computational and Mathematical Methods in Medicine, 2021(1):1250129, 2021. [12] Edward Acheampong, Eric Okyere, Samuel Iddi, Joseph HK Bonney, Joshua Kiddy K Asamoah, Jonathan ADWattis, and Rachel L Gomes. Mathematical modelling of ear- lier stages of covid-19 transmission dynamics in ghana. Results in Physics, 34:105193, 2022. [13] AIK Butt, W Ahmad, M Rafiq, N Ahmad, and M Imran. Optimally analyzed fractional coronavirus model with atangana–baleanu derivative. Results in Physics, N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 33 of 34 53:106929, 2023. [14] Syeda Alishwa Zanib, Tamour Zubair, Sehrish Ramzan, Muhammad Bilal Riaz, Muhammad Imran Asjad, and Taseer Muhammad. A conformable fractional finite difference method for modified mathematical modeling of sar-cov-2 (covid-19) disease. PloS one, 19(10):e0307707, 2024. [15] Azhar Iqbal Kashif Butt, Waheed Ahmad, Hafiz Ghulam Rabbani, Muhammad Rafiq, Shehbaz Ahmad, Naeed Ahmad, and Saira Malik. Exploring optimal control strategies in a nonlinear fractional bi-susceptible model for covid-19 dynamics using atangana- baleanu derivative. Scientific Reports, 14(1):31617, 2024. [16] Hamadjam Abboubakar and Reinhard Racke. Mathematical modeling of the coro- navirus (covid-19) transmission dynamics using classical and fractional derivatives. 2022. [17] Komal Sharma, Vivek Srivastava, and Ravi Kant Singh. From data to cures: lever- aging machine learning, deep learning and pharmacore modelling for targeted ther- apies. In AIP Conference Proceedings, volume 3254, page 020008. AIP Publishing LLC, 2025. [18] Our World in Data. Covid-19 dataset. https://covid.ourworldindata.org/data/owid- covid-data.csv. [19] Roshdi Khalil, Mohammed Al Horani, Abdelrahman Yousef, and Mohammad Sabab- heh. A new definition of fractional derivative. Journal of computational and applied mathematics, 264:65–70, 2014. [20] Nadeem Abbas, Syeda Alishwa Zanib, Sehrish Ramzan, Aqsa Nazir, and Wasfi Shatanawi. A conformable mathematical model of ebola virus disease and its sta- bility analysis. Heliyon, 10(16), 2024. [21] Pauline Van den Driessche and James Watmough. Reproduction numbers and sub- threshold endemic equilibria for compartmental models of disease transmission. Math- ematical biosciences, 180(1-2):29–48, 2002. [22] Helena Sofia Rodrigues, M Teresa T Monteiro, and Delfim FM Torres. Seasonality ef- fects on dengue: basic reproduction number, sensitivity analysis and optimal control. Mathematical Methods in the Applied Sciences, 39(16):4671–4679, 2016. [23] SM Ashrafur Rahman. Study of infectious diseases by mathematical models: predic- tions and controls. 2016. [24] Carlos Castillo-Chavez, Sally Blower, Pauline van den Driessche, Denise Kirschner, and Abdul-Aziz Yakubu. Mathematical approaches for emerging and reemerging infec- tious diseases: models, methods, and theory, volume 126. Springer Science & Business Media, 2002. [25] Getachew Teshome Tilahun, Oluwole Daniel Makinde, and David Malonza. Modelling and optimal control of typhoid fever disease with cost-effective strategies. Computa- tional and mathematical methods in medicine, 2017(1):2324518, 2017. [26] Getachew Teshome Tilahun, Oluwole Daniel Makinde, and David Malonza. Co- dynamics of pneumonia and typhoid fever diseases with cost effective optimal control analysis. Applied Mathematics and Computation, 316:438–459, 2018. [27] Lev Semenovich Pontryagin. Mathematical theory of optimal processes. Routledge, N. Abbas, W. Shatanawi, S. A. Zanib / Eur. J. Pure Appl. Math, 18 (3) (2025), 6347 34 of 34 2018. [28] K Renee Fister, Suzanne Lenhart, and Joseph Scott McNally. Optimizing chemother- apy in an hiv model. 1998.