EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6001 ISSN 1307-5543 – ejpam.com Published by New York Business Global Mathematical Modeling with Two Strains to Investigate Evolution of COVID-19 from Computational and Theoretical Perspectives Muhammad Riaz1, Zeeshan Ali2,∗, Rahim Ud Din1, Sadique Ahmad3, Kamaleldin Abodayeh4, Muhammad Sarwar1,4 1Department of Mathematics, University of Malakand, Chakdara Dir(L), 18000, Khyber Pakhtunkhwa, Pakistan 2Department of Information Management, National Yunlin University of Science and Technology Douliu Taiwan, R. O. C 3EIAS Data Science and Block Chain Laboratory, College of Computer and Information Sciences, Prince Sultan University, Riyadh 11586, Saudi Arabia 4Department of Mathematics and Sciences, Prince Sultan University, P.O. Box 66833, 11586 Riyadh, Saudi Arabia Abstract. The respiratory disease COVID-19 is brought on by a mutagenic ribonucleic acid (RNA) virus. Variants with various traits that could potentially impact transmissibility started to appear globally in December 2020. We develop and examine a computational model of two-strain COVID-19 transmission behaviors in order to deal with this novel aspect of the disease. After a theoretical analysis, enough criteria are derived for the stability of the model’s equilibrium. The model contains single-strain 1 and variant 2 endemic equilibria as well as disease-free and endemic equilibria. The global stability of the model equilibrium is demonstrated using the Lyapunov function. For each strain, we separately compute the basic reproductive numbers R01 < 1, and R02 < 1. To determine the realistic values of the model parameters, the actual data from South Africa for the period of three months, are taken into consideration. Sensitivity study employing the partial rank correlation coefficient technique (PRCC) to look into the key variables that affect R01 < 1, or R02 < 1 decrease or increase. We show the numerical simulation of the model using non-standard finite difference scheme (NSFDS) and Caputo fractional derivative for distinct non- integer order. The world health organizations (WHO) recommends standard operating procedures (SOPs) to minimize infection in the population. Based on these recommendations, some graphical results for the model with sensitive parameters are provided. 2020 Mathematics Subject Classifications: 26A33, 34A08, 03C65 Key Words and Phrases: SEI1I2HRD model, Two strain model, NSFD scheme, Sensitivity analysis, Effective reproduction number, Fractional numerical analysis ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6001 Email addresses: muhammadriaz84@gmail.com (M.Riaz), zeeshan@yuntech.edu.tw (Z. Ali), rahimmaths24@gmail.com (R. Ud Din), saahmad@psu.edu.sa (S. Ahmad), kamal@psu.edu.sa (K.Abodayeh), sarwarswati@gmail.com (M.Sarwar) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 2 of 27 1. Introduction In December 2019, it was reported that Wuhan, China, has been severely impacted by a new type of virus called Corona [1, 2]. The International Committee on Taxonomy of Viruses designated it as SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) in February 2020, as stated by [3–6]. First, the aforementioned virus struck Wuhan, China, and eventually spread to practically every city in the world as a result of the outbreak. Millions of people died as a result of it worldwide. It is difficult to pinpoint a single theory regarding the origin of the aforementioned virus. Some theories include the seafood market, human migration from one location to another brought on by animal transmission, or human-to-human contact. This virus has taken all the countries of the world in its grip. Virtually every aspect of life has been negatively impacted, including social life, health, the economy and education. Health researchers, national policymakers, and the medical community are at a loss as to how to contain this lethal outbreak. They are all watching the scenario from different points of view. They are making a lot of effort to at least reduce the number of fatalities brought on by this outbreak. Individuals who contract this pandemic often have mild respiratory symptoms. The signs of this condition include fever, fatigue, a dry cough, and a throat infection. A sore throat, pain, and a sinus infection are possible symptoms. In response to COVID-19, World Health Organization (WHO) played a vital role in the disease control. During pandemic WHO stay focused on surveillance, vaccination, and strengthening healthcare system. Their efforts provide a crucial framework for ongoing global cooperation to manage and control COVID-19 pandemic. For more detail regarding efforts of WHO for COVID-19 control, see [7, 8]. Mathematicians played a vital role during COVID-19 pandemic. They provide best predictive mathematical models to help medical researcher and policy makers in controlling the spread of the disease. For more detailed study, see [9–12]. Recently, non-integer order and nonlocal derivatives of fractional order have been added to the modelling domain [13– 16]. The dynamics of the spread of a particular disease are understood using mathematical models, which can be useful and provide a hand to save individuals from getting lost. By using the instruments indicated and this technique, you may create a plan to stop the disease from spreading or completely eradicate it. A number of models for COVID-19 have recently been constructed by numerous authors under various fractional calculus concepts. Very helpful models have been developed in this area; some of these are referred to as [17–20]. Since fractional-order systems are more thorough in examining real-world issues related to global dynamics, finding precise analytical solutions to fractal-order systems is becoming a challenging task. Some recent work cited as [21] where the authors developed a fractional order math- ematical model for Tuberculosis (TB) with incomplete treatment using the Atangana- Baleano-Caputo (ABC) derivate. Another useful work where the authors developed a fractional order mathematical model for the investigating the dynamics of tumor-immune cell interactions using fractional calculus, is cited as [22]. A model for dengue fever trans- mission incorporating the key factors like reinfection, vaccination, therapy, and asymptotic M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 3 of 27 carriers analyzed by the authors in [23] using fixed point theorem to prove the Ulam- Hyers stability. The authors analyzed a fractional order mathematical model for chronic myelogenous leukemia (CML) by proving existence, uniqueness, and Ulam-Hyers stability employing Banach and Schaefer fixed point theorems in [24]. For fractional order mathe- matical model of Human Immunodeficiency Virus (HIV) infection using Caputo-Fibrizio derivative the article [25] is cited hare for more information about application of fractional calculus to study mathematical models infectious disease. Incidence rates, a unique function included in epidemic models, describe the method of the disease’s transmission within a community or among a host. Numerous forms of mathematical models of infectious diseases, such as linear, bilinear, concave incidence rate, etc., were extensively employed. A convex incidence rate is one type of incidence rate that is defined as f(S, I) = KS(βS + 1). Where K and β are some constants that are positive, and S represents the susceptible compartment and I = I1 + I2 the infected compartment. The greatest rate of disease transmission brought on by two exposures over a brief period of time is known as the convex incidence rate. Whereas KSI2 results in a double infection, KSI leads to a single occurrence. The aforementioned idea will also be applied to our work. While cross-immunity and co-infection have been taken into account in studies on the dynamics of two viral infections [26], other investigations characterized the impact of two competing strains with cross-immunity [27]. According to reports, the number of strains actually grows logarithmically. The mean mortality rate and the maximum number of afflicted persons [28]. Among the main obstacles in the current COVID-19 pandemic are reinfection (perhaps brought on by declining immunity) and numerous virus strains [29]. The significance of integrating mutations and evolutionary adaptations into epidemic models was highlighted by Yagan et al. [30] due to the belief that certain COVID-19 strains are more transmissible than the original strain. Vaccination was not included in other research using the multi-strain epi- demic model for COVID-19 [31–33]. For numerical simulation, the NSFD scheme has been employed. This approach is a practical and effective way to solve many nonlinear problems numerically. For example, this technique has been applied by numerous researchers to the numerical simulation of numerous bathetical models [34, 35]. On the other hand [36] used fractional numerical simulation to discussed COVID-19 dynamics. Among the approaches discussed, the NSFD method is one of the most effective at estimating solutions to the fractional-order models; see [37, 38] for more information. The idea about two strains also originate from the article [39], where the authors de- velop a biologically realistic deterministic models of the transmission dynamics for multi- strain pathogens like influenza. Based on the genetic distance between strains, the authors modeled cross immunity and showed this affects the disease dynamics. A temporary broad (strain transcending) immunity which reduces infection prevalence and strain diversity were added to extend the model. The authors concluded that stochastic effects are impor- tant to fully capture the influenza evolution. Another article cited as [40] also plays a very M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 4 of 27 important role in our model. Where the authors formulated multi-strain model that links immunological and epidemiological dynamics across scales. Using a viral load dependent parameter, they showed that the strain with the highest reproduction number dominates the scales and Lyapunov function were used to prove the global stability. The article [41] provide useful information for optimizing testing strategies. About the key hosts in the evolution, resentment of swine influenza virus, and more about multi-strain pandemics, see [42–44]. About extension of model to account for multiple pathogen mutations and intervention policies, see [45]. In this paper, we constructed a mathematical SEI1I2HRD model to discussed the dynamics of COVID-19 under convex incidence rate. We formulated the model in seven equations, where S represent susceptible, E exposed, I1 strain one infected, I2 strain two infected, H hospitalized R recovered and D death compartment individual. For basic reproduction number R0 sensitivity analysis is done to find out the most sensitive variables which increase or decrease basic reproduction number R0. This paper aims to model the evolution of COVID-19 pandemic for two strains and complex transmission patterns using both integer-order and fractional-order derivatives using especially the convex incidence rate. The model capture memory and hereditary effects in the disease transmission which provide more realistic representation of how in- fection and immunity evolve over time compared to integer-order model. Using two strains allow the model to simulate variant competition which is very important for understanding the appearance of new variants COVID-19 virus. Using the convex incidence rate provide better saturation effects as compere to simple linear incidence rate. It is also evident that using convex incidence rate provide more realistic responses in case when infection level rises. Furthermore, integer and fractional-order model provide flexibility in theoretical frame work, that is, integer-order part where memory is negligible but fractional-order part is used in case where long term memory is important. No doubt the fractional- order mathematical model is more realistic but incorporating fractional calculus we may be facing with complex mathematics and computation. Moreover, simulating fractional- order models especially using convex incidence rate is computationally expensive and time consuming. This study provides a frame work for the policy makers to understand the dynamics of the disease transmission, propagation and recovery under various intervention scenarios. By incorporating and understanding the behavior of each parameter, the model can better simulate the potential impact of public health strategies. This results identifies the critical parameters that are responsible for the spread and control of the disease and also evaluates the effectiveness of the intervention strategies by looking at the values of the recovery rates γ1 and γ2. This study offers a decision support tool and helps the authorities to prepare for the reinfection due to loss of immunity. The results in this paper are organized as follows: Model formulation is discussed in Section 2. Mathematical results, local and global stability is discussed in Section 3. Numerical simulation through NSFD scheme and fractional order is in Section 4. The last Section is summarization of the results. M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 5 of 27 2. Model Formulation Owing to the novel COVID-19 variations strain1 and strain 2, individuals there were once again subject to a number of limitation, including the need of stick to SOPs religiously, keep social distance, wear facemasks, etc. A number of nations have prohibited travel to other countries due to the novel COVID-19 strains that is rapidly spreading. We develop a new mathematical model to comprehend the dyanimics of infection using strains I1 and I2. The whole population is taken to be M(t) = S + E + I1 + I2 + H + R + D where I = I1 + I2. We formulate the model using the following equations: dS dt = b−KSI(βI + 1)− d0S + rR, dE dt = KSI(βI + 1)− (d0 + η + ϕ)E, dI1 dt = ϕE − (γ1 + h1 + d0 + c1)I1, dI2 dt = ηE − (γ2 + h2 + d0 + c2)I2, dH dt = h1I1 + h2I2 − h0H, dR dt = γ1I1 + γ2I2 + h0H − (r + d0)R, dD dt = c1I1 + c2I2. (1) Where I = I1 + I2. A flowchart known as the compartment diagram illustrates how the population dynamics of various compartment have changed through time. Figure 1: Flow chart of parameters of the model (1). In the above flow chart, the compartment S represent a group of individuals that are not infected with a disease but are at risk to be infected. The individuals that are infected but not yet infectious moved from susceptible to exposed by a force of infection M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 6 of 27 KSI(βI+1) and are represented by E. I1 and I2 represents two strains of infection. Group of individuals that are died due to infection are represented by D. Moreover, H represents the individual that are admitted to a hospital who have developed severe symptoms from infection. The individuals that gaining immunity both temporary or permanent and so recovered from the disease is represented by R. The birth rate b sustains the susceptible population and so prolonging the epidemic. The natural death rate d0 reduce the population size and indirectly reduces the transmis- sion by shrinking the population at risk and the group of individuals that contributing in spreading the infection. The convex incidence rate KSI(βI + 1) shows that the infection force depends on the contact rate β and infectious individuals I. The progression rate ϕ from E to I1 accelerate the emergence of new infection while η the alternate infectious pathway basically making changes in the disease distribution. c1 nad c2 reduces the size of the infected compartments I1 and I2. h1 and h2 represents the movement of the infected individuals from I1 and I2 to hospitalize compartment. High hospitalize rates may indi- cates severe disease burden but a good isolation. The recovery perimeter h0 represent the rate of recover individuals from hospitalized compartment to recovered compartment R which reduce the pressure on healthcare system and policy makers. The most important parameters are the recovery parameters γ1 and γ2 which reduces the infectious population. Larger the value of these parameters shift individuals to the recover class which lead to controlling and limiting the outbreak of the disease. Moreover, the value of r is most important for the reinfection dynamics. Greater value of r leads to the multiple waves of infection. 3. Fundamental analysis of model (1) This section investigates the properties of the equilibrium point and the local and global stability of the system (1). We also investigate the general analysis of the system and the properties of the equilibrium. Lemma 1. The plane S+E+ I1+ I2+H +R+D = b d0 has invariant manifold of model (1), which is in first octant. Proof. By adding all vectors of model (1) and assumeM(t) = S+E+I1+I2+H+R+D, we have dM dt = b− d0S − d0E − d0I1 − d0I2 − d0H − d0R = b− d0(S + E + I1 + I2 +H +R). Implies that dM dt = b− d0M (2) M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 7 of 27 From (2), we get the solution is M(t) = b d0 . For any M(t0) ≥ 0, the general solution of (2) is M(t) = 1 d0 [b− (b− d0M(t0))e −d0(t−t0)]. Thus lim t→∞ M(t) = b d0 Which follow conclusion. Form above lemma, we can see that model (1) is well-posed mathematically and the following is feasible region for dynamical analysis. Φ = (S,E, I1, I2, H,R,D) : 0 ≤ S,E, I1, I2, H,R,D, S + E + I1 + I2 +H +R+D ≤ b d0 . 3.1. Uniqueness and Existence We look to find uniqueness and existence of system (1). For more detailed information, see [46, 47]. To accomplish this, we provide the model (6) in the manner specified by dS dt = f1 (t, S,E, I1, I2, H,R,D) , dE dt = f2 (t, S,E, I1, I2, H,R,D) , dI1 dt = f3 (t, S,E, I1, I2, H,R,D) , dI2 dt = f4 (t, S,E, I1, I2, H,R,D) , dH dt = f5 (t, S,E, I1, I2, H,R,D) , dR dt = f6 (t, S,E, I1, I2, H,R,D) , dD dt = f7 (t, S,E, I1, I2, H,R,D) . (3) Norm is define as ||χ||∞ = sup t∈x . (4) Where x ∈ [0, T ]. We supposed that S,E, I1, I2, H,R and D are bounded in [0, T ] and there exist k1, k2, k3, k4, k5, k6 and k7 for any t belong to [0, T ], such that: ||S||∞ < k1, ||E||∞ < k2, ||I1||∞ < k3, ||I2||∞ < k4, ||H||∞ < k5, ||R||∞ < k6, ||D||∞ < k7, (5) M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 8 of 27 and k = k1 + k2. Now, we have to prove (1) is bonded |f1 (t, S,E, I1, I2, H,R,D) | =|b−KSI(βI + 1)− d0S + rR|, ≤b+ |S|K|I|(β|I|+ 1) + d0|S|+ r|R|, ≤b+ sup M |S|K sup M |I|(β sup M |I|+ 1) + d0 sup M |S|+ r sup M |R|, ≤b+ ||S||∞K||I||∞(β||I||∞ + 1) + d0||S||∞ + r||R||∞, ≤b+ k1Kk(βk + 1) + d0k1 + rk6 < ∞. (6) Where M = t ∈ [0, T ]. Similarly with the same procedure we have |f2 (t, S,E, I1, I2, H,R,D) | = |KSI(βI + 1)− (d0 + η + ϕ)E|, ≤ Kk1k(βk + 1) + (d0 + η + ϕ)k2 < ∞. (7) |f3 (t, S,E, I1, I2, H,R,D) | = |ϕE − (γ1 + h1 + d0 + c1)I1|, ≤ ϕk2 + (γ1 + h1 + d0 + c1)k3 < ∞. (8) |f4 (t, S,E, I1, I2, H,R,D) | = |ηE − (γ2 + h2 + d0 + c2)I2|, ≤ ηk2 + (γ2 + h2 + d0 + c2)k4 < ∞. (9) |f5 (t, S,E, I1, I2, H,R,D) | = |h1I1 + h2I2 − h0H|, ≤ h1k3 + h2k4 + h0k5 < ∞. (10) |f6 (t, S,E, I1, I2, H,R,D) | = |γ1I1 + γ2I2 + h0H − (r + d0)R|, ≤ γ1k3 + γ2k4 + h0k5 + (r + d0)k6 < ∞. (11) |f7 (t, S,E, I1, I2, H,R,D) | = |c1I1 + c2I2|, ≤ c1k3 + c2k4 < ∞. (12) Thus S,E, I1, I2, H,R and D are bounded and there exist K1,K2,K3,K4,K5 and K6 such that  supM |f1 (t, S,E, I1, I2, H,R,D) |< K1, supM |f2 (t, S,E, I1, I2, H,R,D) |< K2, supM |f3 (t, S,E, I1, I2, H,R,D) |< K3, supM |f4 (t, S,E, I1, I2, H,R,D) |< K4, supM |f5 (t, S,E, I1, I2, H,R,D) |< K5, supM |f6 (t, S,E, I1, I2, H,R,D) |< K6. supM |f7 (t, S,E, I1, I2, H,R,D) |< K7. (13) Next, we have to prove |f1 (t, S1, E, I1, I2, H,R,D)− f1 (t, S2, E, I1, I2, H,R,D) | =  |b−KSI(βI + 1)− d0S + rR −(b−KSI(βI + 1)− d0S + rR)| < (KI(βI + 1) + d0)|S2 − S1| < p1|S2 − S1|. (14) M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 9 of 27 Where p1 = (KI(βI + 1) + d0). Similarly with the same method, we have |f2 (t, S,E1, I1, I2, H,R,D)− f2 (t, S,E2, I1, I2, H,R,D) | < p2|E2 − E1|, |f3 (t, S,E, I11 , I2, H,R,D)− f3 (t, S,E, I12 , I2, H,R,D) | < p3|I12 − I11 |, |f4 (t, S,E, I1, I21 , H,R,D)− f4 (t, S,E, I1, I22 , H,R,D) | < p4|I22 − I21 |, |f5 (t, S,E, I1, I2, H1, R,D)− f5 (t, S,E, I1, I2, H2, R,D) | < p5|H2 −H1|, |f6 (t, S,E, I1, I2, H,R1, D)− f6 (t, S,E, I1, I2, H,R2, D) |< p6|R2 −R1|. (15) Where p2 = d0 + ϕ+ η, p3 = γ1 + h1 + d0 + c1, p4 = γ2 + h2 + d0 + c2, p5 = h0, p6 = r + d0. Thus model (1) has a unique set of solutions as indicated by the uniqueness and existence. 3.2. Disease Free Equilibrium Points DFEP: DFEP of model (1) is represented by P0 = (S0, 0, 0, 0, 0, 0). Where P0 = ( b d0 , 0, 0, 0, 0, 0 ) . 3.3. Endemic Equilibrium points EEP We obtained EEP for our model (1) as S∗ = b(r+d0)+r[(γ1+h1)I∗1+(γ2+h2)I∗2 ] (r+d0)(KI∗(βI∗+1)+d0) , E∗= KI∗(βI∗+1)(b(r+d0)+r[(γ1+h1)I∗1+(γ2+h2)I∗2 ]) (r+d0)(d0+ϕ+η) , I∗1 = ϕKI∗(βI∗+1)(b(r+d0)+r[(γ1+h1)I∗1+(γ2+h2)I∗2 ]) (r+d0)(d0+ϕ+η)(γ1+h1+d0+c1) , I∗2 = ηKI∗(βI∗+1)(b(r+d0)+r[(γ1+h1)I∗1+(γ2+h2)I∗2 ]) (r+d0)(d0+ϕ+η)(γ2+h2+d0+c2) , H∗= h1I∗1+h2I∗2 h0 , R∗= (γ1+h1)I∗1+(γ2+h2)I∗2 r+d0 . (16) 3.4. Effective Reproduction number R0 The effective reproduction number, or R0 is a terminology in epidemiology that de- scribes how diseases are managed and spread. R0 provides information about the preva- lence of the disease in the population as well as the most effective measures for defending the local population against the fatal virus. By using the concept of next generation matrix for finding R0. M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 10 of 27 Let χ = (E, I1, I2), then form system (1), we have dχ dt = W −U, where W =  KSI(βI+ 1) 0 0  and U =  (d0 + η + ϕ)E −ϕE+ (γ1 + h1 + d0 + c1)I1 −ηE+ (γ2 + h2 + d0 + c2)I2  . Jacobian of W for the disease-free equilibrium is W =  0 kS0 kS0 0 0 0 0 0 0  and for the DFE Jacobian of U is given U =  (d0 + η + ϕ) 0 0 −ϕ γ1 + h1 + d0 + c1 0 −η 0 γ2 + h2 + d0 + c2  . WU−1 =  0 KS0 KS0 0 0 0 0 0 0  q(γ1 + h1 + d0 + c1) 0 0 ϕ(γ2 + h2 + d0 + c2) q(γ2 + h2 + d0 + c2) 0 η(γ1 + h1 + d0 + c1) 0 q(γ1 + h1 + d0 + c1)  (q(γ1 + h1 + d0 + c1)(γ2 + h2 + d0 + c2) . Where q = d0 + η + ϕ, and WU−1 is the next generation matrix. The non-negative eigenvalues are R01 = Kb d0(d0 + ϕ+ η)(γ2 + h2 + d0 + c2) R02 = Kbϕ d0(d0 + ϕ+ η)(γ1 + h1 + d0 + c1) . R01 represent the effective reproduction of first strain of COVID-19 and R02 represent the effective reproduction of second strain of COVID-19. Hence R0 = max [R01,R02] . (17) M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 11 of 27 4. Stability Analysis This section investigates the properties of the equilibrium point and the local and global stability of the system (1). We provide references to relevant papers as [48, 49] where a known or standard result has been used as a basis for our proofs. 4.1. Local Stability In this section of our study, we go over the system’s local stability (1). The follow- ing theorem is present to obtain the necessary conclusions while employing the Jacobian matrix. Theorem 1. [rgb]1.00,1.00,1.00“ When effective basic reproduction number R01 < 1 or R02 < 1, the DFEP “P” is stable locally asymptotically. ” Proof. Jacobian Matrix as J 0 =  −d0 0 0 0 0 0 r 0 −(d0 + ϕ+ η) 0 0 0 0 0 0 ϕ −(γ1 + h1 + d0 + c1) 0 0 0 0 0 η 0 −(γ2 + h2 + d0 + c2) 0 0 0 0 0 h1 h2 −h0 0 0 0 0 γ1 γ2 h0 −(r − d0) 0 0 0 c1 c2 0 0 0  . and |λ− J 0(Q)| = 0. Hence λ1 = −d0, λ2 = −(d0 + ϕ+ η), λ3 = −(γ2 + h2 + d0 + c2), λ4 = −(γ1 + h1 + d0 + c1), λ5 = −(r + d0), λ6 = (γ2 + h2 + d0 + c2)(R01 − 1), λ7 = (γ1 + h1 + d0 + c1)(R02 − 1). Hence λ1, λ2, λ3, λ4, λ5 are strictly negative, but λ6 and λ7 are negative whenever R01 < 1 and R02 < 1. Which conclude conclusion. 4.2. Global Stability In this section of our paper, we calculated global stability of system (1) at DFEP “P”. To study global stability, we constricted a function called Lyapunov. We need to show Ḟ = −kF , where k is a constant. M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 12 of 27 Theorem 2. [rgb]1.00,1.00,1.00“When effective basic reproduction number R01 < 1 or R02 < 1, the DFEP “P” is stable globally asymptotically.” Proof. Consider the function F = k1I1 + k1I2, (18) where k1 and k2 are some constants. Since I1, I2 > 0, then F > 0 and F attains zero at I1 = I2 = 0. Now, our goal is Ḟ < 0. dF dt = k1 dI1 dt + k2 dI2 dt . (19) Putting values of I1 and I2 from (1) in (19), we get dF dt = k1 (ϕE− (γ1 + h1 + d0 + c1)I1) + k2 (ηE− (γ2 + h2 + d0 + c2)I2) . (20) By putting k1 = k2 = kb d0(d0+ϕ+η) , and rearrange (20), we get dF dt = (γ1 + h1 + d0 + c1)I1 (R01 − 1) + (γ2 + h2 + d0 + c2)I2 (R01 − 1) . (21) By using “LaSalle’s invariant principle”, dF dt = 0 iff I1 = I2 = 0. Furthermore dF dt < 0 whenever R01 < 1 and R02 < 1, which follow conclusion. 5. Sensitivity Analysis We present the sensitivity analysis in this section to take into consideration the vari- ables that have a major impact on the fundamental reproduction number. For more de- tailed information about sensitivity analysis, see [50]. Sensitivity analysis is recommended in order to determine the relative significance of the various factors impacting the preva- lence and transmission of sickness. It is imperative to regulate the variations in the model (1) parameters in order to get R01 and R02 < 1, as well as to stop the spread of diseases. The ratio of a variable’s change to a parameter’s change is known as the sensitivity index, and it may be calculated using the formula. The model’s parameters’ sensitivity indices for R01 < 1 and R02 < 1 are given in Table 1 using normalized sensitivity index defined by f [y] = y R0 ∂R0 ∂y . M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 13 of 27 Parameters Sensitivity Index Parameters Sensitivity Index f(b) 1 f(γ1) 0.8737 f(K) 1 f(γ2) 0.6978 f(h1) −0.00124 f(c1) −0.00103 f(h2) −0.00124 f(c2) −0.0015 f(d0) 0.9674 f(η) −0.1044 f(ϕ) −0.9895 Table 1: Sensitivity indices of R0 with respect to model parameters . -1.5 -1 -0.5 0 0.5 1 1.5 P a ra m te rs v a lu e s Axis Title Figure 2: Sensitivity indices of R0 with respect to model (1) parameters. (1) and (2) make it evident that the disease transmission rate of strain 2 η, contact rate K, death rate d0, infection rate of strain 1 ϕ and rate of recovery individuals γ1 and γ2 are the parameters of the model (1) that are most sensitive to R01 and R02. As the value of infection rate of strain 1 ϕ decrease, so does the value of R01. When infectious rate K, death rate d0 decrease. As a result, R01 and R02 rises infection rate of strain 1 ϕ and infection rate strains 2 η decline. Clearly, fig. 1 shows that contact rate K, death rate d0, birth b rate, and rate of recovery of individuals γ1 and γ2 are positive, so an increase in their values increased R01 and R02. On the other hand, the values of infection rate of strain 1 ϕ, infection rate of strain 2 η, strain 1 hospitalization rate h1, strain 2 hospitalization rate h2, death rate of strain 1 c1, and death rate of strain 1 c2 are negative, so their increased show decreased in R01 and R02. It is clear that a lower value of R01 and R02 is a contributing factor to the decline in infection incidence. We need to preserve parameters with negative indices while lowering the effectiveness of parameters with positive indices to a basic reproduction number in order to completely remove the infection from the system. Therefore, health officials need to carefully consider any preventative steps that would lower the infection burden in order to reduce future outbreaks. We found that the control parameters, such as rate of recovery individuals γ1 and γ2, etc., which are adversely linked with R01 and R02, should be implemented with adequate hygiene and efficient health care services. M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 14 of 27 6. Numerical Results and Discussion Here, we present simulations of model (1) to confirm the analytical analysis. We used NSFD [51, 52] for simulation. The first equation of model (1) is rewrite in difference as dS dt = b−KSI(βI + 1)− d0S + rR. (22) By using the none-standard finite difference method, we decomposed as Sj+1 − Sj h = b−KSjIj(βIj + 1)− d0Sj + rRj . (23) Like (23), we decomposed other equations of model (1) by using none-standard finite difference method Ej+1 = Ej + h ( KSjIj(βIj + 1)− (d0 + η + ϕ)Ej ) , I1(j+1) = I1(j) + h ( ϕEj − (γ1 + h1 + d0 + c1)I1j ) , I2(j+1) = I2(j) + h ( ηEj − (γ2 + h2 + d0 + c2)Ij ) , Hj+1 = Hj + h ( h1I1j + h2I2j − h0Hj ) , Rj+1 = Rj + h ( γ1I1j + γ2I2j + h0Hj − (r + d0)Rj ) , Dj+1 = Dj + h ( c1I1j + c2I2j ) . By using NSFD scheme, we simulate the model (1) by using numerical values given in Table 2 taking the values from [4]. Variable numerical values Variable numerical S 60140000 ϕ 0.9566 E 62000 γ1 0.8447 I1 8000 γ2 0.6746 I2 100 c1 0.0012 H 50 c2 0.0015 R 0 η 0.01010 D 0 b 2560 d0 39% β 0.00231% K 12% r 0.046 h1 0.0010 h2 0.0010 h0 0.0200 d0 0.00004 Table 2: The parameters of the system are described and specified, with approximate real values (1). M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 15 of 27 Here in figures 3-9, we simulate the results using the numerical values of Table 2. In steady state, the number of stage one carcinoma patients with an initial condition drops from 20 to 3 in three to four months of time. At stage two, it decreased to 4 from 14 in three to four months of time. Stage 3 decreased to 5 from 30 in three to four months of time. In Stage 4, it decreased to 5 from 20 in three to four months of time. In the recovered population, a rapid increase was seen, and cardiotoxic patients showed a great increase from 10 to 29 in the first month under steady-state conditions. Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f S us ce pt ib le 0 10 20 30 40 50 60 70 Figure 3: Plot of numerical solutions for the susceptible class of the proposed model (1). Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f E xp os ed 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Figure 4: Plot of numerical solutions for the exposed class of the proposed model (1). M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 16 of 27 Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f s tr ai n 1 co m pa rt m en t 0 0.005 0.01 0.015 0.02 0.025 0.03 Figure 5: Plot of numerical solutions for the strain 1 class of the proposed model (1). Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f s tr ai n 2 co m pa rt m en t ×10-4 0 0.5 1 1.5 2 2.5 3 3.5 Figure 6: Plot of numerical solutions for the strain 2 class of the proposed model (1). Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty H os pi ta liz ed c om pa rt m en t ×10-4 0 1 2 3 4 5 6 Figure 7: Plot of numerical solutions for the hospitalized class of the proposed model (1). M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 17 of 27 Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f R ec ov er ed c om pa rt m en t 0 0.01 0.02 0.03 0.04 0.05 0.06 Figure 8: Plot of numerical solutions for the recovered class of the proposed model (1). Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty d ea th c om pa rt m en t ×10-4 0 1 2 3 4 5 6 Figure 9: Plot of numerical solutions for the death class of the proposed model (1). Here, the dynamical behaviors corresponding to different initial data of each compart- ment of the proposed model is presented in figures 3 to 9. With the help of treatment, social distancing, and other intervention strategies, the susceptible population decreases over time and finally reaching to an equilibrium position which also reduce the rate of infection in other compartments, eventually the disease is going to be decline. Moreover, the decline in susceptible class is shown in figure 3. Further, the decline in exposed com- partment is shown in 4. Dynamical behaviors of strain one and two are shown in figures 5 and 6 respectively. In addition, the dynamical growth in, hospitalization, recovered and death class has shown in figure 7, 8 and 9 respectively. Table 2 provides the parameter values that were taken into consideration for this nu- merical simulation. We select the time unit per day for this numerical simulation. figures 3 to 9 display the simulation results for model (1). The dynamics of the infected compart- ments with varying values of the parameter K are depicted in figure 3 and its subgraphs. As contact decreases, so does the number of infected people in each compartment. The M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 18 of 27 numerical simulation indicates that minimizing contact between healthy individuals in- fected with strains 1 and 2 is necessary to reduce the number of infected cases within the community. Both strain 1 and strain 2 infected people can spread the infection to other exposed people who show any symptoms of the disease. Owing to the strain 1 virus’s easy and high rate of transmission from an infected person, it is thought to be more harm- ful. The WHO’s recommendations can help prevent infection by reducing the number of instances of disease. The infected compartments and different values of the parameter ϕ are displayed in figure 5. As the ϕ value declines, the number of afflicted individuals also gradually declines. To stop other healthy people from limiting the disease’s spread, those who exhibit clinical signs should be admitted to the hospital. The spread of strain 1 infections among those who do not exhibit any clinical symptoms is the most critical, and limiting the contacts of those who are infected with other members of the community should be the first goal in order to prevent this spread. Figures 6 and 7 illustrate how the parameters η and h0 affect the model’s infected compartments. From the data shown in Figs. 8 and 9, it is evident that as the parameters’ values decreased, so did the individuals’ infected compartments. 7. Numerical Analysis through Fractional Calculus In this section of our work, we use fractional calculus to explore the COVID-19 model (1). Here, we recall fractional derivative and integral [52] as: Definition 1. A function S : [0,∞) → R is define with fractional order integral χ > 0 as Iφt S(t) = 1 β(γ) ∫ t 0 S(α) (t− α)1−γ dα, provided the integral exists at the right sides. Moreover, by definition of Caputo, we have Dγ 0+S(t) =  1 β(1− γ) ∫ t 0 (t− α)−γS ′ (α)dη, 0 < γ ≤ 1, dS dt , γ = 1. For fractional order Dγ tX (t) = H(t,X (t)), X (0) = X0. (24) Where X = (S,E, I1, I2,H,R,D) and X0 = (S0,E0, I10 , I20 ,H0,R0,D0) . M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 19 of 27 Then using the generalized form of Taylor series is follow X (t+ h) = X (t) + hγ α(γ + 1) Dγ tX (t) + h2γ χ(2γ + 1) D2γ t t(t) + ... (25) By using fractional calculus results, we have D2γ t X = Dγ tH(t,X ) + H(t,X )Dγ XH(t,X ). Used (25), and some rearrangements, we have, X (t+ h) = X (t) + hγ α(γ + 1) H(t,X (t)) + hγ 2α(γ + 1) [ P1 +P2 ] , (26) where P1 = H(tj ,Xj(tj)), P2 = H ( tj + 2hγα(γ + 1) α(2γ + 1) , X (tj) + 2hγα(γ + 1) α(2γ + 1) H (tj ,Xj(tj)) ) . Here we present the numerical illustration against fractional orders taken as 0.6, 0.62, ..., 1.0 in figures 10-16 respectively Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f S us ce pt ib le 0 10 20 30 40 50 60 70 Figure 10: Plot of numerical solutions for the susceptible class of the proposed model (1) for various fractional orders as 0.6, 0.62, ..., 1.0 . M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 20 of 27 Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f E xp os ed 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Figure 11: Plot of numerical solutions for the exposed class of the proposed model (1)for various fractional orders as 0.6, 0.62, ..., 1.0 . Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f s tr ai n 1 co m pa rt m en t 0 0.005 0.01 0.015 0.02 0.025 0.03 Figure 12: Plot of numerical solutions for the strain 1 class of the proposed model (1)for various fractional orders as 0.6, 0.62, ..., 1.0 . M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 21 of 27 Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f s tr ai n 2 co m pa rt m en t ×10-4 0 0.5 1 1.5 2 2.5 3 3.5 Figure 13: Plot of numerical solutions for the strain 2 class of the proposed model (1)for various fractional orders as 0.6, 0.62, ..., 1.0 . Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty H os pi ta liz ed c om pa rt m en t ×10-4 0 1 2 3 4 5 6 Figure 14: Plot of numerical solutions for the hospitalize class of the proposed model (1)for various fractional orders as 0.6, 0.62, ..., 1.0 . M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 22 of 27 Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty o f R ec ov er ed c om pa rt m en t 0 0.01 0.02 0.03 0.04 0.05 0.06 Figure 15: Plot of numerical solutions for the recovered class of the proposed model (1)for various fractional orders as 0.6, 0.62, ..., 1.0 . Time (Days) 0 50 100 150 200 250 P op ul at io n D en si ty d ea th c om pa rt m en t ×10-4 0 1 2 3 4 5 6 Figure 16: Plot of numerical solutions for the death class of the proposed model (1) for various fractional orders as 0.6, 0.62, ..., 1.0 . Using the numerical values of the model parameters provided in 2, we consider the fractional order model considered under the Caputo derivative to generate the graphical results by the above-mentioned scheme. The results of applying the numerical technique for the Caputo fractional order, as stated in equation (25) are displayed in figures 10 - 16. By taking into account the fractional-order values 0.6, 0.62, ..., 1.0, we plot graphical results, where a rise was noted for the model’s infected compartments and a drop in the fractional-order parameters was observed for the number of susceptible and recovered people. SEIR epidemic model consists only on Susceptible, Exposed, Infected, and Recovered classes, while in our model, we split I into two strains I1 and I2, and added D (death due to corona) and H (hospitalized individuals). Incorporating these new compartments allows the model to simulate disease progression more accurately. The hospitalized com- partment is crucial for the analysis of hospitalized rate of infected individual and load on M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 23 of 27 the healthcare departments. The compartment D accurately separates the death due to corona from the natural death and allows evaluation of how treatment and other interven- tion strategies reduce the number of deaths due to infection. Comparing to SEIR epidemic model, our model offers a more realistic and flexible framework to better understand the transmission and burden of COVID-19. 8. Conclusion In this work, a novel mathematical model has been introduced to comprehend the two strain variations of COVID-19 and show how exposed, asymptomatic, symptomatic, and hospitalized persons interact with COVID-19. We developed the model and provided a thorough explanation of the biological procedure used to model the two strains. In the disease-free scenario, we found that the two-strain model is globally asymptotically stable if R01 ≤ 1 or R02 ≤ 1 and locally asymptotically stable if R01 < 1 or R02. The two strain model’s global stability suggests that the only way to contain the infection is to lower its basic reproduction number R01 or R02. Basic reproduction can be easily managed with less than unity through social distancing, face mask use, frequent hand washing, staying at home, and other practices. To identify the parameters that are most susceptible to differences in R0, a global sensitivity analysis is conducted. According to the sensitivity analysis, interaction between healthy individuals and strains 1 and 2 can spread the virus throughout the nation. Therefore, in order to stop the strains from spreading to other people, it is imperative that the WHO’s instructions are properly followed. Using the values of the real parameters that are received from the data fitting, a unique numerical approach is applied to generate the graphical result for the model. There were several graphical representations of the model parameters and how they affected the elimination of the diseases. It is evident from the stated outcomes that the disease can be eradicated according to the World Health Organization’s (WHO) guidelines. The model that investigated here in this article is useful for studying the dynamics of two circulating COVID-19 strains both for integer and fractional order. The reason of incorporating the convex incidence rate is to capture nonlinear transmission effect such as behavioral changes and saturation, offering more realistic disease progression than tra- ditional bilinear model. Moreover, the use of fractional calculus adds memory effects, making the model specifically relevant for long-term behavior, delayed responses or linger- ing immunity. Analytical and numerical results provide insight into strain competition, persistence, and potential resurgence patterns. Moreover, there are some limitations that affects its real world applicability. Critical factors like vaccination, stochastic effects, ongoing strain evolution, testing, and demo- graphic structure. These limitations provide direction for future research where these factors can be incorporated to enhance the model accuracy. Acknowledgements Authors are thankful to Prince Sultan University for paying the APC and support. M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 24 of 27 References [1] Coronavirus, Ethiopia. 13,968 cases and 223 deaths. Worldometers, 2020. Accessed on 27 August 2020, https://www.worldometers.info/coronavirus/ country/ethiopia. [2] J. T. Wu, K. Leung, and G. 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. [3] S. Zhao, Q. Lin, J. Ran, S. S. Musa, G. Yang, W. Wang, Y. Lou, D. Gao, L. Yang, D. He, and M. H. Wang. Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, from 2019 to 2020: a data-driven analysis in the early phase of the outbreak. International Journal of Infectious Diseases, 92:214–217, 2020. [4] Q. Li, X. Guan, P. Wu, X. Wang, L. Zhou, Y. Tong, R. Ren, K. S. Leung, E. H. Lau, J. Y. Wong, and X. Xing. Early transmission dynamics in Wuhan, China, of novel coronavirus–infected pneumonia. New England Journal of Medicine, 382(13):1199– 1207, 2020. [5] A. B. Gumel, S. Ruan, T. Day, J. Watmough, F. Brauer, P. Van den Driessche, D. Gabrielson, C. Bowman, M. E. Alexander, S. Ardal, and B. M. Sahai. Modelling strategies for controlling SARS outbreaks. Proceedings of the Royal Society of London. Series B: Biological Sciences, 271(1554):2223–2232, 2004. [6] G. Rahman, K. Shah, F. Haq, and N. Ahmad. Host vector dynamics of pine wilt disease model with convex incidence rate. Chaos, Solitons & Fractals, 113:31–39, 2018. [7] World Health Organization. 2019 novel coronavirus (2019-nCoV): strategic preparedness and response plan. World Health Or- ganization, 2020. https://www.who.int/publications/i/item/ strategic-preparedness-and-response-plan-for-the-new-coronavirus. [8] P. M. Nischal. WHO announces end of the COVID-19 global health emergency. The National Medical Journal of India, 36(3):137–138, 2023. [9] R. Ud Din and E. A. Algehyne. Mathematical analysis of COVID-19 by using SIR model with convex incidence rate. Results in Physics, 23:103970, 2021. [10] Q. Lin, S. Zhao, D. Gao, Y. Lou, S. Yang, S. S. Musa, M. H. Wang, Y. Cai, W. Wang, L. Yang, and D. He. A conceptual model for the coronavirus disease 2019 (COVID- 19) outbreak in Wuhan, China with individual reaction and governmental action. International Journal of Infectious Diseases, 93:211–216, 2020. [11] D. S. Hui, E. I. Azhar, T. A. Madani, F. Ntoumi, R. Kock, O. Dar, G. Ippolito, T. D. Mchugh, Z. A. Memish, C. Drosten, and E. Petersen. The continuing 2019- nCoV epidemic threat of novel coronaviruses to global health—the latest 2019 novel coronavirus outbreak in Wuhan, China. International Journal of Infectious Diseases, 91:264–266, 2020. [12] R. Ud Din, K. Shah, I. Ahmad, and T. Abdeljawad. Study of transmission dynamics of novel COVID-19 by using mathematical model. Advances in Difference Equations, M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 25 of 27 2020:323, 2020. [13] K. Shah, M. A. Alqudah, F. Jarad, and T. Abdeljawad. Semi-analytical study of Pine Wilt Disease model with convex rate under Caputo–Fabrizio fractional order derivative. Chaos, Solitons & Fractals, 135:109754, 2020. [14] K. Shah, A. Khan, B. Abdalla, T. Abdeljawad, and K. A. Khan. A mathematical model for Nipah virus disease by using piecewise fractional order Caputo derivative. Fractals, 32(2):2440013, 2024. [15] S. T. Thabet, M. S. Abdo, K. Shah, and T. Abdeljawad. Study of transmission dynamics of COVID-19 mathematical model under ABC fractional order derivative. Results in Physics, 19:103507, 2020. [16] E. Atangana and A. Atangana. Facemasks simple but powerful weapons to protect against COVID-19 spread: can they have side effects? Results in Physics, 19:103425, 2020. [17] A. Atangana. Fractional discretization: the African’s tortoise walk. Chaos, Solitons & Fractals, 130:109399, 2020. [18] D. Kumar, J. Singh, M. Al-Qurashi, and D. Baleanu. A new fractional SIRS-SI malaria disease model with application of vaccines, antimalarial drugs, and spraying. Advances in Difference Equations, 2019(1):278, 2019. [19] A. J. Arenas, G. González-Parra, and B. M. Chen-Charpentier. Construction of nonstandard finite difference schemes for the SI and SIR epidemic models of fractional order. Mathematics and Computers in Simulation, 121:48–63, 2016. [20] M. T. Hoang and O. F. Egbelowo. Dynamics of a fractional-order hepatitis B epidemic model and its solutions by nonstandard numerical schemes. pages 127–153, 2020. [21] M. U. Rahman, M. Arfan, Z. Shah, P. Kumam, and M. Shutaywi. Nonlinear frac- tional mathematical model of tuberculosis (TB) disease with incomplete treatment under Atangana-Baleanu derivative. Alexandria Engineering Journal, 60(3):2845– 2856, 2021. [22] T. Q. Tang, Z. Shah, R. Jan, and E. Alzahrani. Modeling the dynamics of tumor– immune cells interactions via fractional calculus. The European Physical Journal Plus, 137(3):367, 2022. [23] T. Q. Tang, R. Jan, E. Bonyah, Z. Shah, and E. Alzahrani. Qualitative analysis of the transmission dynamics of dengue with the effect of memory, reinfection, and vaccination. Computational and Mathematical Methods in Medicine, 2022:7893570, 2022. [24] T. Q. Tang, R. Jan, Z. Ur Rehman, Z. Shah, N. Vrinceanu, and M. Racheriu. Modeling the dynamics of chronic myelogenous leukemia through fractional-calculus. Fractals, 30(10):2240262, 2022. [25] Z. Shah, R. Jan, P. Kumam, W. Deebani, and M. Shutaywi. Fractional dynamics of HIV with source term for the supply of new CD4+ T-cells depending on the viral load via Caputo–Fabrizio derivative. Molecules, 26(6):1806, 2021. [26] L. J. Allen, M. Langlais, and C. J. Phillips. The dynamics of two viral infections in a single host population with applications to hantavirus. Mathematical Biosciences, 186(2):191–217, 2003. M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 26 of 27 [27] M. Nuno, G. Chowell, X. Wang, and C. Castillo-Chavez. On the role of cross-immunity and vaccines on the survival of less fit flu-strains. Theoretical Population Biology, 71(1):20–29, 2007. [28] T. Lazebnik and S. Bunimovich-Mendrazitsky. Generic approach for mathematical model of multi-strain pandemics. PLoS ONE, 17(4):e0260683, 2022. [29] E. F. Arruda, S. S. Das, C. M. Dias, and D. H. Pastore. Modelling and optimal control of multi strain epidemics, with application to COVID-19. PLoS ONE, 16(9):e0257512, 2021. [30] O. Yagan, A. Sridhar, R. Eletreby, S. Levin, J. B. Plotkin, and H. V. Poor. Modeling and analysis of the spread of COVID-19 under a multiple-strain model with mutations. Harvard Data Science Review, 4(1), 2021. [31] M. Massard, R. Eftimie, A. Perasso, and B. Saussereau. A multi-strain epidemic model for COVID-19 with infected and asymptomatic cases: application to French data. Journal of Theoretical Biology, 545:111117, 2022. [32] O. Khyar and K. Allali. Global dynamics of a multi-strain SEIR epidemic model with general incidence rates: application to COVID-19 pandemic. Nonlinear Dynamics, 102(1):489–509, 2020. [33] Y. Lou and R. B. Salako. Control strategies for a multi-strain epidemic model. Bulletin of Mathematical Biology, 84(10):110, 2022. [34] E. A. Algehyne and R. Ud Din. On global dynamics of COVID-19 by using SQIR type model under non-linear saturated incidence rate. Alexandria Engineering Journal, 60(1):393–399, 2021. [35] K. Shah, T. Abdeljawad, and R. Ud Din. To study the transmission dynamic of SARS- CoV-2 using nonlinear saturated incidence rate. Physica A: Statistical Mechanics and its Applications, 604:127915, 2022. [36] P. M. Manning and G. F. Margrave. Introduction to non-standard finite-difference modelling. CREWES Research Report, 18:1–10, 2006. [37] M. Arfan, M. M. Lashin, P. Sunthrayuth, K. Shah, A. Ullah, K. Iskakova, M. R. Gorji, and T. Abdeljawad. On nonlinear dynamics of COVID-19 disease model corre- sponding to nonsingular fractional order derivative. Medical & Biological Engineering & Computing, 60(11):3169–3185, 2022. [38] M. Arfan, K. Shah, and A. Ullah. Some theoretical and computational results about COVID-19 by using a fractional-order mathematical model. In Fractional-Order Mod- eling of Dynamic Systems with Applications in Optimization, Signal Processing and Control, pages 37–68. Academic Press, 2022. [39] P. Minayev and N. Ferguson. Improving the realism of deterministic multi-strain models: implications for modelling influenza A. Journal of the Royal Society Interface, 6(35):509–518, 2009. [40] Y. X. Dang, X. Z. Li, and M. Martcheva. Competitive exclusion in a multi-strain immuno-epidemiological influenza model with environmental transmission. Journal of Biological Dynamics, 10(1):416–456, 2016. [41] L. Shami and T. Lazebnik. Economic aspects of the detection of new strains in a multi- strain epidemiological–mathematical model. Chaos, Solitons & Fractals, 165:112823, M. Riaz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6001 27 of 27 2022. [42] I. H. Brown. The pig as an intermediate host for influenza A viruses between birds and humans. International Congress Series, 1219:173–178, 2001. [43] A. Alexi, A. Rosenfeld, and T. Lazebnik. Multi-species prey–predator dynamics during a multi-strain pandemic. Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(7):073111, 2023. [44] M. Martcheva. A non-autonomous multi-strain SIS epidemic model. Journal of Biological Dynamics, 3(2-3):235–251, 2009. [45] T. Lazebnik and G. Blumrosen. Advanced multi-mutation with intervention policies pandemic model. IEEE Access, 10:22769–22781, 2022. [46] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo. Theory and applications of fractional differential equations. 204, 2006. [47] R. P. Agarwal, M. Meehan, and D. O’Regan. Fixed point theory and applications. 141, 2001. [48] H. K. Khalil and J. W. Grizzle. Nonlinear systems. 3, 2002. [49] A. Korobeinikov and P. K. Maini. A Lyapunov function and global properties for SIR and SEIR epidemiological models with nonlinear incidence. Mathematical Biosciences and Engineering, 1(1):57–60, 2004. [50] N. Chitnis, J. M. Hyman, and J. M. Cushing. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bulletin of Mathematical Biology, 70(5):1272–1296, 2008. [51] Z. J. Fu, Z. C. Tang, H. T. Zhao, P. W. Li, and T. Rabczuk. Numerical solutions of the coupled unsteady nonlinear convection-diffusion equations based on generalized finite difference method. The European Physical Journal Plus, 134(6):272, 2019. [52] R. Ud Din, K. A. Khan, A. Aloqaily, N. Mlaiki, and H. Alrabaiah. Using non-standard finite difference scheme to study classical and fractional order SEIVR model. Fractal and Fractional, 7(7):552, 2023.