SEIR Mathematical Model for Influenza-Corona Co-Infection with Treatment and Hospitalization Compartments and Optimal Control Strategies

1  Introduction

Coronavirus or COVID-19 first reported in Wuhan, China, in late 2019 [1,2]. The World Health Organization officially declared the COVID-19 pandemic in March 2020 [3,4], following a series of devastating coronaviruses, such as SARS-CoV and MERS-CoV, that had already spread globally [5,6]. The coronavirus can lead to a deadly cytokine storm, similar to previous coronaviruses [7–9]. The available evidence indicates that the primary mode of transmission of the coronavirus among individuals is through the release of respiratory droplets during activities such as sneezing and coughing [10]. The incubation period for the coronavirus disease typically ranges from 2 to 14 days, with an average of 5 to 6 days [11]. Influenza, commonly called the flu, is an infectious respiratory disease caused by the influenza virus [12,13]. During seasonal pandemics and epidemics, influenza spreads widely throughout the world. Apart from pandemic influenza, there is a global influenza outbreak every year that causes between 3 and 5 million cases of severe disease and between 250,000 and 500,000 fatalities [14]. The time it takes for the signs and symptoms of influenza to appear is about 2 days but might range from 1 to 4 days [15].

Concerns about the possibility of a twin epidemic of influenza and corona have persisted since the start of the pandemic. Many cases have been reported where individuals had both SARS-CoV-2 and influenza [16]. Up to 20% of corona patients, according to some research [17], can have co-infections with other respiratory viruses. The influenza season and the corona pandemic pose a significant threat to public health. Both viruses share similar transmission characteristics and clinical symptoms, causing respiratory infections [17]. The interplay between influenza and corona viruses is a major concern, with fatalities among people infected with both viruses being twice as high as those infected with the new coronavirus [16,17].

Co-infection between corona and other respiratory pathogens is more common in the USA and China, with 7% of SARS-CoV-2-positive patients sharing the burden of co-infection [17]. The clinical signs of corona and influenza include coughing, runny noses, sore throats, fevers, headaches, and fatigue [16,17]. Those who are infected with either virus can experience varied levels of illness, including some who show no symptoms, only minor symptoms, or severe disease [17]. Both the flu and corona can have severe consequences, including pneumonia, respiratory failure, acute respiratory distress syndrome, sepsis, heart attacks or strokes, multiple organ failure, severe inflammation, and even death [18,19].

Both coronavirus and influenza are single-stranded encapsulated RNA viruses; the former contain a positive sense RNA strand and the latter a negative one [20]. Both have comparable infection sites in the upper respiratory tracts (URT) and lower respiratory tracts (LRT). While influenza URT infections are highly transmissible but have a low pathogenicity, influenza LRT infections can produce more severe symptoms [20]. Furthermore, both viruses can spread directly from person to person or indirectly through close contact. Their routes of transmission are comparable and include droplet, aerosol, and self-inoculation by hand contamination [16,17,20]. Co-infection can affect the disease prognosis, increase therapeutic intolerance, and severely compromise the immune system of the host. One of the most important area of epidemiology has been the study of the co-existence and co-infection of multiple diseases [21,22]. To assess the dynamics of co-infection and estimate the treatment facilities, epidemiological models are very helpful. Numerous research works focus on the coexistence of two infectious agents in vulnerable hosts [22–26].

There is a considerable body of work in the literature focusing on coronavirus, influenza, and their co-infection. Epidemiologists are consistently striving to provide comprehensive insights into these epidemics and develop effective control strategies with minimal costs. In [9], the authors formulated an extended susceptible-exposed-symptomatic-asymptomatic-superspreader-quarantine-recovered (SEIAPQR) mathematical model, incorporating two non-pharmaceutical optimal control strategies: quarantine and self-precautions. In [21], the authors constructed a model for co-infection of coronavirus and influenza, introducing vaccination compartments for COVID-19, Influenza, and corona-influenza co-infection. They analyzed the impact of vaccinations, including booster doses for COVID-19, and offered a detailed analysis of influenza-corona co-infection, exploring the effects of transmission, cure, or treatment dynamics. Similarly, in [22], the authors proposed a mathematical model for corona-influenza co-infection, dividing the entire model into three susceptible-exposed-infected-recovered (SEIR) submodels with a vaccination compartment. The dynamics of each submodel are analyzed individually to understand the complete model. Finally, an optimal control strategy is defined in [22] with three optimal controls: non-pharmaceutical strategy (minimization of infectious interactions), vaccine for corona, and vaccination for influenza, respectively.

To investigate the co-dynamics of influenza and corona co-infection, we divide the total population into eight distinct time-dependent classes or compartments: susceptible S(t), exposed E(t), infectious with influenza II(t), infectious with corona IC(t), infectious with co-infection IIC(t), under-treatment T(t), hospitalized H(t), and recovered R(t). First we prove fundamental properties (i.e., positive, bounded, and well-posedness) of the proposed model. Then, we determine the important equilibrium points (disease-free equilibrium (DFE) and endemic equilibrium (EE)) of the proposed model, the reproduction number ℛ0, and demonstrate that DFE and EE points are locally and globally asymptotically stable. To determine the most important critical parameters for ℛ0, we perform sensitivity analysis and mark the critical parameters for disease control analysis. In the second phase of this study, we will conduct an optimal control analysis aimed at minimizing disease spread [25,26]. To achieve this, we will develop two optimal control problems, each using a distinct control strategy. The first strategy focuses on the effect of treatment for infected individuals, aiming to find the optimal treatment rate that reduces infection within the population. In the second strategy, we adjust the model to incorporate non-pharmaceutical interventions (i.e., self-precautionary measures and government-provided resources) along with treatment to assess their impact on preventing disease transmission.

The format of this manuscript is as follows: We explore the development of the influenza-corona co-infection model (Section 2). The essential characteristics of the solution are discussed, including the presence of uniqueness and positive bounded solutions (Section 3). We analyze the equilibrium points and the reproduction number (ℛ0) (Section 4). Specifically, we explore the asymptotic stability of the disease-free equilibrium (DFE) and the endemic equilibrium (EE) points at the local and global levels. A sensitivity analysis of the reproduction number with respect to different parameters is included (Section 5). The study discusses several controls and presents their effects at various control levels and optimal control strategies including numerical simulation (Section 6). Numerical simulations are also included in this section to further highlight the findings. Finally, we provide a concise summary of the manuscript’s key findings (Section 7).

2  Mathematical Model

We divide the total population, denoted as N(t), into eight mutually distinct and time-dependent sub-classes or compartments, which we refer to as follows: susceptible S(t), exposed E(t), influenza-infectious II(t), corona-infectious IC(t), influenza-corona infectious IIC(t), under treatment T(t), hospitalized H(t), and those who have recovered or have been removed from the population R(t). The susceptible population S(t) comprises individuals who are at risk of contracting influenza and corona infections following interactions with infectious individuals carrying either influenza, corona, or both viruses. After such interactions, a portion of the susceptible population transitions to the exposed state E(t), making them vulnerable to coronavirus, influenza, or co-infection with both viruses. Following an incubation period of 2 to 5 days for influenza and 2 to 14 days for corona, those in the exposed category progress to become influenza-infectious II(t), corona-infectious IC(t), or influenza-corona co-infectious IIC(t). Individuals who are currently infected and receiving medical treatment in any form are placed in the treatment compartment T(t). Patients with severe conditions, such as those requiring oxygen beds, are considered for hospitalization H(t). Eventually, individuals who have successfully recovered from the diseases are categorized to move in the R(t) compartment. The total population N(t) can be written as:

N(t)=S(t)+E(t)+II(t)+IC(t)+IIC(t)+T(t)+H(t)+R(t),(1)

and flow of co-infection through compartments is shown in Fig. 1. The transmission and translation rates (α′s, β′s, and γ′s, etc.) along with their complete description are given below:

•   Π: Birth/Recruitment rate of susceptible population.

•   α1: Transmission rate of influenza infection.

•   α2: Transmission rate of corona infection.

•   α3: Transmission rate of influenza-corona co-infection.

•   β1: Translation rate of influenza exposed to influenza infectious.

•   β2: Translation rate of corona exposed to corona infectious.

•   β3: Translation rate of influenza-corona exposed to influenza-corona infectious.

•   γ1: Translation rate of influenza infectious to influenza-corona infectious.

•   γ2: Treatment rate of influenza infectious.

•   γ3: Recovery rate of influenza infectious.

•   τ1: Translation rate of corona infectious to influenza-corona infectious.

•   τ2: Treatment rate of corona infectious.

•   τ3: Recovery rate of corona infectious.

•   δ1: Treatment rate of influenza-corona infectious.

•   δ2: Recovery rate of influenza-corona infectious.

•   ϕ1: Hospitalization rate of under-treatment patients.

•   ϕ2: Recovery rate of under-treatment patients.

•   ϕ3: Recovery rate of hospitalized patients.

•   μ: Natural death rate.

•   d1: Death rate due to influenza disease.

•   d2: Death rate due to corona disease.

•   d3: Death rate due to influenza-corona disease.

•   d4: Death rate due to disease in under-treatment patients.

•   d5: Death rate due to disease in hospitalized patients.

Figure 1: Flow diagram of the proposed influenza-corona co-infection model (Eq. (2))

We model the transmission of influenza-corona co-infection by the following system of non-linear differential equations.

dSdt=Π−(α1II+α2IC+α3IIC)S−μS,(2a)

dEdt=(α1II+α2IC+α3IIC)S−(β1+β2+β3+μ)E,(2b)

dIIdt=β1E−(γ1+γ2+γ3+μ+d1)II,(2c)

dICdt=β2E−(τ1+τ2+τc+μ+d2)IC,(2d)

dIICdt=β3E+γ1II+τ1IC−(δ1+δ2+μ+d3)IIC,(2e)

dTdt=γ2II+τ2IC+δ1IIC−(ϕ1+ϕ2+μ+d4)T,(2f)

dHdt=ϕ1T−(ϕ3+μ+d5)H,(2g)

dRdt=γ3II+τcIC+δ2IIC+ϕ2T+ϕ3H−μR,(2h)

with the non-negative initial conditions:

S(0)>0, E(0)≥0, II(0)≥0, IC(0)≥0, IIC(0)≥0, T(0)≥0, H(0)≥0, R(0)≥0.(2i)

The above autonomous system can be written in compact form as:

dydt=F(y(t)),y(0)=y0,0<t<T<+∞,(3)

where y:[0,T]→R+8 and F:R+8→R+8 are vector valued functions such that

y(t)=(S(t)E(t)II(t)IC(t)IIC(t)T(t)H(t)R(t)),y(0)=(S(0)E(0)II(0)IC(0)IIC(0)T(0)H(0)R(0)),

and

F(y)=(F1F2F3F4F5F6F7F8)=(Π−(α1II+α2IC+α3IIC)S−μS(α1II+α2IC+α3IIC)S−(β1+β2+β3+μ)Eβ1E−(γ1+γ2+γ3+μ+d1)IIβ2E−(τ1+τ2+τc+μ+d2)ICβ3E+γ1II+τ1IC−(δ1+δ2+μ+d3)IICγ2II+τ2IC+δ1IIC−(ϕ1+ϕ2+μ+d4)Tϕ1T−(ϕ3+μ+d5)Hγ3II+τcIC+δ2IIC+ϕ2T+ϕ3H−μR),

respectively.

3  Theoretical Properties and Proofs

The fundamental characteristics of the proposed mathematical model, such as the existence of a unique solution and a positive, bounded solution in a feasible region is demonstrated in this section.

Theorem 3.1. The co-infection model (Eq. (2)) has a bounded solution, y(t)=(S(t),E(t),II(t),IC(t),IIC(t),T(t),H(t),R(t)).

Proof. By differentiating Eq. (1) w.r.t time t and then adding right hand side of ODEs of the system (Eq. (2)), we get an equation of the form

dNdt=Π−μN−(d1II+d2IC+d3IIC+d4T+d5H),(4)

with

N0=N(0)=S(0)+E(0)+II(0)+IC(0)+IIC(0)+T(0)+H(0)+R(0)≤Πμ.(5)

Since, d1II+d2IC+d3IIC+d4T+d5H≥0, hence Eq. (4) can be written as:

dNdt≤Π−μN.(6)

By using the properties of the Laplace transformation and partial fraction on inequality (6), we get

L{dN(t)dt}≤L{Π}−μL{N(t)},(7)

N(s)≤Πμs−Πμ(s+μ)+N0(s+μ).(8)

The Laplace inverse transformation yield us

L−1{N(s)}≤L−1{Πμs}−L−1{Πμ(s+μ)}+L−1{N0(s+μ)},(9)

N(t)≤Πμ−e−μt(Πμ−N0).(10)

Thus,

limt→∞N(t)≤Πμ.(11)

Thus, y(t) is a bounded solution ∀ t≥0. ◼

Theorem 3.2. With the positive initial conditions, the co-infection model (Eq. (2)) has a positive solution y(t)=(S(t),E(t),II(t),IC(t),IIC(t),T(t),H(t),R(t)) for all t≥0.

Proof. Since, y(0) represent the total population at t=0, hence, y(0)≥0. Let us prove the positivity of one state variable and the left can be proved similarly. Let us take Eq. (2a)

dSdt=Π−(α1II+α2IC+α3IIC+μ)S.

As we already proved all the state variables are bounded by Πμ, therefore, ∃ χ>0 such that

χ=sup [α1II+α2IC+α3IIC+μ].

Thus,

dSdt≥Π−χS(t).

Applying the Laplace transformation gives

L{dSdt}≥ L{Π−χS(t)},

S(s)≥Πs(s+χ)+S0(s+χ).

Now, using the inverse Laplace transformation

L−1{S(s)}≥L−1{Πs(s+χ)}+L−1{S0(s+χ)},

S(t)≥Πχ−Πχe−χt+S0e−χt.(12)

Since 0≤e−χt≤1,and S0e−χt>0. Thus, it is self evident that S(t) ≥ 0 ∀ t≥0. Furthermore, we can verify this for other state variables. Thus, the feasible region of the proposed influenza-corona co-infection model (Eq. (2)) is given by

Ψ={(S(t),E(t),II(t),IC(t),IIC(t),T(t),H(t),R(t))∈R+8:0<N(t)≤Πμ,∀t≥ 0}.(13)

◼

Before proving the existence and uniqueness properties of the solution of the proposed model, we define following fundamental theorems:

Theorem 3.3. [27] Suppose 𝒮={(t,u1,u2,u3,…,um)|a≤t≤b,−∞<ui<∞,for each i=1,2,3…,m}, and let fi(t,u1,u2,u3,…,um) for each i = 1, 2, 3…, m be continuous on 𝒮 and satisfies a Lipchitz condition there. The system of first order differential equations duidt=fi(t,ui) with ui(a)=ωi for each i = 1, 2, 3…, m, has a unique solution for all a≤t≤b.

Theorem 3.4. The model (Eq. (3)) has a unique solution on [0,T] if F(y(t)) is Lipschitz continuous on [0,T].

Proof. Since, all the state variables are assumed to be continuously differentiable on [0,T], therefore, F(y(t)) and first partial derivatives of F(y(t)) are also continuous on [0,T]. Theorems 3.1 and 3.2 guarantee that F(y(t)) is bounded on [0,T], i.e., 0<F(y(t))≤Πμ. Thus, if the partial derivatives of the F(y(t)) are bounded, then Theorem 3.3 implies existence of a unique solution to the problem (Eq. (2)).

To prove Lipschitz condition, we need to prove that F(y(t)) has partial derivatives with respect to all state variables and are bounded.

|∂F1∂S|=|−(α1II+α2IC+α3IIC+μ)|≤Πμ<∞, |∂F1∂E|=0, |∂F1∂II|=|−α1S|≤Πμ, |∂F1∂IC|=|−α2S|≤Πμ, |∂F1∂IIC|=|−α3S|≤Πμ, |∂F1∂T|=|∂F1∂H|=|∂F1∂R|=0, |∂F2∂S|=|α1II+α2IC+α3IIC|≤Πμ<∞, |∂F2∂E|=|β1+β2+β3+μ|, |∂F2∂II|=|α1S|≤Πμ,

|∂F2∂IC|=|α2S|≤Πμ, |∂F2∂IIC|=|α3S|≤Πμ, |∂F2∂T|=|∂F2∂H|=|∂F2∂R|=0, |∂F3∂S|=0, |∂F3∂E|=β1,|∂F3∂II|=|(γ1+γ2+γ3+μ+d1)|, |∂F3∂IC|=|∂F3∂IIC|=|∂F3∂T|=|∂F3∂H|=|∂F3∂R|=0, |∂F4∂S|=0,|∂F4∂E|=β2, |∂F4∂II|=0, |∂F4∂IC|=|−(τ1+τ2+τc+μ+d2)|, |∂F4∂IIC|=|∂F4∂T|=|∂F4∂H|=|∂F4∂R|=0, |∂F5∂S|=0, |∂F5∂E|=β3, |∂F5∂II|=γ1, |∂F5∂IC|=|τ1|, |∂F5∂IIC|=|−(δ1+δ2+μ+d3)|,|∂F5∂T|=|∂F5∂H|=|∂F5∂R|=0, |∂F6∂S|=|∂F6∂E|=0, |∂F6∂II|=|γ2|, |∂F6∂IC|=|τ2|, |∂F6∂IIC|=|δ1|, |∂F6∂T|=|−(ϕ1+ϕ2+μ+d4)|, |∂F6∂H|=|∂F6∂R|=0, |∂F7∂S|=|∂F7∂E|=|∂F7∂II|=|∂F7∂IC|=|∂F7∂IIC|=0, |∂F7∂T|=|ϕ1|, |∂F7∂H|=|−(ϕ3+μ+d5)|, |∂F7∂R|=0, |∂F8∂S|=0|∂F8∂E|=0, |∂F8∂II|=|γ3|, |∂F8∂IC|=|τc|, |∂F8∂IIC|=|δ2|, |∂F8∂T|=|ϕ2|, |∂F8∂H|=|ϕ3|, |∂F8∂R|=|−μ|.

It is clear that all the partial derivatives, ∂Fi∂yj, i,j=1,2,3,…,8, exist and are bounded on [0,T]. Thus, Theorem 3.3 guarantees that the model (Eq. (3)) has a unique solution. ◼

Thus, we have proved that the proposed model (Eq. (2)) has a unique, positive, and bounded solution (Fig. 2).

Figure 2: Numerical simulations illustrating theoretical results proved analytically, that our proposed model is bounded and positive within a feasible region and has a unique solution

4  Equilibrium Points and Associated Properties

In this section, we determine the important equilibrium points and the properties associated with these equilibrium points. Specifically, we drive the reproduction number ℛ0 and perform stability analysis of the system at the equilibrium points. All the analytical results of stability analysis are validated through numerical simulations and given at the end of each section (Figs. 3 and 4).

Figure 3: Numerical simulation illustrating the our proposed model (Eq. (2)) is stable at disease free equilibrium point when ℛ0<1 (proved in Theorem 4.2). We numerically solve the system with different initial values, and all the solution curves approaches same point which is the disease free equilibrium point

Figure 4: We numerically solve the system by using different initial values, but all the solutions approach same point which is the endemic point. These simulations demonstrate the state theorem (Theorem 4.4) that the proposed model (Eq. (2)) is stable at endemic equilibrium point (P∗∗)

4.1 Equilibrium Points and Reproduction Number

To find equilibrium points, we consider the rate of change of all the state variables equal to zero, i.e., dSdt=dEdt=dIIdt=dICdt=dIICdt=dTdt=dHdt=dRdt=0. There are two main equilibrium points that exist for an epidemic model, DFE point and EE point. To find DFE point of an epidemic model, we assume that there is no infection present, i.e., E=IC=II=IIC=0. Thus, the DFE point of the corona-influenza co-infection model (Eq. (2)) is given as:

P∗=(S∗,E∗,IC∗, II∗, IIC∗, T∗, H∗, R∗)=(Πμ, 0, 0, 0, 0, 0, 0, 0).(14)

A critical epidemiological metric that quantifies the average number of secondary infections produced by a single infectious person in a fully susceptible (disease free) population is the reproduction number ℛ0. The ℛ0 can be computed mathematically using the next-generation matrix method [28–30]. To calculate ℛ0, we sub-divide infection carrier compartments into ℱ and 𝒢:

ℱ=((α1II+α2IC+α3IIC)S000),𝒢=(k1E−β1E+k2II−β2E+k3IC−β3E−γ1II−τ1IC+k4IIC),

where, k1=β1+β2+β3+μ, k2=γ1+γ2+γ3+μ+d1, k3=τ1+τ2+τc+μ+d2, and k4=δ1+δ2+μ+d3.

The Jacobian of the matrices ℱ and 𝒢 evaluated at DFE point P∗ are given as:

F=(0α1Πμα2Πμα3Πμ000000000000),G=(k1000−β1k200−β20k30−β3−γ1−τ1k4).

Thus, the spectral radius of FG−1=ℛ0=ℛI+ℛC+ℛIC,

where

ℛI=Πα1β1μk1k2, ℛC=Πα2β2μk1k3, and ℛIC=α3Πμ(β3k1k4+β1γ1k1k2k4+β2τ1k1k2k3).

Here ℛI, ℛC, and ℛIC represent the reproduction number of influenza, corona and influenza-corona co-infection, respectively. Thus, the reproduction number of the whole model (Eq. (2)) is given as:

ℛ0=Πμ[α1β1k3k4+α2β2k2k4+α3β3k2k3+α3γ1β1k3+α3τ1β2k2k1k2k3k4].(15)

Local Stability at DFE

Theorem 4.1. The proposed model (Eq. (2)) is locally asymptotically stable (LAS) at the disease free equilibrium point P∗ if ℛ0<1 and unstable for ℛ0>1.

Proof. This theorem can be proved by using the Jacobian matrix approach. For this purpose, we need to find Jacobian matrix at DFE point, i.e.,

J(P∗)=(−μ0−α1Πμ−α2Πμ−α3Πμ0000−k1α1Πμα2Πμα3Πμ0000β1−k2000000β20−k300000β3γ1τ1−k400000γ2τ2δ1−k50000000ϕ1−k6000γ3τcδ2ϕ2ϕ3−μ).(16)

As we know that, a system is said to be stable at an equilibrium if all the eigenvalues are negative at this equilibrium point and unstable otherwise. Therefore, we calculate the eigenvalues of the above matrix with the help of Maple software:

λ1∗=−μ,(17a)

λ2∗=−μ,(17b)

λ3∗=−k6,(17c)

λ4∗=−k5,(17d)

λ5∗=−k1,(17e)

λ6∗=−k2(1−β1α1Πμk1k2),(17f)

λ7∗=[(1−ℛ0)+(α3β3k2k3+α3β3k2k3+α3γ1β1k3+α3τ1β2k2)k1k4μλ6∗],(17g)

λ8∗=−[1−ℛ0]k1μλ6∗λ7∗.(17h)

Since all the eigenvalues are negatives when ℛ0<1, β1α1Πμk1k2<1 (true) and at least Eq. (17h) is positive otherwise. Thus, our proposed co-infection model (Eq. (2)) is locally asymptotically stable at DFE point (Eq. (14)). ◼

4.2 Global Stability

In order to show that the DFE point P∗ is globally stable, we use the approach given by Castillo-Chavez [31].

Theorem 4.2. The disease free equilibrium (DFE) point P∗ of the model (Eq. (2)) is globally asymptotically stable (GAS) if ℛ0<1 and the conditions (ℋ1) and (ℋ2) are fulfilled:

(ℋ1)dXdt=K(X,0)=0, X0 is GAS,(ℋ2)dYdt=N(X,Y)=BY−N¯(X,Y),

where N¯(X,Y)≥0  ∀ t≥0 and B=DYN(X0,0) is an M-matrix.

Proof. Let X=(S) represents non-infectious class and Y=(E,II,IC,IIC,T) represent infectious classes, and P∗=(X∗,0) is the DFE point. So,

dXdt=K(X,Y)=Π−(α1II+α2IC+α3IIC)S−μS.(18)

If X=X∗, then K(X,0)=0, i.e.,

dXdt=Π−μS∗=0.(19)

From Eq. (19), as t→∞, X→X∗. Therefore, X∗ is GAS. Now,

dYdt=BY−N¯(X,Y)=[−k1α1Πμα2Πμα3Πμ00β1−k20000β20−k3000β3γ1τ1−k4000γ2τ2δ1−k500000ϕ1−k6][EIIICIICTH]−[κ00000],(20)

where

B=[−k1α1Πμα2Πμα3Πμ00β1−k20000β20−k3000β3γ1τ1−k4000γ2τ2δ1−k500000ϕ1−k6],Y=[EIIICIICTH],N¯(X,Y)=[κ00000],

and κ=(α1II+α2IC+α3IIC)(S∗−S).

It is evident that B is an M-matrix. Since at the DFE point N=S∗, matrix N¯(X,Y)≥0. So, the DFE point P∗ is GAS. ◼

4.2.1 Endemic Equilibrium Point

To find endemic equilibrium point (P∗∗), we set E≠0,II≠0,IC≠0 and IIC≠0 in steady state equation of model (Eq. (2)).

P∗∗=(S∗∗,E∗∗,IC∗∗, II∗∗, IIC∗∗, T∗∗, H∗∗, R∗∗),(21)

where

S∗∗=Πμk1ℛ0,E∗∗=Π(ℛ0−k1k1ℛ0),II∗∗=β1k2E∗∗,IC∗∗=β2k3E∗∗,IIC∗∗=β3k2k3+γ1β1k3+τ1β2k2k2k3k4E∗∗,T∗∗=γ2II∗∗+τ2IC∗∗+δ1IIC∗∗k5,H∗∗=ϕ1T∗∗k6,R∗∗=γ3II∗∗+τcIC∗∗+δ2IIC∗∗+ϕ2T∗∗+ϕ3H∗∗μ.

4.2.2 Local Stability at EE

Theorem 4.3. The Model (Eq. (2)) is locally asymptotically stable (LAS) at endemic equilibrium point P∗∗ if ℛ0>1 and is unstable for ℛ0<1.

Proof. Now, we linearized the proposed model (Eq. (2)) at endemic equilibrium point (P∗∗), i.e.,

J(P∗∗)=(−(μ+α1II∗∗+α2IC∗∗+α3IIC∗∗)0−α1ℛ0−α2ℛ0−α3ℛ0000α1II∗∗+α2IC∗∗+α3IIC∗∗−k1α1ℛ0α2ℛ0α3ℛ00000β1−k2000000β10−k300000β3γ1τ1−k400000γ2τ2δ1−k50000000ϕ1−k6000γ3τcδ2ϕ2ϕ3−μ).

λ1∗∗=−μ,(22a)

λ2∗∗=−k5,(22b)

λ3∗∗=−k6,(22c)

λ4∗∗=−μℛ0,(22d)

λ5∗∗=−k1,(22e)

λ6∗∗=−k2(1−α1β1Πℛ02k1k2μ),(22f)

λ7∗∗=−k2k3λ6∗∗ℛ0[1−ℛ0−α3β3k2k3+α3γ1β1k3+α3τ1β2k2ℛ0k1k2k3k4μ],(22g)

λ8∗∗=−(k2k3k4λ6∗∗λ7∗∗ℛ0)(ℛ0−1).(22h)

Since all the eigenvalues are negative when ℛ0>1, α1β1Πℛ02k1k2μ<1 (true) and unstable otherwise, the proposed co-infection model is LAS at EE point (P∗∗). ◼

4.2.3 Global Stability at EE

Theorem 4.4. The endemic equilibrium (EE) point P∗∗ of the model (Eq. (2)) is globally asymptotically stable (GAS) provided ℛ0>1 and unstable when ℛ0<1.

Proof. Let ℒ be a Volterra-type Lyapunov function and defined as:

ℒ(S,E,II,IC,IIC,T,H,R)=[S−S∗∗−S∗∗log⁡SS∗∗]+[E−E∗∗−E∗∗log⁡EE∗∗]+[II−II∗∗−II∗∗log⁡IIII∗∗]+[IC−IC∗∗−IC∗∗log⁡ICIC∗∗]+[IIC−IIC∗∗−IIC∗∗log⁡IICIIC∗∗]+[T−T∗∗−T∗∗log⁡TT∗∗]+[H−H∗∗−H∗∗log⁡HH∗∗]+[R−R∗∗−R∗∗log⁡RR∗∗].(23)