EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6005 ISSN 1307-5543 – ejpam.com Published by New York Business Global Numerical Analysis of Deterministic and Stochastic Model of COVID-19 Co-infection with Influenza Rahim Ud Din1,∗, Muhammad Shoaib Arif2 1 Department of Mathematics, Air University, PAF Complex, Islamabad 44000, Pakistan, 2 Department of Mathematics and Sciences, College of Humanities and Sciences, Prince Sultan University, Riyadh, 11586, Saudi Arabia. Abstract. The dynamics of co-infection model of SARS-Cov-2 and influenza is presented in this paper. A detailed analysis is conducted on the possible effects of the influenza vaccination alone as well as the combined effect of both vaccinations on the co-infection dynamics. The two diseases’ basic reproduction numbers utilizing the next-generation matrix method. Endemic equilibrium point (EEP), and disease free equilibrium points (DFEP) are calculated for deterministic model. The global stability of the model equilibrium is demonstrated using the Lyapunov function function, and the local stability is discussed for deterministic as well as stochastic. There are supplied simulations to highlights the theoretical outcomes of the model. For each infection, we compute the basic reproductive numbers R01 < 1, and R02 < 1. To determine the realistic values of the model parameters, data regarding COVID-19 instances was obtained from China National Health Commission. Sensitivity study employing the PRCC technique to look into the key variables that affect R01, or R02 < 1 decrease or increase. We show the numerical simulation of the model using non-standard finite difference scheme (NSFD) and stochastic numerical method. On the basis of these mentioned numerical methods, some graphical results for the model with sensitive parameters are provided. 2020 Mathematics Subject Classifications: 92D30, 60H10, 65C30, 34F05 Key Words and Phrases: Co-infection model, NSFD scheme, Sensitivity analysis, Stochastic model, Stochastic numerical analysis 1. Introduction In December 2019, a city in China known as Wuhan identified the coronavirus for the first time [1, 2]. As of November 2022, the new infection has caused over six hun- dred million cases and over seven million fatalities worldwide in a short period of time [3]. Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) is the cause of COVID-19 disease, the common cold, and Middle East Respiratory Syndrome (MERS) ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6005 Email addresses: rahimmaths24@gmail.com (Rahim ud din), marif@psu.edu.sa (Muhammad Shoaib Arif) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 2 of 27 [4]. However, COVID-19 is more contagious and deadly than the aforementioned illnesses. The World Health Organization (WHO) considered the COVID-19 outbreak to be a pan- demic in March 2020 [4, 5]. Direct transmission of the virus can occur from one person to another via breathing in harmful particles generated when an infected individual meets each other, talks, coughs, or sneezes. If a susceptible individual comes into contact with somebody who is infected through touching, holding, kissing, or other behaviors, they could become infected with the virus. The indirect transmission mechanism involves con- tacting a contaminated surface with the hands to touch the eyes, nose, and mouth [6, 7]. Prior to the introduction of vaccines, government officials and decision-makers advised the public to focus on non-pharmaceutical means. Good personal hygiene, using the right sneezing mode, physical distancing, closing schools, isolations, hospitalization of exposed people and infected individuals, contact tracing, and donning face masks are some of these non-pharmaceutical interventions [8, 9]. To reduce the number of cases, other measures included screening exposed persons, providing rapid immunization, and providing intense medical care [10, 11]. In many countries throughout the world, infection controls and preventive measures were in place, but the disease continued to spread. It started with alpha, went on to beta and delta, and is still running strong with the latest omicron version [12, 13], which significantly raises the number of cases and death rate. Public health professionals have long been concerned about the bacterial influenza disease often known as the flu. Seasonal influenza outbreaks affect millions of people each year, with an estimated 500,000 deaths from the illness [14, 15]. It has been demonstrated that only close interaction with droplets and interaction with an infected individual can carry in- fluenza dynamics, despite mounting evidence suggesting aerosols play a significant role in the spread of the illness [16, 17]. Muscle aches, loss of hunger, temperature, sore head, and sore throats are typical symptoms. The cyclical nature of outbreaks of influenza and the persistent frequency of COVID-19 may cause influenza and SARS-Cov2 to revolve, greatly increasing the likelihood of co- infection. The research that is currently applicable indicates that co-infection with the SARS-Cov2 increases influenza infectivity in an extensive range of cell types [18] and that people with co-infection seem to have same clinical signs and imaging results compared with individuals infected with SARS-Cov-2 only [19]. However, there is limited informa- tion available regarding the biological and medical consequences of co-infection. Seasonal influenza vaccination is considered to play a major impact in influenza prevention as well as offer extra advantages during the difficult COVID-19 pandemic period. The influenza vaccine considerably reduces the risk of respiratory infectious disorders such as COVID-19, which also relieves pressure on medical facilities and save materials for treating more seri- ous infection [20]. Furthermore, it is believed that vaccination against influenza enhances COVID-19 surveillance’s specificity and accuracy [21]. The co-infection of COVID-19 and influenza virus was investigated by [22] using two mathematical models of COVID-19 and the co-infection idea (COVID-19 and influenza virus). A unique characteristic of epidemic models, incidence rates explain how a disease spreads within a population or among a host. There are many types of interaction rates of infectious disease in mathematical models, including concave, bilinear, square root, monotonic, saturated and linear models, R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 3 of 27 the following is used in our work, which is called convex incidence rate. Follow as, f(S, I) = KS(βS + 1). Here, β > 0 and K > 0 are belong to positive real constant. S stand for susceptible population and infection is I = If + Ic. Hence If is represent the influenza infection rate and Ic for COVID infection rate. Due to the KSI2 term, squaring on I gives a double rate of infection; furthermore, the term KSI gives a single infection rate. This type of interaction will be considered in this research work. For numerical simulation, the non-standard finite difference (NSFD) scheme is used for deterministic model. Numerous nonlinear problems can be effectively and practi- cally solved numerically with this method. For instance, many scholars have used this technique to numerically simulate a variety of effusive systems [23, 24]. However, [25] discussed COVID-19 dynamics using fractional numerical simulation. One of the best methods for estimating solutions to fractional-order models among those discussed is the NSFD method; for further details, see [26, 27]. This study introduces a new vaccination-based deterministic co-infection approach for COVID-19 and influenza. The purpose of this study is to examine the dynamic and stabi- lization characteristics of the model equations as well as the connection between influenza prevalence and the COVID-19 vaccine. To discuss the dynamics of co-infection under convex incidence rate, we examine a mathematical SEIcIfVRC model. S stands for sus- ceptible, E for exposed, If for influenza infected, Ic for covid infected, V for vaccinated, and R for recovered, and C for combined infection compartment individual. These seven equations comprised the model. To determine the most vulnerable variables that raise or lower the fundamental reproduction number R0, sensitivity analysis is performed. The organization of the work is as follows: model formulation is discussed in Section 2. Mathe- matical results, local and global stability are discussed in Section 3. Numerical simulation through the NSFD scheme and the stochastic model are in Section 4, and the last section is a summarization of the results. 2. Model Formulation In this section, we present the dynamic of the co-infection of influenza and COVID-19 in the following model form. Basically the infection is divided into two way influenza If and COVID Ic, the whole population is represent by M which is equal to the sum of all compartments, M = S + E + If + Ic + C + V + R. We formulate the model using the R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 4 of 27 following equations: dS dt = b−KSI(βI + 1)− d0S + rR, dE dt = KSI(βI + 1)− (d0 + η + ϕ)E , dIf dt = ϕE − (γ1 + vf + d0 + c1)If , dIc dt = ηE − (γ2 + vc + d0 + c2)Ic, dC dt = c1If + c2Ic − (d0 + v3)C, dV dt = vfIf + vcIc + v3C − (γ3 + d0)V, dR dt = γ1If + γ2Ic + γ3V − (r + d0)R. (1) Here is the model flow chart which shows the dynamics of the flow with respect to time of different classes of the co-infection. The variables used in the system (1) with their physical meaning are, new recruitment 1. Flow chart S E R V If Ic C b Figure 1: Flow chart of parameters of the model (1). represent by b, death rate is represented by d0, the co-infection rate from influenza is c1 and COVID is c2. vf is vaccination rate of influenza, and vc is vaccination rate of covid-19. The recovery rate from influenza is represented by γ1, the recovery rate from COVID is by γ2 and recovery from vaccination is represented by γ3. K stands for infection constant, the influenza infection rate is represented by ϕ, the COVID infection rate by η, and the co-infection rate by β. The model recovered individuals lose immunity is represent by r. R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 5 of 27 3. Mathematical Analysis of Deterministic model (1) In this part of our work, we present the well posedness of the model, and feasible re- gion where solution exist. Furthermore, the uniqueness and existence are discussed. The equilibrium points are calculated for the model, disease free as well as endemic equilibrium points of the deterministic model (1). Lemma 1. The deterministic model (1) have invariant manifold of plane S + E + If + Ic + C + V +R = b d0 , which exist in first quadrant. Proof. Sum all equation of the system (1), and let suppose M represent the whole population of the community, then we have M = S + E + If + Ic + V +R+ C, which can be calculated as, dM dt = b− d0S − d0E − d0If − d0Ic − d0C − d0V − d0R = b− d0(S + E + If + Ic + C + V +R). By simple calculation, we get dM dt = b− d0M (2) Hence, the solution M = b d0 , is obtained from (2) Also, the general solution for M(t0) ≥ 0 of (2) is M = 1 d0 [b− (b− d0M(t0))e −d0(t−t0)]. Thus lim t→∞ M = b d0 . Hence, the proof is complete. This above result shows that the system (1) is mathematically well-posed. The possible region where solution lies is the following, Φ = (S, E , If , Ic,V,R, C) : 0 ≤ S, E , If , Ic,V,R, C,S + E + If + Ic + V +R+ C ≤ b d0 . R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 6 of 27 3.1. Existence and Uniqueness of Deterministic model (1) We search for the existence and uniqueness of the system (1). To do this, we supply the model (1) in the way that is dictated by S ′ = f1 (t,S, E , If , Ic, C,V,R, ) , E ′ = f2 (t,S, E , If , Ic, C,V,R, ) , I ′ f = f3 (t,S, E , If , Ic, C,V,R, ) , I ′ c = f4 (t,S, E , If , Ic, C,V,R, ) , V ′ = f5 (t,S, E , If , Ic, C,V,R, ) , R′ = f6 (t,S, E , If , Ic, C,V,R, ) , C′ = f7 (t,S, E , If , Ic, C,V,R, ) . (3) Norm is define as ||χ||∞ = sup t∈x . (4) Let for some positive constant h1, h2, h3, h4, h5, h6 and h7, where t ∈ [0, T ], we assume that the each compartment of deterministic model (1) is bounded in [0, T ] ||S||∞ < h1, ||E||∞ < h2, ||If ||∞ < h3, ||Ic||∞ < h4, ||V||∞ < h5, ||R||∞ < h6, ||C||∞ < h7, and k = h1 + h2. To investigated the bondness of model 1, we have |g1 (t,S, E , If , Ic,V,R, C) | =|b−KSI(βI + 1)− d0S + rR|, ≤b+ |S|K|I|(β|I|+ 1) + d0|S|+ r|R|, ≤b+ sup X |S|K sup X |I|(β sup X |I|+ 1) + d0 sup X |S|+ r sup X |R|, ≤b+ ||S||∞K||I||∞(β||I||∞ + 1) + d0||S||∞ + r||R||∞, ≤b+ h1Kk(βk + 1) + d0h1 + rh6 < ∞. (5) R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 7 of 27 Where M is equal to t, and belong to [0, T ]. Using the same procedure, we have |g2 (t,S, E , If , Ic,V,R, C) | =|KSI(βI + 1)− (d0 + η + ϕ)E|, ≤Kh1k(βk + 1) + (d0 + η + ϕ)h2 < ∞. (6) |g3 (t,S, E , If , Ic,V,R, C) | =|ϕE − (γ1 + vf + d0 + c1)If |, ≤ϕh2 + (γ1 + vf + d0 + c1)h3 < ∞. (7) |g4 (t,S, E , If , Ic,V,R, C) | =|ηE − (γ2 + vc + d0 + c2)Ic|, ≤ηh2 + (γ2 + vc + d0 + c2)h4 < ∞. (8) |g5 (t,S, E , If , Ic,V,R, C) | =|vfIf + vcIc − γ3V|, ≤vfh3 + vch4 + γ3h5 < ∞. (9) |g6 (t,S, E , If , Ic,V,R, C) | =|γ1If + γ2Ic + γ3V − (r + d0)R|, ≤γ1h3 + γ2h4 + γ3h5 + (r + d0)h6 < ∞. (10) |g7 (t,S, E , If , Ic,V,R, C) | =|c1If + c2Ic|, ≤c1h3 + c2h4 < ∞. (11) Thus S, E , If , Ic,V,R and C all are bounded in the given region, and there exist h1, h2, h3, h4, h5 and h6 such that sup X |g1 (t,S, E , If , Ic,V,R, C) | < h1, sup X |g2 (t,S, E , If , Ic,V,R, C) | < h2, sup X |g3 (t,S, E , If , Ic,V,R, C) | < h3, sup X |g4 (t,S, E , If , Ic,V,R, C) | < h4, sup X |g5 (t,S, E , If , Ic,V,R, C) | < h5, sup X |g6 (t,S, E , If , Ic,V,R, C) | < h6. sup X |g7 (t,S, E , If , Ic,V,R, C) | < h7. Now, we are going to prove that, |g1 (t,S2, E , If , Ic,V,R, C)− g1 (t,S1, E , If , Ic,V,R, C) | = |b−KS2I(βI + 1)− d0S2 + rR− (b−KS1I(βI + 1)− d0S1 + rR)| < (KI(βI + 1) + d0)|S2 − S1| < p1|S2 − S1|. (12) R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 8 of 27 Where p1 = (KI(βI + 1) + d0). Follow the similar method, we get |g2 (t,S, E2, If , Ic,V,R, C)− g2 (t,S, E1, If , Ic,V,R, C) | < p2|E2 − E1|, |g3 ( t,S, E , Ig2 , I⌋,V,R, C ) − f3 (t,S, E , Ig1 , Ic,V,R, C) | < p3|Ig2 − Ig1 |, |g4 (t,S, E , If , Ic2 ,V,R, C)− f4 (t,S, E , If , Ic1 ,V,R, C) | < p4|Ic2 − Ic1 |, |g5 (t,S, E , If , Ic,V,R, C2)− f5 (t,S, E , If , Ic,V,R, C1) | < p5|C2 − C1|. |g6 (t,S, E , If , Ic,V2,R, C)− f6 (t,S, E , If , Ic,V2,R, C) | < p6|V2 − V1|, |g7 (t,S, E , If , Ic,V,R2, C)− f7 (t,S, E , If , Ic,V,R1, C) | < p7|R2 −R1|. Where p2 = d0 + ϕ+ η, p3 = γ1 + vf + d0 + c1, p4 = γ2 + vc + d0 + c2, p5 = d0 + v3 p6 = γ3 + d0, p7 = r + d0. This uniqueness and existence show that the system (1) has a unique set of solutions. 3.2. Disease Free Equilibrium Points DFEP: P0 represent the DFEP of model (1), which is equal to (S0, 0, 0, 0, 0, 0, 0). May be write in the following form, P0 = ( b d0 , 0, 0, 0, 0, 0, 0 ) . R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 9 of 27 3.3. Endemic Equilibrium points EEP We obtained EEP for our model (1) as S∗ = b(r + d0) + r [ (γ1 + vf )I∗ f + (γ2 + vc)I∗ c ] (r + d0)(KI∗(βI∗ + 1) + d0) , E∗ = KI∗(βI∗ + 1)(b(r + d0) + r[(γ1 + vf )I∗ f + (γ2 + vc)I∗ c ]) (r + d0)(d0 + ϕ+ η) , I∗ f = ϕKI∗(βI∗ + 1)(b(r + d0) + r[(γ1 + vf )I∗ f + (γ2 + vc)I∗ c ]) (r + d0)(d0 + ϕ+ η)(γ1 + vf + d0 + c1) , I∗ c = ηKI∗(βI∗ + 1)(b(r + d0) + r[(γ1 + vf )I∗ f + (γ2 + vc)I∗ c ]) (r + d0)(d0 + ϕ+ η)(γ2 + vc + d0 + c2) , V∗ = vfI∗ f + vcI∗ c γ3 , R∗ = (γ1 + vf )I∗ f + (γ2 + vc)I∗ c r + d0 . (13) 3.4. Basic Reproduction number R0 In epidemiology, the term R0, or basic reproduction number, refers to the way in which infections are controlled and propagated. Information regarding the disease’s prevalence in the community and the best ways to protect the local people from the deadly virus are provided by R0. We utilize the next generation matrix notion to determine R0. Let Φ = (E , If , Ic), then form system (1), we have dΦ dt = X − Y, such that, X =  KSI(βI + 1) 0 0  and Y =  (d0 + η + ϕ)E −ϕE + (γ1 + vf + d0 + c1)If −ηE + (γ2 + vc + d0 + c2)Ic  . At DFEP the Jacobian of X is, we have X =  0 kS0 kS0 0 0 0 0 0 0  R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 10 of 27 where at DFEP the Jacobian of X is, we have Y =  (d0 + η + ϕ) 0 0 −ϕ γ1 + vf + d0 + c1 0 −η 0 γ2 + vc + d0 + c2  . XY−1 =  0 KS0 KS0 0 0 0 0 0 0  q(γ1 + vf + d0 + c1) 0 0 ϕ(γ2 + vc + d0 + c2) q(γ2 + vc + d0 + c2) 0 η(γ1 + vf + d0 + c1) 0 q(γ1 + vf + d0 + c1)  (q(γ1 + vf + d0 + c1)(γ2 + vc + d0 + c2) . Assume that q = d0 + η + ϕ, where XY−1 have the following positive radius, which is known as reproduction rate of influenza and covid respectively, R01 = Kb d0(d0 + ϕ+ η)(γ1 + vf + d0 + c1) R02 = Kbϕ d0(d0 + ϕ+ η)(γ2 + vc + d0 + c2) . May some one consider as, R0 = max [R01,R02] . (14) 4. Stability Analysis In this part of the manuscript consist of stability analysis of the deterministic model (1). Actually, we check here local and global stability of the system at disease free and endemic equilibrium points. 4.1. Local Stability Local stability analysis is investigate though applying the Jacobian matrix of the de- terministic model (1. As in the next theorem, we get condition under which the system is locally asymptotically stable. Theorem 1. “ When effective basic reproduction number R01 < 1 or R02 < 1, the DFEP “P” is stable locally asymptotically. ” Proof. At disease free equilibrium state P0 = ( b d0 , 0, 0, 0, 0, 0, 0 ) the Jacobian matrix of deterministic model (1 as, T 0 = −d0 0 0 0 0 0 r 0 −(d0 + ϕ+ η) 0 0 0 0 0 0 ϕ −(γ1 + vf + d0 + c1) 0 0 0 0 0 η 0 −(γ2 + vc + d0 + c2) 0 0 0 0 0 c1 c2 0 0 0 0 0 vf vc −γ3 0 0 0 0 γ1 γ2 γ3 −(r − d0) 0  R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 11 of 27 and |χ− T 0(P)| = 0. Hence χ1 = −d0, χ2 = −(d0 + ϕ+ η), χ3 = −(γ2 + vc + d0 + c2), χ4 = −(γ1 + vf + d0 + c1), χ5 = −(r + d0), χ6 = (γ2 + vc + d0 + c2)(R01 − 1), χ7 = (γ1 + vf + d0 + c1)(R02 − 1). From last two equation, we see that χ6 is negative whenever R01 < 1, χ7 has condition R02 < 1, and the rest are also strictly negative which follow the conclusion. 4.2. Global Stability Global stability of the deterministic model is highlight through a function “F”, known is Laypaunov function. By using [25] at disease free equilibrium sate “P”, we investigate the stability analysis in next theorem. Theorem 2. “When effective basic reproduction number R01 < 1 or R02 < 1, the DFEP “P” is globally asymptotically stable.” Proof. Let define the Laypaunov function F as, F = h1If + h1Ic, (15) where h1 and h2 are some positive constants. Since If , Ic > 0, then F . Also whenever If = Ic = 0 then F = 0. We need to show that Ḟ < 0, for this follow Now, our goal is Ḟ < 0. dF dt = h1 dIf dt + h2 dIc dt . (16) Putting values of If and Ic from (1) in (16), we get dF dt = h1 (ϕE − (γ1 + vf + d0 + c1)If ) + h2 (ηE − (γ2 + vc + d0 + c2)Ic) . (17) Where h1 = kb d0(d0 + ϕ+ η) also same as h2 = kb d0(d0 + ϕ+ η) R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 12 of 27 . From (17), we have dF dt = (R01 − 1) (γ1 + d0 + vf + c1)If + (R02 − 1) (γ2 + d0 + vc + c2)Ic. (18) From the principle of “LaSalle’s invariant”, dF dt is equal to zero whenever If = Ic = 0. Also, we see that if R01 < 1 and R02 < 1, then dF dt must be less than zero, which follows the conclusion. 5. Sensitivity Analysis of the Deterministic Model (1) In order to account for the factors that significantly affect the basic reproduction number R0, we present the sensitivity analysis in this section. To figure out the significance of the many factors influencing the incidence and spread of disease, sensitivity analysis is advised. Controlling the system (1) variables is essential to achieving R01 < 1 or R02 < 1, in addition to to prevent the spread of illnesses. The formula can be used to determine the sensitivity index, which represents the proportion of changes in a variable to a parameter’s change. In the following 1, we present the sensitivity index of each parameters, involve in the basic reproduction numbers R01 < 1. The formula used for determining the the sensitivity index is h[y] = X R0 ∂R0 ∂X , where X mean parameter of R0. Parameters Sensitivity Index h(b) 1 h(K) 1 h(d0) -0.81 h(ϕ) 0.69 h(η) -0.21 h(γ1) -0.78 h(γ2) -0.50 h(c1) -0.78 h(c2) -0.13 h(vf ) -0.78 h(vc) -0.67 Table 1: R01 and R02 sensitivity index relative to system (1) variables.) It is clear from Table (1) and Fig. (2) that the variables of birth rate b, constant of infection and influenza infection rate are positive so there increase while effect increasing in the value of R01 and R02. On the other hand, the values of d0, η, ϕ, γ1, γ2, c1, c2, vf , and vc are negative so their increased which follow decrease in R01 and R02. R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 13 of 27 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 b K d_0 Ø � �_1 �_1 c1 c2 vf vc V a lu e s R_0 Sensitivity index Figure 2: R01 and R02 sensitivity index relative to system (1) variables. 6. Numerical simulation through NSFD of deterministic model By using non-standard finite difference scheme (NSFD) [28? , 29], we plot the numer- ical simulation of the deterministic model (1) in this section. First, the model equations are written in the following way, as consider the first equation of model (1), dS dt = b−KSI(βI + 1)− d0S + rR. (19) From the non-standard finite difference method, we decomposed as Sj+1 − Sj h = b−KSjIj(βIj + 1)− d0Sj + rRj . (20) Like (20), we decomposed other equations of model (1) by using none-standard finite difference method Ej+1 = Ej + h ( KSjIj(βIj + 1)− (d0 + η + ϕ)Ej ) , If(j+1) = I1(j) + h ( ϕEj − (γ1 + d0 + vf + c1)I1j ) , Ic(j+1) = 2(j) + h ( ηEj − (γ2 + vc + d0 + c2)Ij ) , Vj+1 = Vj + h ( vfI1j + vcI2j + v3Cj − γ3Vj ) , Rj+1 = Rj + h ( γ1I1j + γ2I2j + γ3Vj − (r + d0)Rj ) , Cj+1 = Cj + h ( c1I1j + c2I2j − (v3 + d0)Cj ) . R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 14 of 27 Variable numerical values Variable numerical S 0.4540197 Million ϕ 0.12 E 0.0030000 Million γ1 0.015 If 0.0001455 Million γ2 0.0667 Ic 0.1428015 Million c1 0.1 V 0.0141480 Million c2 0.005 R 0.2700000 Million η 0.3425 C 0.0170000 Million b 0.002500 d0 0.01 β 0.12 K 0.1175 r 0.000167 vf 0.003 vc 0.0015 γ3 0.036 Table 2: The parameters of the system are described and specified, with approximate real values (1). t 0 5 10 15 20 25 30 35 40 S u s c e p ti b le c la s s 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 Region 1 Region 2 Figure 3: Graph of numerical solutions using the NSFD technique for the deterministic model (1) susceptible compartment. By using NSFD scheme, we simulate the system (1) utilizing the numerical values listed in the 2 taking the values from [4]. It is apparent in Fig.3 to Fig.9, the findings that the recovered populations grow slowly and eventually stable at the top, while the susceptible and exposed populations progressively decline toward zero throughout the first 50 days, by using a Table.2 and graph the numerical data using the NSFD method. As predicted, the covid-19 infected and influenza infected populations steadily declines after one week. Initially co-infected populations increased significantly in first two week, but after two week it start decline and become stable in 50 days. Also, when infection decreased the vaccinated individuals and recovered become stabilize. These findings support Theorem 3.1, which states that the disease dies out whenever the basic reproduction number less then 1. Additionally, it indicates that the SARS-Cov-2 virus will be eliminate from the population in 50 days. R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 15 of 27 t 0 5 10 15 20 25 30 35 40 E x p o s e d c la s s ×10-3 2 4 6 8 10 12 14 Region 1 Region 2 Figure 4: Graph of numerical solutions using the NSFD technique for the deterministic model (1) exposed compartment. t 0 5 10 15 20 25 30 35 40 In fl u e n z a i n fe c te d ×10-3 0 1 2 3 4 5 6 7 8 Region 1 Region 2 Figure 5: Graph of numerical solutions using the NSFD technique for the deterministic model (1) influenza infected compartment. 7. Stochastic model Here, we introduce the influence of environmental white noise to transform the pre- viously proposed deterministic system (1) into a stochastic model. We accomplish this by introducing nonlinear perturbations into each of the system’s equations, in this case merely the rate of each class, as illustrated for each class below. S(t) : −β −→ −β + (χ11S + χ12)dA1(t), E(t) : −ϕ −→ −ϕ+ (Π21S + χ22)dA2(t), If (t) : −c1 −→ −c1 + (χ31S + χ32)dA3(t), R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 16 of 27 t 0 5 10 15 20 25 30 35 40 C O V ID i n fe c te d 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Region 1 Region 2 Figure 6: Graph of numerical solutions using the NSFD technique for the deterministic model COVID (1) infected compartment. t 0 5 10 15 20 25 30 35 40 C O -i n fe c te d 0.015 0.02 0.025 0.03 0.035 0.04 Region 1 Region 2 Figure 7: Graph of numerical solutions using the NSFD technique for the deterministic model (1) co-infection compartment. Ic(t) : −c2 −→ −c2 + (χ41S + χ42)dA4(t), C(t) : −v3 −→ −v3 + (χ51S + χ52)dA5(t), V(t) : γ3 −→ −γ3 + (χ61S + χ62)dA6(t), R(t) : −d0 −→ −d0 + (χ71S + χ72)dA7(t). The modified stochastic model (1) is thus represented by the system of equations that follows: dS = [b−KSI(βI + 1)− d0S + rR]dt+ (χ11S + χ12)A1(t), dE = [KSI(βI + 1)− (d0 + η + ϕ)E ]dt+ (χ21S + χ22)A2(t), R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 17 of 27 t 0 5 10 15 20 25 30 35 40 R e c o v e re d 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 Region 1 Region 2 Figure 8: Graph of numerical solutions using the NSFD technique for the deterministic model (1). recovered compartment t 0 5 10 15 20 25 30 35 40 V a c c in a te d 0.004 0.006 0.008 0.01 0.012 0.014 0.016 Region 1 Region 2 Figure 9: Graph of numerical solutions using the NSFD technique for the deterministic model (1) vaccinated compartment. dIf = [ϕE − (γ1 + vf + d0 + c1)If ]dt+ (χ31S + χ32)A3(t), dIc = [ηE − (γ2 + vc + d0 + c2)Ic]dt+ (χ41S + χ42)A4(t), dC = [c1If + c2Ic − (d0 + v3)C]dt+ (χ51S + χ52)A5(t), dV = [vfIf + vcIc + v3C − (γ3 + d0)V]dt+ (χ61S + χ62)A6(t), dR = [γ1If + γ2Ic + γ3V − (r + d0)R]dt+ (χ71S + χ72)A7(t). R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 18 of 27 8. Numerical simulation of stochastic model The following equation is used to determine the trajectories or estimated solutions of a stochastic model equations using the Euler-Maruyama method: χti+1 = χti + α(ti, χti) + β(ti, χti)∆Ai. (21) Where i = 0, 1, ..., n − 1. It is required to understand how to compute ∆Ai in order to implement the procedure computationally. The distribution of the differences ∆Ai, i = 0, 1, ..., n− 1 has the same value,∆Ai ∼ M(0,∆t), because the partition is composed of equal intervals. Consider a random variable η which has a normal η ∼ M(0, 1). Then√ ∆tη1 has zero mean and variance ∆t, indicating a normal distribution; that is √ ∆tη1 ∼ M(0,∆t). For our proposed model, to implement the Euler-Maruyama method algorithm similarly (21), we must perform the appropriate separating of the system of stochastic differential equations (21), which is given by Sti+1 = Sti + [b−KSI(βI + 1)− d0S + rR] ∆t+ √ ∆tη1, Eti+1 = Eti + [KSI(βI + 1)− (d0 + η + ϕ)E ] ∆t+ √ ∆tη2, Ift i+1 = Ifti + [ϕE − (γ1 + vf + d0 + c1)If ] ∆t+ √ ∆tη3, Ict i+1 = Icti + [ηE − (γ2 + vc + d0 + c2)Ic]∆t+ √ ∆tη4, Cti+1 = Cti + [c1If + c2Ic − (d0 + v3)C]∆t+ √ ∆tη5, Vti+1 = Vti + [vfIf + vcIc + v3C − (γ3 + d0)V]∆t+ √ ∆tη6, Rti+1 = Rti + [γ1If + γ2Ic + γ3V − (r + d0)R]∆t+ √ ∆tη7. (22) t 0 5 10 15 20 25 30 35 40 S u s c e p ti b le c la s s s to c h a s ti c m o d e l 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 Figure 10: Graph of numerical solutions for the stochastic model (21) susceptible compartment by using Euler- Maruyama. R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 19 of 27 t 0 5 10 15 20 25 30 35 40 E x p o s e d c la s s s to c h a s ti c m o d e l 0 0.05 0.1 0.15 0.2 0.25 0.3 Figure 11: Graph of numerical solutions for the stochastic model (21) exposed compartment by using Euler- Maruyama.. t 0 5 10 15 20 25 30 35 40 In fl u e n z a i n fe c te d s to c h a s ti c m o d e l 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 Figure 12: Graph of numerical solutions for the stochastic model (21) influenza compartment by using Euler- Maruyama. Using the algorithm mentioned above and data from Table.2, we next apply the model and get the results in Fig.(10-16). This confirms the analytical result done in Theorems 5.1 and 5.2, as well as the effect of vaccination and the degree of white noise on the movement of the novel coronavirus disease. Mortality analysis clearly shows that, in the case that RS 1 < 1, there will always be susceptible and recovered people and the infected individual disappears, as above plot illustrates. However, if RS 1 > 1, there will still be sick individuals and the disease will continue, as mention above figures illustrates. Furthermore, Figs (15) illustrate how vaccination and the level of white noise affect the dynamics of the compartmental population. It is evident that immunization and the white noise degree variable are important factors in both situations. We found that the amount R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 20 of 27 t 0 5 10 15 20 25 30 35 40 C O V ID i n fe c te d s to c h a s ti c m o d e l 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Figure 13: Graph of numerical solutions for the stochastic model (21) covid compartment by using Euler- Maruyama. t 0 5 10 15 20 25 30 35 40 C O -i n fe c te d s to c h a s ti c m o d e l 0 0.05 0.1 0.15 0.2 0.25 0.3 Figure 14: Graph of numerical solutions for the stochastic model (21) co-infected compartment by using Euler- Maruyama. of white noise and vaccination had a significant impact on the spread of new coronavirus diseases. It should be highlighted that, as Figs. 6 and 7 demonstrate, increasing the values of these two parameters would accelerate the extinction of the disease. Consequently, as the values of vf and vc rise, the number of susceptible and infected individuals falls while the numbers of recovered population rises. Similarly, as shown in above figures, in the situation of disease persistence, vf is directly proportional to the number of recovered individuals while inversely related to the amount of susceptible and infected individuals. However, vf the parameter is directly related to the number of susceptible and infected persons and inversely proportional to the amount of infected individuals. To see the figures (17-23) the dynamics of the stochastic and deterministic follow real data, which R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 21 of 27 t 0 5 10 15 20 25 30 35 40 V a c c in a te d s to c h a s ti c m o d e l 0 0.005 0.01 0.015 0.02 0.025 0.03 Figure 15: Graph of numerical solutions for the stochastic model (21) vaccination compartment by using Euler-Maruyama. t 0 5 10 15 20 25 30 35 40 R e c o v e re d s to c h a s ti c m o d e l 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Figure 16: Graph of numerical solutions for the stochastic model (21) recovered compartment by using Euler- Maruyama. show accuracy of our model. 9. Conclusion Cases of co-infection with SARS-CoV-2 and influenza have been investigated in recent studies [30–34]. Mathematical models can provide valuable insights into the dynamics of mutual infection within a host. In this work, we developed and analyzed a system of differ- ential equations (DDEs) and stochastic differential equations (SDEs) to characterize the progression of influenza and SARS-CoV-2 co-infection. We established the system’s funda- mental properties, including boundedness and non-negativity. Furthermore, we identified R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 22 of 27 Time (Days) 0 5 10 15 20 25 30 35 40 S u s c e p ti b le 0.25 0.3 0.35 0.4 0.45 0.5 Susceptible Stochastic Deterministic Real Data Figure 17: Comparison of susceptible compartment. Time (Days) 0 5 10 15 20 25 30 35 40 E x p o s e d 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 Exposed Stochastic Deterministic Real Data Figure 18: Comparison of exposed compartment. the system’s equilibria and examined their existence and stability, locally and globally. Our numerical simulations were consistent with the theoretical results, reinforcing the validity of our analysis. As part of future work, we aim to enhance our model by integrating neural networks and artificial intelligence (AI) techniques [35]. These tools will enable more accurate prediction, parameter estimation, and real-time adaptation of the model, ultimately ad- vancing our understanding of co-infection dynamics and improving decision-making in public health interventions. Declarations All authors have read and approved the final manuscript. Authors contributions R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 23 of 27 Time (Days) 0 5 10 15 20 25 30 35 40 In fl u e n z a I n fe c te d ×10-3 0 1 2 3 4 5 6 7 8 Influenza Infected Stochastic Deterministic Real Data Figure 19: Comparison of influenza infected compartment. Time (Days) 0 5 10 15 20 25 30 35 40 C O V ID I n fe c te d 0 0.05 0.1 0.15 COVID Infected Stochastic Deterministic Real Data Figure 20: Comparison of COVID infected compartment. Conceptualization: [Author A, Author B] Methodology: [Author A] Formal Analysis: [Author A] Investigation: [Author A, Author B] Data Correction: [Author B] Writing Original Draft: [Author A] Writing Review & Editing: [Author A, Author B, ] Visualization: [ Author B] Supervision: [Author B] Project Administration: [Author B] Funding Acquisition: [Author B] R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 24 of 27 Time (Days) 0 5 10 15 20 25 30 35 40 C o -I n fe c te d 0.015 0.02 0.025 0.03 0.035 0.04 Co-Infected Stochastic Deterministic Real Data Figure 21: Comparison of co-infection compartment. Time (Days) 0 5 10 15 20 25 30 35 40 V a c c in a te d ×10-3 2 4 6 8 10 12 14 16 Vaccinated Stochastic Deterministic Real Data Figure 22: Comparison of vaccination compartment. • Acknowledgements The authors would like to acknowledge the support of Prince Sultan University for paying the Article Processing Charges (APC) of this publication. R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 25 of 27 Time (Days) 0 5 10 15 20 25 30 35 40 R e c o v e re d 0.26 0.28 0.3 0.32 0.34 0.36 0.38 Recovered Stochastic Deterministic Real Data Figure 23: Comparison of recovered compartment. References [1] World Health Organization. WHO Timeline–COVID-19, 2021. [2] World Health Organization (WHO). WHOwww.who.int/emergencies/diseases/novel- coronavirus-2019/interactivetimeline/category-information, 2022. [3] World Health Organization (WHO). WHOwww.worldometers.info/coronavirus/nov2022, 2022. [4] World Health Organization (WHO). WHOwww.who.int/emergencies/diseases/novel- coronavirus-2019, 2022. [5] Rohit C Khanna, Maria Vittoria Cicinelli, Suzanne S Gilbert, Santosh G Honavar, and Gudlavalleti VS Murthy. Covid-19 pandemic: Lessons learned and future directions. Indian journal of ophthalmology, 68(5):703–710, 2020. [6] Sean Wei Xiang Ong, Yian Kim Tan, Po Ying Chia, Tau Hong Lee, Oon Tek Ng, Michelle Su Yen Wong, and Kalisvar Marimuthu. Air, surface environmental, and personal protective equipment contamination by severe acute respiratory syndrome coronavirus 2 (sars-cov-2) from a symptomatic patient. Jama, 323(16):1610–1612, 2020. [7] Rutu Karia, Ishita Gupta, Harshwardhan Khandait, Ashima Yadav, and Anmol Ya- dav. Covid-19 and its modes of transmission. SN comprehensive clinical medicine, 2(10):1798–1801, 2020. [8] Michael O Adeniyi, Segun I Oke, Matthew I Ekum, Temitope Benson, and Matthew O Adewole. Assessing the impact of public compliance on the use of non-pharmaceutical intervention with cost-effectiveness analysis on the transmission dynamics of covid- 19: Insight from mathematical modeling. Modeling, control and drug development for COVID-19 outbreak prevention, pages 579–618, 2022. [9] Saif Ullah and Muhammad Altaf Khan. Modeling the impact of non-pharmaceutical interventions on the dynamics of novel coronavirus with optimal control analysis with R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 26 of 27 a case study. Chaos, Solitons & Fractals, 139:110075, 2020. [10] Legesse Lemecha Obsu and Shiferaw Feyissa Balcha. Optimal control strategies for the transmission risk of covid-19. Journal of biological dynamics, 14(1):590–607, 2020. [11] Zhong-Hua Shen, Yu-Ming Chu, Muhammad Altaf Khan, Shabbir Muhammad, Omar A Al-Hartomy, and M Higazy. Mathematical modeling and optimal control of the covid-19 dynamics. Results in Physics, 31:105028, 2021. [12] Jiahui Chen, Rui Wang, Nancy Benovich Gilby, and Guo-Wei Wei. Omicron variant (b. 1.1. 529): infectivity, vaccine breakthrough, and antibody resistance. Journal of chemical information and modeling, 62(2):412–422, 2022. [13] Petra Mlcochova, Steven A Kemp, Mahesh Shanker Dhar, Guido Papa, Bo Meng, Isabella ATM Ferreira, Rawlings Datir, Dami A Collier, Anna Albecka, Sujeet Singh, et al. Sars-cov-2 b. 1.617. 2 delta variant replication and immune evasion. Nature, 599(7883):114–119, 2021. [14] Mohamed Derouich and Abdesslam Boutayeb. An avian influenza mathematical model. Applied mathematical sciences, 2(36):1749–1760, 2008. [15] World Health Organization (WHO). WHO statistics, 2022. [16] Timo Smieszek, Gianrocco Lazzari, and Marcel Salathé. Assessing the dynamics and control of droplet-and aerosol-transmitted influenza using an indoor positioning system. Scientific reports, 9(1):2185, 2019. [17] Mayowa M Ojo, Temitope O Benson, Olumuyiwa James Peter, and Emile Franc Doungmo Goufo. Nonlinear optimal control strategies for a mathematical model of covid-19 and influenza co-infection. Physica A: Statistical Mechanics and its Applications, 607:128173, 2022. [18] Lei Bai, Yongliang Zhao, Jiazhen Dong, Simeng Liang, Ming Guo, Xinjin Liu, Xin Wang, Zhixiang Huang, Xiaoyi Sun, Zhen Zhang, et al. Coinfection with influenza a virus enhances sars-cov-2 infectivity. Cell research, 31(4):395–403, 2021. [19] Xuan Xiang, Zi-hao Wang, Lin-lin Ye, Xin-liang He, Xiao-shan Wei, Yan-ling Ma, Hui Li, Long Chen, Xiao-rong Wang, and Qiong Zhou. Co-infection of sars-cov-2 and influenza a virus: a case series and fast review. Current Medical Science, 41(1):51–57, 2021. [20] John Paget, Saverio Caini, Ben Cowling, Susanna Esposito, Ann R Falsey, Angela Gentile, Jan Kyncl, C MacIntyre, Richard Pitman, and Bruno Lina. The impact of influenza vaccination on the covid-19 pandemic? evidence and lessons for public health policies. Vaccine, 38(42):6485, 2020. [21] Deus Thindwa, Maria Garcia Quesada, Yang Liu, Julia Bennett, Cheryl Cohen, Maria Deloria Knoll, Anne von Gottberg, Kyla Hayford, and Stefan Flasche. Use of seasonal influenza and pneumococcal polysaccharide vaccines in older adults to reduce covid-19 mortality. Vaccine, 38(34):5398, 2020. [22] Suman Bhowmick, Igor M Sokolov, and Hartmut HK Lentz. Decoding the dou- ble trouble: A mathematical modelling of co-infection dynamics of sars-cov-2 and influenza-like illness. Biosystems, 224:104827, 2023. [23] Ebrahem A Algehyne and Rahim ud Din. On global dynamics of covid-19 by using sqir type model under non-linear saturated incidence rate. Alexandria Engineering R. Ud Din, M. S. Arif / Eur. J. Pure Appl. Math, 18 (2) (2025), 6005 27 of 27 Journal, 60(1):393–399, 2021. [24] Kamal Shah, Thabet Abdeljawad, and Rahim 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. [25] Peter M Manning and Gary F Margrave. Introduction to non-standard finite- difference modelling. CREWES Res. Rep, 18:1–10, 2006. [26] Muhammad Arfan, Maha MA Lashin, Pongsakorn Sunthrayuth, Kamal Shah, Aman Ullah, Kulpash Iskakova, MR Gorji, and Thabet Abdeljawad. On nonlinear dynamics of covid-19 disease model corresponding to nonsingular fractional order derivative. Medical & Biological Engineering & Computing, 60(11):3169–3185, 2022. [27] Muhammad Arfan, Kamal Shah, and Aman Ullah. Some theoretical and compu- tation results about covid-19 by using a fractional-order mathematical model. In Fractional-Order Modeling of Dynamic Systems with Applications in Optimization, Signal Processing and Control, pages 37–68. Elsevier, 2022. [28] Zhuo-Jia Fu, Zhuo-Chao Tang, Hai-Tao Zhao, Po-Wei Li, and Timon Rabczuk. Nu- merical solutions of the coupled unsteady nonlinear convection-diffusion equations based on generalized finite difference method. The European Physical Journal Plus, 134(6):272, 2019. [29] Rahim Ud Din, Khalid Ali Khan, Ahmad Aloqaily, Nabil Mlaiki, and Hussam Alraba- iah. Using non-standard finite difference scheme to study classical and fractional order seivr model. Fractal and Fractional, 7(7):552, 2023. [30] Resat Ozaras, Rasim Cirpin, Arif Duran, Habibe Duman, Ozgur Arslan, Yasin Bak- can, Metin Kaya, Huseyin Mutlu, Leyla Isayeva, Fatih Kebanlı, et al. Influenza and covid-19 coinfection: report of six cases and review of the literature. Journal of medical virology, 92(11):2657–2665, 2020. [31] Qiang Ding, Panpan Lu, Yuhui Fan, Yujia Xia, and Mei Liu. The clinical charac- teristics of pneumonia patients coinfected with 2019 novel coronavirus and influenza virus in wuhan, china. Journal of medical virology, 92(9):1549–1555, 2020. [32] Xiaojuan Zhu, Yiyue Ge, Tao Wu, Kangchen Zhao, Yin Chen, Bin Wu, Fengcai Zhu, Baoli Zhu, and Lunbiao Cui. Co-infection with respiratory pathogens among covid-2019 cases. Virus research, 285:198005, 2020. [33] Mhammed Mediani, Abdeldjalil Slama, Ahmed Boudaoui, and Thabet Abdeljawad. Analysis of a stochastic seiuirr epidemic model incorporating the ornstein-uhlenbeck process. Heliyon, 10(16), 2024. [34] Israr Ahmad, Zeeshan Ali, Babar Khan, Kamal Shah, and Thabet Abdeljawad. Ex- ploring the dynamics of gumboro-salmonella co-infection with fractal fractional anal- ysis. Alexandria Engineering Journal, 117:472–489, 2025. [35] Hasib Khan, Jehad Alzabut, DK Almutairi, Haseena Gulzar, and Wafa Khalaf Alqurashi. Data analysis of fractal-fractional co-infection covid-tb model with the use of artificial intelligence. Fractals, 2025.