Global dynamics of IAV/SARS-CoV-2 coinfection model with eclipse phase and antibody immunity

1.   Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and influenza A virus (IAV) are two respiratory RNA viruses with high pandemic potential. SARS-CoV-2 causes the coronavirus disease 2019 (COVID-19). According to the update provided by the World Health Organization (WHO) on 2 October 2022 ^[1], over 593 million confirmed cases and over 6.4 million deaths were reported globally. Influenza viruses infect about 20% of the world's population in annual epidemics, resulting in 3-5 million severe illnesses and 290, 000–650, 000 deaths each year ^[2].

Both SARS-CoV-2 and IAV infect the uninfected epithelial cells of the host respiratory tract ^[3,4], and have analogous transmission ways. Moreover, they have common clinical manifestations including dyspnea, cough, fever, headache, rhinitis, myalgia and sore throat ^[5]. Viral shedding usually takes place 5 to 10 days in influenza, whereas it does 2 to 5 weeks in COVID-19 ^[5]. Acute respiratory distress is less common in influenza than in COVID-19 ^[5]. Deaths in influenza cases are less than 1%, while in cases of COVID-19 it ranges from 3 to 4% ^[5]. The precautionary measures implemented by governments to limit the transmission of SARS-CoV-2 can play a role in reducing the transmission of the IAV ^[6].

Eleven vaccines for COVID-19 were approved by WHO for emergency use. These include Novavax/Nuvaxovid, Bharat Biotech/Covaxin, CanSino/Convidecia, Pfizer/BioNTech/Comirnaty, Moderna/Spikevax, Serum Institute of India COVOVAX (Novavax formulation), Janssen (Johnson & Johnson)/Jcovden, Oxford/AstraZeneca/Vaxzevria, Serum Institute of India Covishield (Oxford/AstraZeneca formulation), Sinopharm (Beijing)/Covilo, and Sinovac/CoronaVac ^[7]. Currently, there are three types of influenza vaccines used worldwide: live attenuated influenza vaccine, inactivated influenza vaccine and recombinant hemagglutinin vaccine ^[8].

It was reported in ^[9] that, 94.2% of individuals with COVID-19 were also coinfected with several other microorganisms, such as fungi, bacteria and viruses. Important viral copathogens include the respiratory syncytial virus (RSV), human enterovirus (HEV), human rhinovirus (HRV), influenza A virus (IAV), influenza B virus (IBV), human metapneumovirus (HMPV), parainfluenza virus (PIV), human immunodeficiency virus (HIV), cytomegalovirus (CMV), dengue virus (DENV), Epstein Barr virus (EBV), hepatitis B virus (HBV) and other coronaviruses (COVs), among which the HRV, HEV and IAV are the most common copathogens ^[10]. Several coinfection cases of COVID-19 and influenza were reported in ^{[5,9,11,12,13]} (see also the review papers ^{[14,15,16,17,18]}). Lansbury et al. ^[14] presented a systematic review and meta-analysis that included 30 studies for evaluating coinfections among patients infected with COVID-19. It was reported that 7% of patients had a bacterial coinfection and 3% of patients had another respiratory virus, with RSV and IAV being among the most common coinfecting viruses. Dao et al. ^[15] conducted a systematic review and meta-analysis that included 54 publications and found that, 7% of COVID-19 patients are co-infected with influenza viruses. Most influenza co-infections were due to the IAV ^[15]. A respective study in Wuhan, China showed that the coinfection rate of IAV and SARS-CoV-2 was 49.8% during the outbreak period of COVID-19 ^[19]. Based on two separate studies presented in ^[11] and ^[12], COVID-19-influenza coinfection did not result in worse clinical outcomes ^[11]. In addition, this condition reduced the mortality rate among COVID-19-influenza coinfected patients. Coinfection with influenza virus in COVID-19 patients might render them less vulnerable to morbidities associated with COVID-19, and therefore, a better prognosis overall ^[12]. In ^[18], it is found that, although patients with IAV and SARS-CoV-2 coinfection did not experience longer hospital stays compared with those SARS COV-2 single-infection, they usually presented with a more severe clinical conditions. In an animal study ^[20], it was found that the disease severity is increased in hamsters with SARS-CoV-2 and IAV coinfection compared with those with SARS-CoV-2 mono-infection.

Viral interference phenomenon can appear in case of infections with multiple competitive respiratory viruses ^[21,22,23]. One virus may be able to suppress the growth of another virus ^[21,24]. In ^[22], it was reported that an H3N2 strain of IAV was inhibited by SARS-CoV-2 coinfection in the hamster model. Oishi et al. ^[23] used the golden Syrian hamster model, and found that, IAV interferes with SARS-CoV-2 replication in the lung, even more than one week after IAV clearance. Disease progression and outcome in SARS-CoV-2 infection are highly dependent on the host immune response, particularly in the elderly in whom immunosenescence may predispose to increased risk of coinfection ^[21].

1.1.   Mathematical models of within-host IAV and SARS-CoV-2 infections

Over the years, mathematical models have demonstrated their ability to provide useful insight to gain a further understanding of the dynamics and mechanisms of the viruses within a host level. These models may assist in the development of viral therapies and vaccines as well as the selection of appropriate therapeutic and vaccine strategies. Moreover, these models are helpful in determining the sufficient number of factors to analyze the experimental results and explain the biological phenomena ^[3]. Stability analysis of the model's equilibria can help researchers to (i) expect the qualitative features of the model for a given set of values of the model's parameters, (ii) establish the conditions that ensure the persistence or deletion of this infection, and (iii) determine under what conditions the immune system is stimulated against the infection.

1.1.1.   Mathematical models of IAV single-infection

Mathematical models of within-host IAV single-infection were developed in several works (see the review papers ^{[25,26,27,28,29]}). Baccam et al. ^[30] presented the following IAV-single-infection with limited target cells and eclipse (or latent) phase:

where $X = X(t)$, $E = E(t)$, $I = I(t)$ and $P = P(t)$ are the concentrations of uninfected epithelial cells, latent IAV-infected epithelial cells, active IAV-infected epithelial cells and free IAV particles, at time $t$, respectively. The model was fitted using real data from six patients infected with influenza ^[30].

Several works were devoted to developing IAV single-infection dynamics models by considering the following:

● Innate immune response: It represents the first line of defense that recognizes the antigens and activate the adaptive immune response. In ^[30], the effect of interferon (IFN) response was included in the IAV infection model. The dynamics of the IFN are given as:

where, $F$ is the concentration of IFN, $\varpi_{F}$ is the IFN production rate constant, $\mu_{F}$ is the IFN decay rate constant, and $\tau$ is the time lag that occurs between the initiation of an IAV infection and the appearance of IFN. IFN can reduce viral replication in an infected cell, the rate of viral production in the presence of IFN was modeled by replacing $\delta_{E}$ by $\frac{\delta_{E} }{1+\epsilon_{E}F}$. The rate that IAV-infected cells in the eclipse phase begin virus production ($\kappa_{P}$) may also be altered and was accounted for by replacing this parameter by $\frac{\kappa_{P} }{1+\epsilon_{P}F}$. The efficiency of these interferon effects is reflected by the parameters $\epsilon_{E}$ and $\epsilon_{P}$ ^[30,31]. Saenz et al. ^[32] presented an influenza model with interferon response. It is reported that the model with interferon response provided better fitting with real data than that without interferon response.

● Adaptive immune response: Cytotoxic T Lymphocytes (CTLs) and antibodies are the two major components of the adaptive immune response. CTLs destroy the viral-infected cells, while the antibodies neutralize the viruses. An influenza dynamics model with different forms of the CTL response was developed in ^[33]. It was shown that slight changes in the virus dynamics was observed when different choices of CTL response were implemented. Both CTL and antibody immunities were included into the IAV model in ^[34].

● Both innate and adaptive immune responses ^{[3,35,36,37,38]}. The model presented in ^[37] predicted that, the level and time of the viral peak are affected by the innate interferon response, while the clearance phase and duration of infection are determined by the CTL response. Handel et al. ^[38] showed that, both the innate and adaptive immune responses are required to give an appropriate explanation of the real data.

● Drug therapy: There are two approved anti-IAV drugs, adamantane antiviral drugs which block infection by reducing the rate of infection, and neuraminidase inhibitors which block the production of newly formed virions ^[31]. Beauchemin et al. ^[39] used model (1.1) to study the effect of the antiviral drug amantadine on IAV infection. Handel et al. ^[40] presented a mathematical model for within-host influenza infection under the effect of neuraminidase inhibitors drugs. Lee et al. ^[34] included the effect of a combination of neuraminidase inhibitors and anti-IAV therapies in the IAV model. The IAV model predicts that the drug therapies are more beneficial when they are administered early.

● Regrowth and death of the uninfected epithelial cells. In ^[34], the first equation of model (1.1) was modified by considering the target cell production and death as:

where $X(0)$ is the initial concentration of the uninfected epithelial cells.

Mathematical analysis of within-host IAV infection model was studied in a few papers ^[33,41,42].

1.1.2.   Mathematical models of SARS-CoV-2 single-infection

Model (1.1) was utilized to characterize the dynamics of SARS-CoV-2 within a host in ^[43,44,45]. Li et al. ^[46], used Eq (1.2) for the uninfected epithelial cell dynamics in case of SARS-CoV-2 infection. The model with target-cell limited and model with regrowth and death of the uninfected epithelial cells presented, respectively, in ^[43,46] were extended and modified by including (i) effect of immune response ^{[44,47,48,49,50,51]}, (ii) effect of different drug therapies ^[52,53], and (iii) effect of time delay ^[54].

Stability analysis of within-host SARS-CoV-2 single-infection models was investigated in ^{[49,50,51,55]}.

Mathematical model of IAV/SARS-CoV-2 coinfection. Recently, mathematical models were developed to characterize within-host co-dynamics of COVID-19 with other diseases, such as: SARS-CoV-2-cancer ^[56], SARS-CoV-2/HIV coinfection ^[57], SARS-CoV-2/malaria coinfection ^[58]. Based on the target cell-limited model (1.1), and the Pinky and Dobrovolny ^[24] developed a model for the within-host dynamics of two respiratory viruses coinfection (SARS-CoV-2 and IAV).

where $L = L(t)$, $Y = Y(t)$ and $V = V(t)$ are the concentrations of latent SARS-CoV-2-infected epithelial cells, active SARS-CoV-2-infected epithelial cells and free SARS-CoV-2 particles, at time $t$, respectively. Model (1.3) describes the competition between two respiratory viruses, SARS-CoV-2 and IAV. However, the effect of the immune response was not modeled. Further, the regrowth and death of the uninfected epithelial cells were not considered. Furthermore, mathematical analysis of the model was not studied.

The objective of the present work is to formulate a mathematical model for within-host IAV/SARS-CoV-2 coinfection with eclipse phase. The model is a generalization of the model (1.3) by taking into account (i) the regrowth and death of the uninfected epithelial cells, (ii) the death of the latent SARS-CoV-2-infected cells and latent IAV-infected cells, (iii) the effect of SARS-CoV-2-specific antibody and IAV-specific antibody. We study the basic qualitative properties of the model, calculate all equilibria, investigate the global stability of equilibria and demonstrate the theoretical results via numerical simulations. We discuss the importance of including the antibody immunity in the IAV/SARS-CoV-2 co-infection model.

Our proposed model can be helpful to characterize the dynamics of coinfection with SARS-CoV-2 strains (Alpha, Beta, Gamma, Delta, Lambda and Omicron), or coinfection of SARS-CoV-2 (or IAV) and other respiratory viruses. Moreover, the model may help to predict new treatment regimens for viral coinfections.

2.   Model formulation

In this section, we present an IAV/SARS-CoV-2 coinfection dynamics model with a latent phase. The dynamics of IAV/SARS-CoV-2 coinfection is presented in the diagram Figure 1. We denote $Z = Z(t)$ and $M = M(t)$ for the concentrations of SARS-CoV-2-specific antibodies and IAV-specific antibodies, at time $t$, respectively. The ODEs that describe the coinfection dynamics are:

where $(X, L, E, Y, I, V, P, Z, M) = (X(t), L(t), E(t), Y(t), I(t), V(t), P(t), Z(t), M(t))$.

In model (2.1) the regrowth death of the uninfected epithelial cells is considered. Further, the death of the latent SARS-CoV-2-infected and latent IAV-infected cells are included, Furthermore, the effect of SARS-CoV-2-specific and IAV-specific antibodies are modeled. First, we start our mathematical analysis of the system by examining the nonnegativity and boundedness of the system's solutions.

3.   Basic qualitative properties

Here, we study the basic qualitative properties of system (2.1).

Lemma 1. The solutions of system (2.1) are nonnegative and bounded.

Proof. We have that

This guarantees that, $(X(t), L(t), E(t), Y(t), I(t), V(t), P(t), Z(t), M(t))\in \mathbb{R}_{\geq0}^{9}$ for all $t\geq0$ when $(X(0), L(0), E(0), Y(0), I(0), V(0), P(0), Z(0), M(0))\in \mathbb{R}_{\geq0}^{9}$. Let us define

Then

where $\phi = \min \{ \alpha, \eta_{L}, \eta_{E}, \frac{\gamma_{Y}}{2}, \frac {\gamma_{I}}{2}, \pi_{V}, \pi_{P}, \mu_{Z}, \mu_{M}\}$. Hence, $0\leq \Psi (t)\leq \Delta_{1}$ if $\Psi(0)\leq \Delta_{1}$ for $t\geq0,$ where $\Delta _{1} = \frac{\lambda}{\phi}.$ Since $X, L, E, Y, I, V, P, Z$ and $M$ are all nonnegative, then $0\leq X(t), L(t), E(t), Y(t), I(t)\leq \Delta_{1}$, $0\leq V(t)\leq \Delta_{2}$, $0\leq P(t)\leq \Delta_{3}$, $0\leq Z(t)\leq \Delta_{4}$, $0\leq M(t)\leq \Delta_{5}$ if $X(0)+L(0)+E(0)+Y(0)+I(0)+\frac{\gamma_{Y} }{2\kappa_{V}}V(0)+\frac{\gamma_{I}}{2\kappa_{P}}P(0)+\frac{\gamma _{Y}\varkappa_{V}}{2\kappa_{V}\sigma_{Z}}Z(0)+\frac{\gamma_{I}\varkappa_{P} }{2\kappa_{P}\sigma_{M}}M(0)\leq \Delta_{1}$, where $\Delta_{2} = \frac {2\kappa_{V}}{\gamma_{Y}}\Delta_{1}$, $\Delta_{3} = \frac{2\kappa_{P}} {\gamma_{I}}\Delta_{1}$, $\Delta_{4} = \frac{2\sigma_{Z}\kappa_{V}} {\varkappa_{V}\gamma_{Y}}\Delta_{1}$ and $\Delta_{5} = \frac{2\sigma_{M} \kappa_{P}}{\varkappa_{P}\gamma_{I}}\Delta_{1}.$ This proves the boundedness of the solutions.

4.   Equilibria

Here, we calculate the system's equilibria and deduce the conditions of their existence. Any equilibrium point $\Xi = \; (X, L, E, Y, I, V, P, Z, M)$ satisfies:

Solving Eqs (4.1)–(4.9), we get eight equilibria.

(i) Infection-free equilibrium, $\Xi_{0} = (X_{0}, 0, 0, 0, 0, 0, 0, 0, 0)$, where $X_{0} = \lambda/\alpha$.

(ii) SARS-CoV-2 single-infection equilibrium without antibody immunity $\Xi_{1} = (X_{1}, L_{1}, 0, Y_{1}, 0, V_{1}, 0, 0, 0),$ where

Therefore, $L_{1} > 0, \; Y_{1} > 0$ and $V_{1} > 0$ when $\frac{X_{0}\kappa_{V} \beta_{V}\delta_{L}}{\gamma_{Y}\pi_{V}\left(\eta_{L}+\delta_{L}\right) } > 1$. We define the basic SARS-CoV-2 single-infection reproductive ratio as:

The parameter $\Re_{1}$ determines whether or not a SARS-CoV-2 single-infection can be established. Thus, we can write

It follows that, $\Xi_{1}\, \ $exists if $\Re_{1} > 1.$

(iii) IAV single-infection equilibrium without antibody immunity, $\Xi_{2} = (X_{2}, 0, E_{2}, 0, I_{2}, 0, P_{2}, 0, 0)$, where

Therefore, $E_{2} > 0, \; I_{2} > 0$ and $P_{2} > 0$ when $\frac{X_{0}\kappa_{P} \beta_{P}\delta_{E}}{\gamma_{I}\pi_{P}\left(\eta_{E}+\delta_{E}\right) } > 1$. We define the basic IAV-infection reproductive ratio as:

The parameter $\Re_{2},$ determines whether or not the IAV single-infection can be established. In terms of $\Re_{2}$, we can write:

Therefore, $\Xi_{2}\, \ $exists if $\Re_{2} > 1$

(iv) SARS-CoV-2 single-infection equilibrium with stimulated SARS-CoV-2-specific antibody immunity, $\Xi_{3} = (X_{3}, L_{3}, 0, Y_{3}, 0, V_{3}, 0, Z_{3}, 0)$, where

We note that $\Xi_{3}$ exists when $\frac{\lambda \beta_{V}\sigma_{Z}\kappa _{V}\delta_{L}}{\gamma_{Y}\pi_{V}\left(\eta_{L}+\delta_{L}\right) \left(\beta_{V}\mu_{Z}+\alpha \sigma_{Z}\right) } > 1$. Let us define the SARS-CoV-2-specific antibody activation ratio in case of SARS-CoV-2 single-infection as:

Thus, $Z_{3} = \frac{\pi_{V}}{\varkappa_{V}}\left(\Re_{3}-1\right)$. The parameter $\Re_{3}$ determines whether or not the SARS-CoV-2-specific antibody immunity is activated in the absence of IAV infection.

(v) IAV single-infection equilibrium with stimulated of IAV-specific antibody immunity, $\Xi_{4} = (X_{4}, 0, E_{4}, 0, I_{4}, 0, P_{4}, 0, M_{4})$, where

We note that $\Xi_{4}$ exists when $\frac{\lambda \beta_{P}\sigma_{M}\kappa _{P}\delta_{E}}{\gamma_{I}\pi_{P}\left(\eta_{E}+\delta_{E}\right) \left(\beta_{P}\mu_{M}+\alpha \sigma_{M}\right) } > 1$. We define the IAV-specific antibody immunity activation ratio for IAV single-infection as:

Thus, $M_{4} = \frac{\pi_{P}}{\varkappa_{P}}\left(\Re_{4}-1\right)$. The parameter $\Re_{4}$ determines whether or not the IAV-specific antibody immunity is activated in the absence of SARS-CoV-2 infection.

(vi) IAV/SARS-CoV-2 coinfection equilibrium with only stimulated SARS-CoV-2-specific antibody immunity, $\Xi_{5} = (X_{5}, L_{5}, E_{5}, Y_{5}, I_{5}, V_{5}, P_{5}, Z_{5}, 0)$, where

We note that $\Xi_{5}$ exists when,

Now, we define the SARS-CoV-2 infection reproductive ratio in the presence of IAV infection as:

The parameter $\Re_{5}$ determines whether or not SARS-CoV-2-infected patients could be coinfected with IAV. Hence,

and then $\Xi_{5}$ exists if $\frac{\Re_{1}}{\Re_{2}} > 1$ and $\Re_{5} > 1.$

(vii) IAV/SARS-CoV-2 coinfection equilibrium with only stimulated IAV-specific antibody immunity, $\Xi_{6} = (X_{6}, L_{6}, E_{6}, Y_{6}, I_{6}, V_{6}, P_{6}, 0, M_{6})$, where

We note that $\Xi_{6}$ exists when

We define the SARS-CoV-2 infection reproductive ratio in the presence of IAV infection as:

Thus,

The parameter $\Re_{6}$ determines whether or not SARS-CoV-2-infected patients could be coinfected with IAV.

(viii) IAV/SARS-CoV-2 coinfection equilibrium with stimulated both SARS-CoV-2-specific and IAV-specific antibody immunities, $\Xi_{7} = (X_{7}, L_{7}, E_{7}, Y_{7}, I_{7}, V_{7}, P_{7}, Z_{7}, M_{7})$, where

It is obvious that $\Xi_{7}$ exists when

Now, we define

Here, $\Re_{7}$ is the SARS-CoV-2-specific antibody activation ratio in case of IAV/SARS-CoV-2 coinfection, $\Re_{8}$ is the IAV-specific antibody activation ratio in case of IAV/SARS-CoV-2 coinfection. Hence, $Z_{7} = \frac{\pi_{V}}{\varkappa_{V}}\left(\Re_{7}-1\right)$ and $M_{7} = \frac {\pi_{P}}{\varkappa_{P}}\left(\Re_{8}-1\right)$. If $\Re_{7} > 1$ and $\Re_{8} > 1$, then $\Xi_{7}$ exists.

In summary, we have eight threshold parameters which determine the existence of the model's equilibria

5.   Global stability

This section is devoted to studying the global asymptotic stability of all equilibria. We configure Lyapunov functions following the way outlined in ^[59]. The following arithmetic-mean-geometric-mean inequality will be utilized:

Let a function $\Lambda_{j}(X, L, E, Y, I, V, P, Z, M)$ and $\tilde{\Omega}_{j}$ be the largest invariant subset of

Define a function

Theorem 1. If $\Re_{1}\leq1$ and $\Re_{2}\leq1$, then $\Xi_{0}$ is globally asymptotically stable (G.A.S).

Proof. Define

We note that, $\Lambda_{0} > 0$ for all $X, L, E, Y, I, V, P, Z, M > 0$ and $\Lambda _{0}(X_{0}, 0, 0, 0, 0, 0, 0, 0, 0) = 0$. We calculate $\frac{d\Lambda_{0}}{dt}$ along the solutions of model (2.1) as:

Using the equilibrium condition, $\lambda = \alpha X_{0}$ we get

Since $\Re_{1}\leq1$ and $\Re_{2}\leq1$, then $\frac{d\Lambda_{0}}{dt}\leq0$ for all $X, Y, I, Z, M > 0$. In addition $\frac{d\Lambda_{0}}{dt} = 0$, when $X = X_{0}$ and $Y = I = Z = M = 0.$ The solutions of system (2.1) tend to $\tilde{\Omega }_{0}$ ^[60] which includes elements with $Y = I = 0$. Thus, $\dot{Y} = \dot{I} = 0$ and from the fourth and fifth equations of system (2.1) we have:

In addition, from the second and third equations of system (2.1) we have:

Therefore, $\tilde{\Omega}_{0} = \left \{ \Xi_{0}\right \}$ and applying Lyapunov-LaSalle Asymptotic Stability Theorem (L-LAST) ^[61,63], we obtain that $\Xi_{0}$ is G.A.S.

Theorem 2. Suppose that $\Re_{1} > 1$, $\Re_{2}/\Re_{1}\leq1$ and $\Re_{3}\leq1$, then $\Xi_{1}$ is G.A.S.

Proof. Let us formulate a Lyapunov function $\Lambda_{1}$ as:

We calculate $\frac{d\Lambda_{1}}{dt}$ as:

Simplifying Eq (5.2), we get

Using the equilibrium conditions for $\Xi_{1}$:

we obtain

Then collecting terms of (5.3), we get

Using inequality (5.1), we get

Since $\Re_{2}/\Re_{1}\leq1$ and $\Re_{3}\leq1,$ then $\frac{d\Lambda_{1}} {dt}\leq0$ for all $X, L, Y, I, V, Z, M > 0$. Moreover, $\frac{d\Lambda_{1}}{dt} = 0$ when $X = X_{1}$, $L = L_{1}, \; Y = Y_{1}$, $V = V_{1}$ and $I = Z = M = 0.$ The solutions of system (2.1) tend to $\tilde{\Omega}_{1}$ where $I = 0$. Hence, $\dot {I} = 0$ and the fifth equation of system (2.1) gives

In addition, from the third equation of system (2.1) we get,

Hence, $\tilde{\Omega}_{1} = \left \{ \Xi \text{ _{1} }\right \}$ and $\Xi_{1}$ is G.A.S by using L-LAST ^[61,62,63].

Theorem 3. Let $\Re_{2} > 1$, $\Re_{1}/\Re_{2}\leq1$ and $\Re_{4}\leq 1$, then $\Xi_{2}$ is G.A.S.

Proof. Consider

We calculate $\frac{d\Lambda_{2}}{dt}$ as:

Then simplifying Eq (5.4), we get

Using the equilibrium conditions for $\Xi_{2}$:

we obtain,

If $\Re_{1}/\Re_{2}\leq1$ and $\Re_{4}\leq1$, then employing inequality (5.1), we get $\frac{d\Lambda_{2}}{dt}\leq0$ for all $X, E, Y, I, P, Z, M > 0$. Further, $\frac{d\Lambda_{2}}{dt} = 0$ when $X = X_{2}, \; E = E_{2}$, $I = I_{2},$ $P = P_{2}$ and $Y = Z = M = 0.$ The solutions of system (2.1) tend to $\tilde{\Omega}_{2}$ which has $Y = 0$, and gives $\dot{Y} = 0$. The fourth equation of system (2.1) gives

In addition, from the second equation of system (2.1) gives

Therefore, $\tilde{\Omega}_{2} = \left \{ \Xi \text{ _{2} }\right \}$. Applying L-LAST, we get $\Xi_{2}$ is G.A.S.

Theorem 4. Let $\Re_{3} > 1$ and $\Re_{5}\leq1$, then $\Xi_{3}$ is G.A.S.

Proof. Define

We calculate $\frac{d\Lambda_{3}}{dt}$ as:

Then simplifying Eq (5.5), we get

Using the equilibrium conditions for $\Xi_{3}$:

we obtain,

Using inequality (5.1) and $\Re_{5}\leq1,$ we get $\frac{d\Lambda_{3} }{dt}\leq0$ for all $X, L, Y, V, P, M > 0$. Further, $\frac{d\Lambda_{3}}{dt} = 0$ when $X = X_{3}$, $L = L_{3}$, $Y = Y_{3}, \; V = V_{3}$ and $P = M = 0.$ Further, the trajectories of system (2.1) tend to $\tilde{\Omega}_{3}$ which has elements with $V = V_{3}$ and $P = 0$. Then $\dot{V} = 0$ and $\dot{P} = 0.$ The sixth and seventh equations of system (2.1), provide

In addition, from the fifth equation of system (2.1) gives

Consequently, $\tilde{\Omega}_{3} = \left \{ \Xi \text{ _{3} }\right \}$. Applying L-LAST, we find that $\Xi_{3}$ is G.A.S.

Theorem 5. If $\Re_{4} > 1$ and $\Re_{6}\leq1$, then $\Xi_{4}$ is G.A.S.

Proof. Define a function $\Lambda_{4}$ as:

Calculating $\frac{d\Lambda_{4}}{dt}$ as:

Equation (5.6) can be written as:

Using the equilibrium conditions for $\Xi_{4}$:

we obtain,

Since $\Re_{6}\leq1$, then employing inequality (5.1), we get $\frac{d\Lambda_{4}}{dt}\leq0$ for all $X, E, I, V, P, Z > 0$. Further, $\frac{d\Lambda_{4}}{dt} = 0$ when $X = X_{4}$, $E = E_{4}, \; I = I_{4}, \; P = P_{4}$ and $V = Z = 0.$ The solutions of system (2.1) tend to $\tilde{\Omega}_{4}$ which contains elements with $P = P_{4}$ and $V = 0$, then $\dot{V} = \dot{P} = 0$. The sixth and seventh equations of system (2.1) imply

In addition, since $Y = 0$, then $\dot{Y} = 0$. The fourth equation of system (2.1) gives

Therefore, $\tilde{\Omega}_{4} = \left \{ \Xi_{4}\right \}.$ Applying L-LAST, we get $\Xi_{4}$ is G.A.S.

Theorem 6. If $\Re_{5} > 1$, $\Re_{1}/\Re_{2} > 1$ and $\Re_{8}\leq1$, then $\Xi_{5}$ is G.A.S.

Proof. Define

Calculating $\frac{d\Lambda_{5}}{dt}$ as:

Equation (5.7) can be simplified as:

Using the equilibrium conditions for $\Xi_{5}$:

we obtain,

Since $\Re_{8}\leq1$, then employing inequality (5.1), we get $\frac{d\Lambda_{5}}{dt}\leq0$ for all $X, L, E, Y, I, V, P, M > 0$. Moreover, we have $\frac{d\Lambda_{5}}{dt} = 0,$ when $X = X_{5}, \; L = L_{5}, \; E = E_{5}$, $Y = Y_{5},$ $I = I_{5}, \; V = V_{5}, \; P = P_{5}$ and $M = 0.$ The trajectories of system (2.1) converge to $\tilde{\Omega}_{5}$ which comprises elements with $Y = Y_{5}$ and $V = V_{5},$ then $\dot{V} = 0$. The sixth equation of system (2.1) implies that

Consequently, $\tilde{\Omega}_{5} = \left \{ \Xi \text{ _{5} }\right \}.$ and by applying L-LAST, we get $\Xi_{5}$ is G.A.S.

Theorem 7. Let $\Re_{6} > 1$, $\Re_{7}\leq1$ and $\Re_{2}/\Re_{1} > 1$, then $\Xi_{6}$ is G.A.S.

Proof. Consider a function $\Lambda_{6}$ as:

Calculating $\frac{d\Lambda_{6}}{dt}$ as:

We collect the terms of Eq (5.8) as:

Using the equilibrium conditions for $\Xi_{6}$:

we obtain,

Since $\Re_{7}\leq1$, then employing inequality (5.1), we get $\frac{d\Lambda_{6}}{dt}\leq0$ for all $X, L, E, Y, I, V, P, Z > 0$. Moreover, $\frac{d\Lambda_{6}}{dt} = 0$ when $X = X_{6}$, $L = L_{6}, \; E = E_{6}, \; Y = Y_{6},$ $I = I_{6}, \; V = V_{6}, \; P = P_{6}$ and $Z = 0.$ The solutions of system (2.1) tend to $\tilde{\Omega}_{6}$ which contains elements with $P = P_{6}$ then, $\dot{P} = 0$. The seven equation of system (2.1) implies that

Consequently, $\tilde{\Omega}_{6} = \left \{ \Xi_{6}\right \}$. Using L-LAST we deduce that $\Xi_{6}$ is G.A.S.

Theorem 8. If $\Re_{7} > 1$ and $\Re_{8} > 1$, then $\Xi_{7}$ is G.A.S.

Proof. Define a function $\Lambda_{7}$ as:

Calculating $\frac{d\Lambda_{7}}{dt}$ as:

We collect the terms of Eq (5.9) as:

Using the equilibrium conditions for $\Xi_{7}$:

we obtain,

Using inequality (5.1), we get $\frac{d\Lambda_{7}}{dt}\leq0$ for all $X, L, E, Y, I, V, P > 0$, where $\frac{d\Lambda_{7}}{dt} = 0$ when $X = X_{7}$, $L = L_{7}, \; E = E_{7}, \; Y = Y_{7}, \; I = I_{7}, \; V = V_{7}$ and $P = P_{7}.$ The solutions of system (2.1) tend to $\tilde{\Omega}_{7}$ which includes element with $V = V_{7}$ and $P = P_{7}$ which gives $\dot{V} = \dot{P} = 0$ and from the sixth and seventh equations of system (2.1) we get

Therefore, $\tilde{\Omega}_{7} = \left \{ \Xi_{7}\right \}$ and by employing L-LAST, we get $\Xi_{7}$ is G.A.S. Based on the above findings, we summarize the existence and global stability conditions for all equilibrium points in Table 1.

6.   Effect of the antibody immunity on the IAV/SARS-CoV-2 coinfection dynamics

We noted that system (2.1) has eight equilibria for which the coexistence case of IAV and SARS-CoV-2 can only be occurred if at least one type of the specific antibody immunities is active. Now, we discuss the importance of considering the antibody immune response in the IAV/SARS-CoV-2 dynamics model. If the antibody immune response is neglected then system (2.1) becomes:

We can see that system (6.1) describes the competition between IAV and SARS-CoV-2 on one source of target cells, epithelial cells. The model admits only three equilibria:

(i) Infection-free equilibrium, $\tilde{\Xi}_{0} = (\tilde{X} _{0}, 0, 0, 0, 0, 0, 0)$, where both IAV and SARS-CoV-2 are cleared,

(ii) SARS-CoV-2 single-infection equilibrium $\tilde{\Xi}_{1} = (\tilde{X}_{1}, \tilde{L}_{1}, 0, \tilde{Y}_{1}, 0, \tilde{V}_{1}, 0)$, where the IAV is blocked,

(iii) IAV single-infection equilibrium, $\tilde{\Xi}_{2} = (\tilde {X}_{2}, 0, \tilde{E}_{2}, 0, \tilde{I}_{2}, 0, \tilde{P}_{2})$, where the SARS-CoV-2 is blocked, where $\tilde{X}_{i} = X_{i}$, $i = 0, 1, 2$, $\tilde{L}_{1} = L_{1}$, $\tilde{Y}_{1} = Y_{1}$, $\tilde{V}_{1} = V_{1}$, $\tilde{E}_{2} = E_{2}$, $\tilde {I}_{2} = I_{2},$ and $\tilde{P}_{2} = P_{2}$.

We note that the case of IAV and SARS-CoV-2 coexistence does not appear. In the recent studies presented in ^[5,9,11,12], it was recorded that some COVID-19 patients were detected to be coinfected with IAV. Therefore, neglecting the immune response may not describe the coinfection dynamics accurately. This supports the idea of including the immune response into the IAV/SARS-CoV-2 coinfection model, where the case of IAV and SARS-CoV-2 coexistence is observed.

7.   Numerical simulations

The global stability of the system's equilibria will be illustrated numerically. We use the values of the parameters presented in Table 2. In addition, we make a comparison between single-infection and coinfection.

7.1.   Stability of the equilibria

Now, we solve system (2.1) with three different initial conditions (states) as:

Selecting the values of $\beta_{V}$, $\beta_{P}$, $\sigma_{Z}$ and $\sigma_{M}$ leads to the following situations:

Situation 1 (Stability of $\Xi_{0}$): $\beta_{V} = 0.001, \; \beta_{P} = 0.001$, $\sigma_{Z} = 0.01$ and $\sigma_{M} = 0.02$. For these values of parameters, we have $\Re_{1} = 0.0455 < 1$ and $\Re_{2} = 0.1 < 1$. Figure 2 shows that the trajectories tend to the equilibrium $\Xi_{0} = (10, 0, 0, 0, 0, 0, 0, 0, 0)$ for all initials C1-C3. This demonstrates that, $\Xi_{0}$ is G.A.S based on Theorem 1. In this situation, both SARS-CoV-2 and IAV will be removed.

Situation 2 (Stability of $\Xi_{1}$): $\beta_{V} = 0.05, \; \beta_{P} = 0.001,$$\sigma_{Z} = 0.002$ and $\sigma_{M} = 0.02$. With such selection we obtain $\Re_{1} = 2.2727 > 1$, $\Re _{3} = 0.0874 < 1$ and hence $\Re_{2}/\Re_{1} = 0.044 < 1$. The equilibrium point $\Xi_{1}$ exists with $\Xi_{1} = (4.4, 2.8, 0, 1.27, 0, 1.27, 0, 0, 0)$. It is clear from Figure 3 that, the trajectories tend to $\Xi_{1}$ for all initials C1-C3. Thus, the numerical results agree with Theorem 2. This case simulates a SARS-CoV-2 single-infection without antibody immunity.

Situation 3 (Stability of $\Xi_{2}$): $\beta_{V} = 0.005, \; \beta_{P} = 0.03$, $\sigma_{Z} = 0.01$ and $\sigma_{M} = 0.001$. This gives $\Re_{2} = 3 > 1$, $\Re_{4} = 0.12 < 1$ and then $\Re_{1}/\Re_{2} = 0.0758 < 1$. The numerical results show that, $\Xi_{2} = \left(3.33, 0, 2.78, 0, 0.83, 0, 3.33, 0, 0\right)$ exists. We can observe from Figure 4 that, the trajectories converge to $\Xi_{2}$ regardless of the initial states. This result supports the result of Theorem 3. This situation represents an IAV single-infection without antibody immunity.

Situation 4 (Stability of $\Xi_{3}$): $\beta_{V} = 0.09, \; \beta_{P} = 0.002$, $\sigma_{Z} = 0.05$ and $\sigma_{M} = 0.05$. This yields $\Re_{3} = 1.461 > 1$ and $\Re_{5} = 0.0714 < 1$. Figure 5 shows that the trajectories tend to $\Xi_{3} = (3.57, 3.21, 0, 1.46, 0, 1, 0, 1.84, 0)$ regardless of the initial states. Therefore, $\Xi_{3}$ is G.A.S and this supports Theorem 4. Hence, a SARS-CoV-2 single-infection with stimulated SARS-CoV-2-specific antibody is attained.

Situation 5 (Stability of $\Xi_{4}$): $\beta_{V} = 0.01, \; \beta_{P} = 0.1$, $\sigma_{Z} = 0.01$ and $\sigma _{M} = 0.02$. The values of $\Re_{4}$ and $\Re_{6}$ are computed as $\Re _{4} = 2 > 1$ and $\Re_{6} = 0.0909 < 1$. Thus, $\Xi_{4}$ exists with $\Xi _{4} = (2, 0, 3.33, 0, 1, 0, 2, 0, 2.5)$. In Figure 6 we see that the trajectories tend to $\Xi_{4}$ regardless of the initial states. It follows that $\Xi_{4}$ is G.A.S according to Theorem 5. Hence, an IAV single-infection with activated IAV-specific antibody is achieved.

Situation 6 (Stability of $\Xi_{5}$): $\beta_{V} = 0.09, \; \beta_{P} = 0.02$, $\sigma_{Z} = 0.095$ and $\sigma_{M} = 0.009$. Then, we calculate $\Re_{5} = 1.027 > 1$, $\Re_{8} = 0.5369 < 1$ and $\Re_{1}/\Re_{2} = 2.0455 > 1$ The numerical results drawn in Figure 7 show that $\Xi_{5} = (5, 2.37, 0.11, 1.08, 0.03, 0.53, 0.13, 4.18, 0)$ exists and it is G.A.S and this is consistent with Theorem 6. As a result, a coinfection with SARS-CoV-2 and IAV is attained where only SARS-CoV-2-specific antibody is stimulated.

Situation 7 (Stability of $\Xi_{6}$): $\beta_{V} = 0.09, \; \beta_{P} = 0.09$, $\sigma_{Z} = 0.03$ and $\sigma _{M} = 0.03$. We compute $\Re_{6} = 1.2032 > 1$, $\Re_{7} = 0.6392 < 1$ and $\Re_{2} /\Re_{1} = 2.2 > 1$. We find that, the equilibrium $\Xi_{6} = \left(2.44, 0.84, 2.44, 0.38, 0.73, 0.38, 1.33, 0, 3\right)$ exists. Further, the numerical solutions outlined in Figure 8 show that, $\Xi_{6}$ is G.A.S and this boosts the result of Theorem 7. In this situation, a coinfection with SARS-CoV-2 and IAV are attained where only IAV-specific antibody is activated.

Situation 8 (Stability of $\Xi_{7}$): $\beta_{V} = 0.09, \; \beta_{P} = 0.09$, $\sigma_{Z} = 0.5$ and $\sigma _{M} = 0.5$. This selection yields $\Re_{7} = 3.0898 > 1$ and $\Re_{8} = 6.7976 > 1$. Figure 9 shows that $\Xi_{7} = (7.55, 0.68, 0.45, 0.31, 0.14, 0.1, 0.08, 8.36, 14.49)$ exists and it is G.A.S based on Theorem 8. In this situation, a coinfection with SARS-CoV-2 and IAV is established regardless of the initial states. In this case, both SARS-CoV-2-specific antibodies and IAV-specific antibodies are working against the coinfection.

For more confirmation, we investigate the local stability of the system's equilibria. Calculating the Jacobian matrix $J = J(X, L, E, Y, I, V, P, Z, M)$ of system (2.1) as:

where $J_{11} = -\left(\alpha+\beta_{V}V+\beta_{P}P\right), \ \ J_{22} = -\left(\eta_{L}+\delta_{L}\right), \ \ J_{33} = -\left(\eta_{E}+\delta_{E}\right), \ \ J_{66} = -\left(\pi_{V}+\varkappa_{V}Z\right)$, $J_{77} = -\left(\pi_{P}+\varkappa_{P}M\right), \ \ J_{88} = \sigma_{Z}V-\mu_{Z}, \ \ J_{99} = \sigma_{M}P-\mu_{M}.$

At each equilibrium, we compute the eigenvalues $\lambda_{j}, \; j = 1, 2, ..., 9$ of $J.$ If $\operatorname{Re}(\lambda_{j}) < 0, \; i = 1, 2, ..., 9$, then the equilibrium point is locally stable. We select the parameters $\beta_{V}$, $\beta_{P}$, $\sigma_{Z}$ and $\sigma_{M}$ as given in situations 1-8, then we compute all nonnegative equilibria and the accompanying eigenvalues. Table 3 outlined the nonnegative equilibria, the real parts of the eigenvalues and whether or not the equilibrium point is stable.

7.2.   Comparison results

In this subsection, we present a comparison between the single-infection and coinfection.

7.2.1.   Influence of IAV infection on the dynamics of SARS-CoV-2 single-infection

Here, we compare the solutions of model (2.1) and the following SARS-CoV-2 single-infection model:

We fix parameters $\beta_{V} = 0.09$, $\beta_{P} = 0.05$, $\sigma_{Z} = 0.5$, and $\sigma_{M} = 0.9$ and select the initial state as:

From Figure 10 we observe that when the SARS-CoV-2 single-infected individual is coinfected with IAV, then the concentrations of uninfected epithelial cells, latent SARS-CoV-2-infected cells, active SARS-CoV-2-infected cells and SARS-CoV-2-specific antibodies are reduced. However, the concentration of free SARS-CoV-2 particles tends to be the same value in both SARS-CoV-2 single-infection and IAV/SARS-CoV-2 coinfection. This result agrees with the observation of Ding et al. ^[11] which said that "IAV/SARS-CoV-2 coinfection did not result in worse clinical outcomes in comparison with SARS-CoV-2 single-infection".

7.2.2.   Influence of SARS-CoV-2 infection on the dynamics IAV single-infection

To examine the impact of SARS-CoV-2 infection on IAV single-infection, we compare the solutions of model (2.1) and the following IAV single-infection model:

We fix parameters $\beta_{V} = \; 0.09$, $\beta_{P} = 0.08, \; \sigma_{Z} = 0.07$ and $\sigma_{M} = 0.05$ and consider the following initial condition:

It can be observed from Figure 11 that, when the IAV single-infected individual is coinfected with SARS-CoV-2 then the concentrations of uninfected epithelial cells, latent IAV-infected cells, active IAV-infected cells and IAV-specific antibodies are decreased. However, the concentration of free IAV particles tend to the same value in both IAV single-infection and IAV/SARS-CoV-2 coinfection.

8.   Discussion

IAV and SARS-CoV-2 coinfection cases were reported in some works (see ^[5,9,11,12]). Therefore, it is important to understand the within-host dynamics of this coinfection. In this paper, we develop and examine a within-host IAV/SARS-CoV-2 coinfection model. The model considered the interactions between uninfected epithelial cells, latent SARS-CoV-2-infected cells, latent IAV-infected cells, active SARS-CoV-2-infected cells, active IAV-infected cells, free SARS-CoV-2 particles, free IAV particles, SARS-CoV-2-specific antibodies and IAV-specific antibodies. We examined the nonnegativity and boundedness of the solutions. We found that the system has eight equilibria and we proved the following:

(I) The infection-free equilibrium $\Xi_{0}$ always exists. It is G.A.S when $\Re_{1}\leq1$ and $\Re_{2}\leq1$. In this case, the patient is recovered from both IAV and SARS-CoV-2.

(II) The SARS-CoV-2 single-infection equilibrium without antibody immunity $\Xi_{1}$ exists if $\Re_{1} > 1$. It is G.A.S when $\Re_{1} > 1$, $\Re_{2} /\Re_{1}\leq1$ and $\Re_{3}\leq1$. This case leads to the situation of the patient only infected by SARS-CoV-2 with an inactive immune response. As we will see below that if both SARS-CoV-2-specific antibody and IAV-specific antibody immunities are not activated against the two viruses, then according to the competition between the two viruses, SARS-CoV-2 may be able to block the IAV.

(III)- The IAV single-infection equilibrium without antibody immunity $\Xi _{2}$ exists if $\Re_{2} > 1$. It is G.A.S when $\Re_{2} > 1$, $\Re_{1}/\Re _{2}\leq1$ and $\Re_{4}\leq1$. This case leads to the situation of the patient only infected by IAV with an unstimulated immune response. Then, IAV may be able to block the SARS-CoV-2.

(IV) The SARS-CoV-2 single-infection equilibrium with stimulated SARS-CoV-2-specific antibody immunity $\Xi_{3}$ exists if $\Re_{3} > 1$. It is G.A.S when $\Re_{3} > 1$ and $\Re_{5}\leq1$. This point represents the situation of SARS-CoV-2 single-infection patient with active SARS-CoV-2-specific antibody immunity.

(V) The IAV single-infection equilibrium with stimulated IAV-specific antibody immunity $\Xi_{4}$ exists if $\Re_{4} > 1$. It is G.A.S when $\Re_{4} > 1$ and $\Re_{6}\leq1$. This point represents the case of IAV single-infection patient with active IAV-specific antibody immunity.

(VI) The IAV/SARS-CoV-2 coinfection equilibrium with only stimulated SARS-CoV-2-specific antibody immunity $\Xi_{5}$ exists if $\Re_{5} > 1$ and $\Re_{1}/\Re_{2} > 1$. It is G.A.S when $\Re_{5} > 1$, $\Re_{8}\leq1$ and $\Re _{1}/\Re_{2} > 1$. Here, the IAV/SARS-CoV-2 coinfection occurs with only stimulated SARS-CoV-2-specific antibody immunity.

(VII) The IAV/SARS-CoV-2 coinfection equilibrium with only stimulated IAV-specific antibody immunity $\Xi_{6}$ exists if $\Re_{6} > 1$ and $\Re _{2}/\Re_{1} > 1$. It is G.A.S when $\Re_{6} > 1$, $\Re_{7}\leq1$ and $\Re_{2} /\Re_{1} > 1$. It means that the IAV/SARS-CoV-2 coinfection occurs with only stimulated IAV-specific antibody immunity.

(VIII) The IAV/SARS-CoV-2 coinfection equilibrium with stimulated both SARS-CoV-2-specific antibodies and IAV-specific antibody immunities $\Xi_{7}$ exists and it is G.A.S if $\Re_{7} > 1$ and $\Re_{8} > 1$. It means that, the IAV/SARS-CoV-2 coinfection occurs with both SARS-CoV-2-specific antibodies and IAV-specific antibody immunities are activated.

The global stability of equilibria was established using the Lyapunov method. We performed numerical simulations and demonstrated that they are in good agreement with the theoretical results. We discussed the influence of IAV infection on SARS-CoV-2 single-infection dynamics and vice versa. We found that the concentration of free IAV or SARS-CoV-2 particles cells tends to be the same value in both single-infection and coinfection. This agrees with the work of Ding et al. ^[11] which reported that IAV/SARS-CoV-2 coinfection did not result in worse clinical outcomes. In addition, the spread of seasonal influenza can increase the likelihood of coinfection in patients with COVID-19 ^[9].

The model developed in this work can be improved by (i) utilizing real data to find a good estimation of the parameters' values, (ii) studying the effect of time delays that occur during infection or production of IAV and SARS-CoV-2 particles, (iii) considering viral mutations ^[64,65], (iv) considering the effect of treatments on the progression of both viruses, and (v) including the influence of CTLs in killing SARS-CoV-2-infected and IAV-infected cells. Memory is an important characteristic of viral infections and immune response. It will be important to address the effect of memory on the dynamics of IAV/SARS-CoV-2 coinfection by formulation of the model via fractional differential equations ^[66,67,68].

The innate immune response is one of the major antiviral responses to explain host-pathogen interaction. Also, it is a trigger to induce adaptive immunity which is the major focus of our proposed model. Model (2.1) can be extended to include the effect of IFN response as:

where $\epsilon_{L}$, $\epsilon_{E}$, $\epsilon_{V}$ and $\epsilon_{P}$ are the efficiencies of the IFN effects. These research points need further investigations so we leave them to future works.

Acknowledgment

This research work was funded by Institutional Fund Projects under grant no. (IFPIP:69-130-1443). The authors gratefully acknowledge technical and financial support provided by the Ministry of Education and King Abdulaziz University, DSR, Jeddah, Saudi Arabia.