# Dynamical Analysis of Secondary Dengue Viral Infection with Multiple Target Cells and Diffusion by Mathematical Model

Authors: Aeshah A. Raezah
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

Dengue is an epidemic disease rapidly spreading throughout many parts of the world, which is a serious public health concern. Understanding disease mechanisms through mathematical modeling is one of the most effective tools for this purpose. The aim of this manuscript is to develop and analyze a dynamical system of PDEs that describes the secondary infection caused by DENV, considering (i) the diffusion due to spatial mobility of cells and DENV particles, (ii) the interactions between multiple target cells, DENV, and antibodies of two types (heterologous and homologous). Global existence, positivity, and boundedness are proved for the system with homogeneous Neumann boundary conditions. Three threshold parameters are computed to characterize the existence and stability conditions of the model’s four steady states. Via means of Lyapunov functional, the global stability of all steady states is carried out. Our results show that the uninfected steady state is globally asymptotically stable if the basic reproduction number is less than or equal to unity, which leads to the disappearance of the disease from the body. When the basic reproduction number is greater than unity, the disease persists in the body with an active or inactive immune antibody response. To demonstrate such theoretical results, numerical simulations are presented.


---

## Body

## 1. Introduction

There are several types of viruses that infect humans and cause many diseases, such as HIV (human immunodeficiency virus), HCV (hepatitis C virus), HBV (hepatitis B virus), HTLV-I (human T lymphotropic virus type I), HPV (Human papilloma virus), Impetigo, EVD (Ebola virus), Chicken pox virus, Monkeypox virus, Rotavirus, Poliomyelitis, Measles, DENV (dengue virus), and COVID-19.The infection with dengue is induced by a mosquito-borne virus called dengue virus (DENV), which belongs to the family Flaviviridae. Dengue first appeared in the Philippines and Thailand during the 1950s. In 1970, only 9 countries experienced a severe dengue epidemic. Now, the disease is endemic in over 100 countries in the WHO regions of Africa, the Americas, the Eastern Mediterranean, South-East Asia, and the Western Pacific. DENV is transmitted to humans by mosquito bites. Aedes aegypti mosquitoes and Aedes albopictus mosquitoes are the common vectors of transmission [1, 2]. Dengue is most often asymptomatic (people may not even realize they are sick at all) or may have mild symptoms, but it can also cause a severe flu-like illness that often occurs with symptoms such as high fever, headache, rashes, nausea, leukopenia, joint pain, pain in the bones or muscles, severe bleeding, plasma leakage and/or organ impairment. In most cases, symptoms last 2–7 days, after the incubation period has passed for 4–10 days after mosquito bites [3, 4]. Many studies have indicated that the types of cells which can be infected by DENV are monocytes, hepatocytes, dendritic cells macrophages, and mast cells [5–27]. Several aspects of the immune system play a crucial role in controlling and eliminating DENV infection, including two important immune responses: cytotoxic T lymphocytes (CTL) and humoral immunity [6].


### 1.1. Primary and Secondary DENV Dynamics Models

Dengue virus (DENV) can cause dengue fever in four distinct but genetically related serotypes (DENV1-4). There is lifelong immunity to one serotype of dengue virus, but not cross-protection to other serotypes and an individual may become infected four times in a lifetime [2, 28, 29]. Obtaining a second infection with DENV can lead to DSS (dengue shock syndrome) or DHF (dengue hemorrhagic fever), which occurs after being infected again with a DENV serotype different from the first time [30]. In this case, the infection is called Secondary infection and the most common symptoms of DHF and DSS are thrombocytopenia, hemorrhage, and plasma leakage [27]. Recently, infection caused by DENV has been the subject of numerous mathematical studies.Primary DENV dynamics models: During primary infections, the body produces antibodies that neutralize the specific serotype of the infection (type-specific responses). Several mathematical models that describe Primary DENV infection within hosts have been studied (see [8–18]). Some related works regarding the immune responses have been investigated in primary DENV infections, including: cytotoxic T lymphocytes (CTLs) responses has been introduced by [8, 11, 14], a study of both humoral and innate immune responses has been published by [16], humoral immune response [30], Models of both humoral and CTL responses [12, 13, 17, 31]. In [30], Gujarati and Ambika introduced a mathematical model which describes the primary DENV infection with humoral immunity as follows:(1)dStdt=θ−βStVt−ξSt,dItdt=βStVt−ϱIt,dVtdt=τIt−ηVt−ψ1VtHt,dUtdt=μ−γUt+φVtUt,dHtdt=ϑUt−ϰ1VtHt−ϕ1Ht.For timet>0, the variable St is the concentration of healthy monocytes which are created at rate θ, and die at rate ξSt. It denotes DENV-infected monocytes which die at rate ϱIt and produce DENV particles Vt at the rate τIt. St become DENV-infected at rate βStVt. Vt are cleared at a rate ηVt. Ut and Ht represent B cells and neutralizing antibodies, respectively. The term ψ1VtHt denotes the neutralization rate of the DENV by antibodies, and μ,γUt,φVtUt are the production rate, die rate and stimulation rate of B-cells, respectively. ϑUt gives the rate of antibody production from the B-cells. Antibody virus complexes affect antibody production at the rate ϰ1VtHt, and the antibodies die at the rate ϕ1Ht.Secondary DENV dynamics model: During secondary infections, the body produces antibodies that neutralize two or more serotypes. Several related studies regarding secondary DENV infections are discussed in [8–27, 30, 32]. There have been various immune responses studied in secondary DENV infections, among them: secondary DENV infection model with humoral immunity by [20, 24, 30], both humoral and CTL immunity [22, 27], the mobility of cells and viruses has been studied by [25], Alshaikh et al. [26] developed a model of secondary DENV infection that takes into account that the DENV can infect multiple classes cell types, Alves Rubio et al. [32] studied and analyzed a mathematical model to explain antibody-mediated enhancement in heterologous of secondary infections and assess the effect of limiting cloning for plasma cells. The following model is a secondary DENV infection model with humoral immunity, which was introduced by [30]:(2)dStdt=θ−βStVt−ξSt,dItdt=βStVt−ϱIt,dVtdt=τIt−ηVt−ψ1VtHt−ψ2VtWt,dUtdt=μ−γUt+φVtUt,dHtdt=ϵϑVt−ϰ1VtHt−ϕ1Ht,dWtdt=ϑUt−ϰ2VtWt−ϕ2Wt.Here, at timet>0, Ht denotes the concentration of heterologous antibodies previously produced during the primary infection, and Wt is the concentration of homologous antibody specific to the new DENV serotype of the secondary infection. Neutralization rate of the DENV by antibodies Wt is represented by the term ψ2VtWt. Wt die at rate ϕ2Wt. The antibodies Ht and Wt are generated by B-cells at rates ϵϑUt and ϑUt, respectively, where ϵ∈0,1 is a correlation factor used to measure the similarity of one serotype to another. Antibody growth is affected by the antibody virus complex at rates ϰ1VtHt and ϰ2VtWt. All other variables and parameters of the model have the same biological concept as explained in the model (1). All parameters of models (1) and (2) are positive.It has been noteworthy that most DENV infection models in the above works did not consider the spatial structure (diffusion), as the models were presented by ordinary or delayed differential equations assuming that all components (DENV and cells) are homogeneous. While there have been other studies of various viral infections that have considered the spatial effects of various viral infections. They hypothesized that the spatial distributions of cells and virus particles are not fully mixed (see, e.g. [25, 33–38]). Among the studies that introduced DENV infection, taking into account diffusion, was Elaiw and Alofi study [25], they modeled a DENV infection with diffusion, by taking into consideration the spatial mobility of cells and DENV particles. Furthermore, [25] considered the DENV can infect one target of cells, but there are several works that have indicated that the dengue virus can infect multiple types of cells (e.g., [8–27]).In view of the foregoing, the DENV infection with both multiple target cells and diffusion did not a subject of mathematical studies before. Therefore, in this work, we present a model to study and analyze secondary DENV infection. Inspired by the study of Elaiw and Alofi [25] and Alshaikh et al. [26], we focus on that DENV can infect multiple target cells considering the mobility of every component of the model.


## 1.1. Primary and Secondary DENV Dynamics Models

Dengue virus (DENV) can cause dengue fever in four distinct but genetically related serotypes (DENV1-4). There is lifelong immunity to one serotype of dengue virus, but not cross-protection to other serotypes and an individual may become infected four times in a lifetime [2, 28, 29]. Obtaining a second infection with DENV can lead to DSS (dengue shock syndrome) or DHF (dengue hemorrhagic fever), which occurs after being infected again with a DENV serotype different from the first time [30]. In this case, the infection is called Secondary infection and the most common symptoms of DHF and DSS are thrombocytopenia, hemorrhage, and plasma leakage [27]. Recently, infection caused by DENV has been the subject of numerous mathematical studies.Primary DENV dynamics models: During primary infections, the body produces antibodies that neutralize the specific serotype of the infection (type-specific responses). Several mathematical models that describe Primary DENV infection within hosts have been studied (see [8–18]). Some related works regarding the immune responses have been investigated in primary DENV infections, including: cytotoxic T lymphocytes (CTLs) responses has been introduced by [8, 11, 14], a study of both humoral and innate immune responses has been published by [16], humoral immune response [30], Models of both humoral and CTL responses [12, 13, 17, 31]. In [30], Gujarati and Ambika introduced a mathematical model which describes the primary DENV infection with humoral immunity as follows:(1)dStdt=θ−βStVt−ξSt,dItdt=βStVt−ϱIt,dVtdt=τIt−ηVt−ψ1VtHt,dUtdt=μ−γUt+φVtUt,dHtdt=ϑUt−ϰ1VtHt−ϕ1Ht.For timet>0, the variable St is the concentration of healthy monocytes which are created at rate θ, and die at rate ξSt. It denotes DENV-infected monocytes which die at rate ϱIt and produce DENV particles Vt at the rate τIt. St become DENV-infected at rate βStVt. Vt are cleared at a rate ηVt. Ut and Ht represent B cells and neutralizing antibodies, respectively. The term ψ1VtHt denotes the neutralization rate of the DENV by antibodies, and μ,γUt,φVtUt are the production rate, die rate and stimulation rate of B-cells, respectively. ϑUt gives the rate of antibody production from the B-cells. Antibody virus complexes affect antibody production at the rate ϰ1VtHt, and the antibodies die at the rate ϕ1Ht.Secondary DENV dynamics model: During secondary infections, the body produces antibodies that neutralize two or more serotypes. Several related studies regarding secondary DENV infections are discussed in [8–27, 30, 32]. There have been various immune responses studied in secondary DENV infections, among them: secondary DENV infection model with humoral immunity by [20, 24, 30], both humoral and CTL immunity [22, 27], the mobility of cells and viruses has been studied by [25], Alshaikh et al. [26] developed a model of secondary DENV infection that takes into account that the DENV can infect multiple classes cell types, Alves Rubio et al. [32] studied and analyzed a mathematical model to explain antibody-mediated enhancement in heterologous of secondary infections and assess the effect of limiting cloning for plasma cells. The following model is a secondary DENV infection model with humoral immunity, which was introduced by [30]:(2)dStdt=θ−βStVt−ξSt,dItdt=βStVt−ϱIt,dVtdt=τIt−ηVt−ψ1VtHt−ψ2VtWt,dUtdt=μ−γUt+φVtUt,dHtdt=ϵϑVt−ϰ1VtHt−ϕ1Ht,dWtdt=ϑUt−ϰ2VtWt−ϕ2Wt.Here, at timet>0, Ht denotes the concentration of heterologous antibodies previously produced during the primary infection, and Wt is the concentration of homologous antibody specific to the new DENV serotype of the secondary infection. Neutralization rate of the DENV by antibodies Wt is represented by the term ψ2VtWt. Wt die at rate ϕ2Wt. The antibodies Ht and Wt are generated by B-cells at rates ϵϑUt and ϑUt, respectively, where ϵ∈0,1 is a correlation factor used to measure the similarity of one serotype to another. Antibody growth is affected by the antibody virus complex at rates ϰ1VtHt and ϰ2VtWt. All other variables and parameters of the model have the same biological concept as explained in the model (1). All parameters of models (1) and (2) are positive.It has been noteworthy that most DENV infection models in the above works did not consider the spatial structure (diffusion), as the models were presented by ordinary or delayed differential equations assuming that all components (DENV and cells) are homogeneous. While there have been other studies of various viral infections that have considered the spatial effects of various viral infections. They hypothesized that the spatial distributions of cells and virus particles are not fully mixed (see, e.g. [25, 33–38]). Among the studies that introduced DENV infection, taking into account diffusion, was Elaiw and Alofi study [25], they modeled a DENV infection with diffusion, by taking into consideration the spatial mobility of cells and DENV particles. Furthermore, [25] considered the DENV can infect one target of cells, but there are several works that have indicated that the dengue virus can infect multiple types of cells (e.g., [8–27]).In view of the foregoing, the DENV infection with both multiple target cells and diffusion did not a subject of mathematical studies before. Therefore, in this work, we present a model to study and analyze secondary DENV infection. Inspired by the study of Elaiw and Alofi [25] and Alshaikh et al. [26], we focus on that DENV can infect multiple target cells considering the mobility of every component of the model.


## 2. Model Formulation

As a result of the biological and mathematical considerations, we propose DENV infection model considering a system of multiple target cellsSj with two different types of antibodies (H1 and H2) and diffusion as follows:(3)∂Sjx,t∂t=djΔSjx,t+θj−βjSjx,tVx,t−ξjSjx,t,j=1,2,…,n,∂Ijx,t∂t=dn+jΔIjx,t+βjSjx,tVx,t−ϱjIjx,t,j=1,2,…,n,∂Vx,t∂t=d2n+1ΔVx,t+∑j=1nτjIjx,t−ηVx,t−ψ1Vx,tH1x,t−ψ2Vx,tH2x,t,∂H1x,t∂t=d2n+2ΔH1x,t+μ1Vx,tH1x,t−ϕ1H1x,t,∂H2x,t∂t=d2n+3ΔH2x,t+μ2Vx,tH2x,t−ϕ2H2x,t,where x∈Λ is the position. μ1VH1 and μ2VH2 are the reproduction rates of heterologous and homologous antibodies. The spatial domain Λ⊂ℝk (where k≥1) is connected and bounded, and the boundary of ∂Λ is smooth. Finally, Δ is the Laplacian operator and di denotes the diffusion coefficient of cells and DENV particles, where i∈1,2,…,n,…,2n+3.For model (3), we consider the following homogeneous Neumann boundary conditions:(4)∂S1∂G⟶=…=∂Sn∂G⟶=∂I1∂G⟶=…=∂In∂G⟶=∂V∂G⟶=∂H1∂G⟶=∂H2∂G⟶=0,t>0,x∈∂Λ.And initial conditions(5)Sjx,0=ℬj1x,j=1,2,…,n,Ijx,0=ℬj2x,j=1,2,…,n,Vx,0=ℬ3x,H1x,0=ℬ4x,H2x,0=ℬ5x,x∈Λ¯,.where ∂/∂G⟶ is the outward normal derivative on the boundary ∂Λ, ℬj1x,ℬj2x,ℬlx≥0, ℓ=3,4,5, are continuous functions. The boundary conditions (4) point out that cells and DENVs cannot pass the isolated boundary [39].


## 3. Biologically Realistic Domain

To ensure that the model (3) is biologically acceptable, we define a bounded domain for the compartment concentrations.Proposition 1.
Assume that di=d˜ for all i=1,2,..,n,…,2n+3. Then, model (3) with any initial satisfying (5) has a unique, nonnegative, and bounded solution defined on Λ¯×0,+∞.Proof.
LetX=BUCΛ¯,ℝ2n+3 be the set of all bounded and uniformly continuous functions from Λ¯ to ℝ2n+3, with the norm ϖX=supx∈Λ¯ϖx, where ⋅ is the Euclidean norm on ℝ2n+3. We define the positive cone X+=BUCΛ¯,ℝ+2n+3⊂X which induces a partial order on X. Therefore, the space X,⋅X is a Banach lattice [40, 41].
Now, for any initial data(6)ℬ=ℬ11,…ℬn1,ℬ12,…,ℬn2,ℬ3,ℬ4,ℬ5T∈X+,
we defineG=G11,…,Gn1,G12,…,Gn2,G3,G4,G5T:X+⟶X by(7)Gj1ℬx=θj−βjℬj1xℬ3x−ξjℬj1x,j=1,2,…,n,Gj2ℬx=βjℬj1xℬ3x−ϱjℬj2x,j=1,2,…,n,G3ℬx=∑j=1nτjℬj2x−ηℬ3x−ψ1ℬ3xℬ4x−ψ2ℬ3xℬ5x,G4ℬx=μ1ℬ3xℬ4x−ϕ1ℬ4x,G5ℬx=μ2ℬ3xℬ5x−ϕ2ℬ5x.
Clearly,ℬ is locally Lipschitz on X+. With boundary conditions (4) and initial conditions (5), the system (3) can be rewritten as an abstract functional differential equation with the following formula(8)dU¯dt=DU¯+GU¯,t>0,U¯0=ℬ∈X+,where(9)U¯=S1,…,Sn,I1,…,In,V,H1,H2T,and(10)DU¯=d1ΔS1,…,dnΔSn,dn+1ΔI1,…,d2nΔIn,d2n+1ΔV,d2n+2ΔH1,d2n+3ΔH2T.
As a result, it can be shown that(11)limq⟶0+1qdistℬ0+qGℬ,X+=0,∀ℬ∈X+.
Hence, for anyℬ∈X+, system (3) with conditions (4) and (5) has a single non-negative solution S1x,t,…,Snx,t,I1x,t,…,Inx,t,Vx,t,H1x,t,H2x,t defined on Λ¯×0,Tk, where the interval 0,Tk is the maximal existence time period where the solution exists [40–42]. Moreover, for the given problem, this solution also constitutes a classical solution.
To prove boundedness, let(12)Tjx,t=Sjx,t+Ijx,t,
then(13)∂Tjx,t∂t=djΔSjx,t+θj+dn+jΔIjx,t−ξjSjx,t+ϱjIjx,t≤djΔSjx,t+θj+dn+jΔIjx,t−σjTjx,t,where σj=minξj,ϱj,j=1,…,n. Since dj=dn+j=d˜, we have(14)∂Tjx,t∂t≤θj+d˜ΔSjx,t+Ijx,t−σjTjx,t=θj+d˜ΔTjx,t−σjTjx,t,
then(15)∂Tjx,t∂t−d˜ΔTjx,t≤θj−σjTjx,t.
Thus,Tjx,t satisfies(16)∂Tjx,t∂t−d˜ΔTjx,t≤θj−σjTjx,t,j=1,2,…,n,Tjx,0=ℬj1x+ℬj2x≥0,j=1,2,…,n,∂Tj∂G⟶=0,j=1,2,…,n,
LetTj˜t be a solution of the following ordinary differential equation (ODE):(17)dTj˜tdt=θj−σjTj˜t,j=1,2,…,n,Tj˜0=maxx∈Λ¯Tjx,0,j=1,2,…,n.
This implies thatT˜jt≤maxθj/σj,maxx∈Λ¯Tjx,0. Based on the principle of comparison [43], we have Tjx,t≤Tj˜t. Hence(18)Tjx,t≤maxθjσj,maxx∈Λ¯Tjx,0=Kj.
Therefore,Sjx,t,Ijx,t are bounded on Λ¯×0,Tk.
Additionally, let(19)Mx,t=Vx,t+ψ1μ1H1x,t+ψ2μ2H2x,t,
then(20)∂Mx,t∂t=d2n+1ΔVx,t+∑j=1nτjIjx,t−ηVx,t+ψ1μ1d2n+2ΔH1x,t−ψ1μ1ϕ1H1x,t+ψ2μ2d2n+3ΔH2x,t−ψ2μ2ϕ2H2x,t≤d2n+1ΔVx,t+ψ1μ1d2n+2ΔH1x,t+ψ2μ2d2n+3ΔH2x,t+∑j=1nτjIjx,t−αMx,t,
withα=minη,ϕ1,ϕ2. Since d2n+1=d2n+2=d2n+3=d˜, and Ijx,t≤Kj, we obtain(21)∂Mx,t∂t−d˜ΔMx,t≤∑j=1nτjKj−αMx,t.
Thus,Mx,t satisfies(22)∂Mx,t∂t−d˜ΔMx,t≤K−αMx,t,Mx,0=ℬ3x+ψ1μ1ℬ4x+ψ2μ2ℬ5x≥0,∂M∂G⟶=0,where K=∑j=1nτjKj. Suppose the following ODE(23)dM˜tdt=K−αM˜t,M˜0=maxx∈Λ¯Mx,0,
has a solutionM˜t, this gives that M˜t≤maxK/α,maxx∈Λ¯Mx,0. As above, we obtain Mx,t≤M˜t. Thus,(24)Mx,t≤maxKα,maxx∈Λ¯Mx,0.
Therefore, it impliesVx,t,H1x,t, and H2x,t are bounded on Λ¯×0,Tm.
According to semilinear parabolic system (standard theory);Tm=+∞ [44]. This proves that the solution S1x,t,…,Snx,t,I1x,t,…,Inx,t,Vx,t,H1x,t,H2x,t is defined for all x∈Λ, t>0 and it is unique and nonnegative.


## 4. Steady States and Thresholds

Through this section, we deduce three threshold parameters which define the existence of all possible steady states of model (3). To calculate the steady states, we solve the following system of algebraic equations(25)0=θj−βjSjV−ξjSj,j=1,2,…,n,(26)0=βjSjV−ϱjIj,j=1,2,…,n,(27)0=∑j=1nτjIj−ηV−ψ1VH1−ψ2H2,(28)0=μ1VH1−ϕ1H1,(29)0=μ2VH2−ϕ2H2,and then, we have the following:(i)
IfH1=H2=0,V=0, there is an uninfected steady state Π0S10,…,Sn0,0,…,0,0︷n+3, with Sj0=θj/ξj.(ii)
IfH1=H2=0,V≠0, there is an infected steady state with inactive immune antibody response Π1S˜1,…,S˜n,I˜1,…,I˜n,V˜,0,0, where(30)S˜j=θjξj+βjV˜,I˜j=βjS˜jV˜ϱj,j=1,2,…,n,andV˜ satisfies(31)∑j=1nτjβjS˜jηϱj−1ηV˜=0.To determine suchV˜, in case of the equilibrium Π1, we have V≠0, then(32)∑j=1nτjβjSjηϱj−1=0.Substituting from equation (25) into equation (32) we get(33)∑j=1nτjθj−ξjSjηϱj−V=0,whereSj=SjV. From equation (33), we can define a function ϜV as:(34)ϜV=∑j=1nτjηϱjθj−ξjθjξj+βjV−V.Now, we need to show that∃V˜>0 such that ϜV=0. Note that(35)Ϝ0=0,Ϝ′0=∑j=1nτjβjθjηξjϱj−1>0if∑j=1nτjβjθjηξjϱj>1,ϜV∗=−∑j=1nτjξjSj∗ηϱj<0, whereV∗=∑j=1nτjθjηϱj>0andSjV∗=Sj∗,Then, if∑j=1nτjβjθj/ηξjϱj>1, ∃V˜∈0,V∗ such that ϜV˜=0. This yields that the steady state Π1S˜1,…,S˜n,I˜1,…,I˜n,V˜,0,0 exists when ∑j=1nτjβjθj/ηξjϱj>1. Therefore, the basic reproduction number of our model (3), is given by(36)ℜ0=∑j=1nτjβjθjηξjϱj.ℜ0 determines whether the infection is present in the host.(iii)
IfH1≠0,H2=0, an infected steady state with only active heterologous antibody Π2S¯1,…,S¯n,I¯1,…,I¯n,V¯,H¯1,0 is obtained, where(37)S¯j=θjξj+βjϕ1/μ1,I¯j=βjθjϕ1ϱjμ1ξj+βjϕ1,V¯=ϕ1μ1,H¯1=ηψ1ℜ1−1,and(38)ℜ1=∑j=1nμ1τjβjθjηϱjμ1ξj+βjϕ1=∑j=1nτjβjθjηϱjξj+βjV¯.It follows thatH¯1>0 only when ℜ1>1. In this case, ℜ1 is the threshold parameter that indicates the initiation of an immune reaction to a heterologous antibody.(iv)
IfH1=0,H2≠0, then an infected steady state with only active homologous antibody Π3S^1,…,S^n,I^1,…,I^n,V^,0,H^2 exists, where(39)S^j=θjξj+βjϕ2/μ2,I^j=βjθjϕ2ϱjμ2ξj+βjϕ2,V^=ϕ2μ2,H^2=ηψ2ℜ2−1,and(40)ℜ2=∑j=1nμ2τjβjθjηϱjμ2ξj+βjϕ2=∑j=1nτjβjθjηϱjξj+βjV^.Thus,H^2>0 only when ℜ2>1. Here, ℜ2 stands for the threshold parameter that indicates the initiation of an immune reaction to a homologous antibody.It is clear from the above that:ℜ0>ℜ1, ℜ0>ℜ2, and if ℜ1>1 and ℜ2>1, Π0, Π1, Π2 and Π3 all exist.Corollary 1.
Suppose that the conditionsℜ0>0,ℜ1>0, and ℜ2>0 exist. Then the model (3) has four steady states when the following conditions are attained:(a)
The uninfected steady stateΠ0 is always exist,(b)
The infected steady state with inactive immune antibody responseΠ1 exists if ℜ0>1,(c)
The infected steady state with only active heterologous antibodyΠ2 exists if ℜ1>1,(d)
The infected steady state with only active homologous antibodyΠ3 exists if ℜ2>1.


## 5. Global Stability

This section investigates the global stability of the steady states by using Lyapunov’s method. We construct Lyapunov functions using the same technique in [25, 26, 38].Before discussing theories of global stability, it is important to consider the following:(i)
We specifygw=w−1−lnw, where w>0, and gw=0⇔w=1.(ii)
For simplicity, we omit the notation for inputs. i.e.S1,…,Sn,I1,…,In,V,H1,H2=S1x,t,…,Snx,t,I1x,t,…,Inx,t,Vx,t,H1x,t,H2x,t, Lℓ=Lℓx,t, and L˜ℓ=L˜ℓt for all ℓ=0,1,2,3.(iii)
Numann boundary conditions (4) and Divergence Theorem imply that(41)0=∫∂Λ∇U⋅G⟶dx=∫Λdiv∇Udx=∫ΛΔUdx,0=∫∂Λ1U∇U⋅G⟶dx=∫Λdiv1U∇Udx=∫ΛΔUU−▽L2U2dx,forU∈S1,…,Sn,I1,…,In,V,H1,H2. Thus, we obtain(42)∫ΛΔUdu=0,∫ΛΔUUdx=∫Λ▽U2U2dx.(iv)
Consider a functionLℓ=LℓS1,⋯,Sn,I1,⋯,In,V,H1,H2 and define L˜ℓt such that(43)L˜ℓt=∫ΛLℓx,tdx,ℓ=0,1,2,3.(v)
LetΥℓ′ be the largest invariant subset of(44)Υℓ=S1,⋯,Sn,I1,⋯,In,V,H1,H2:dL˜ℓdt=0,ℓ=0,1,2,3.Theorem 1.
For model (3), suppose that ℜ0≤1, then steady state Π0 is globally asymptotically stable (G.A.S).Proof.
Define a Lyapunov functionL0 as follows:(45)L0=∑j=1nτjϱjSj0gSjSj0+Ij+V+ψ1μ1H1+ψ2μ2H2.
Clearly,L˜0S1,⋯,Sn,I1,⋯,In,V,H1,H2>0 for all S1,…,Sn,I1,…,In,V,H1,H2>0 and L˜0S10,⋯,Sn0,0,⋯,0,0︷n+3=0. The partial derivative ∂L0/∂t can be calculated along the solutions of model (3) as:(46)∂L0∂t=∑j=1nτjϱj1−Sj0SjdjΔSj+ξjSj0−βjSjV−ξjSj+dn+jΔIj+βjSjV−ϱjIj+d2n+1ΔV+∑j=1nτjIj−ηV−ψ1VH1−ψ2VH2+ψ1μ1d2n+2ΔH1+μ1VH1−ϕ1H1+ψ2μ2d2n+3ΔH2+μ2VH2−ϕ2H2.
After collecting the terms of equation (46), we have(47)∂L0∂t=−∑j=1nξjτjϱjSj−Sj02Sj+∑j=1nτjϱj1−Sj0SjdjΔSj+βjSj0V+dn+jΔIj+d2n+1ΔV−ηV+ψ1μ1d2n+2ΔH1−ψ1μ1ϕ1H1+ψ2μ2d2n+3ΔH2−ψ2μ2ϕ2H2=−∑j=1nξjτjϱjSj−Sj02Sj+ℜ0−1ηV−ψ1μ1ϕ1H1−ψ2μ2ϕ2H2+∑j=1nτjϱj1−Sj0SjdjΔSj+∑j=1nτjϱjdn+jΔIj+d2n+1ΔV+ψ1μ1d2n+2ΔH1+ψ2μ2d2n+3ΔH2.
Therefore, the derivativedL˜0/dt will be as follows:(48)dL˜0dt=−∑j=1nξjτjϱj∫ΛSj−Sj02Sjdx+ℜ0−1η∫ΛVdx−ψ1μ1ϕ1∫ΛH1dx−ψ2μ2ϕ2∫ΛH2dx+∑j=1nτjdjϱj∫Λ1−Sj0SjΔSjdx+∑j=1nτjϱjdn+j∫ΛΔIjdx+d2n+1∫ΛΔVdx+ψ1μ1d2n+2∫ΛΔH1dx+ψ2μ2d2n+3∫ΛΔH2dx.
Using equalities (42), (48) is transformed to(49)dL˜0dt=−∑j=1nξjτjϱj∫ΛSj−Sj02Sjdx+ℜ0−1η∫ΛVdx−ψ1μ1ϕ1∫ΛH1dx−ψ2μ2ϕ2∫ΛH2dx−∑j=1nτjdjSj0ϱj∫Λ▽Sj2Sj2dx.
Hence,dL˜0/dt≤0 for all S1,…,Sn,V,H1,H2>0. Additionally, dL˜0/dt=0 when S1,…,Sn,V,H1,H2=S10,…,Sn0,0,0,0. The solutions of model (3) converge to Υ0′. The elements of Υ0′ satisfy S1,…,Sn,V,H1,H2=S10,…,Sn0,0,0,0 and then ∂V/∂t=ΔV=0. Third equation of (3) reduces to(50)0=∂V∂t=∑j=1nτjIj.
This yieldsIjt=0. Hence, Υ0′=Π0 and by applying LaSalle’s invariance principle [45–47] we get that Π0 is G.A.S.Theorem 2.
For model (3), suppose that ℜ1≤1<ℜ0 and ℜ2≤1, then Π1 is G.A.S.Proof.
DefineL1 as:(51)L1=∑j=1nτjϱjS˜jgSjS˜j+I˜jgIjI˜j+V˜gVV˜+ψ1μ1H1+ψ2μ2H2.
Calculating∂L1/∂t as:(52)∂L1∂t=∑j=1nτjϱj1−S˜jSjdjΔSj+θj−βjSjV−ξjSj+1−I˜jIjdn+jΔIj+βjSjV−ϱjIj+1−V˜Vd2n+1ΔV+∑j=1nτjIj−ηV−ψ1VH1−ψ2VH2+ψ1μ1d2n+2ΔH1+μ1VH1−ϕ1H1+ψ2μ2d2n+3ΔH2+μ2VH2−ϕ2H2=∑j=1nτjϱj1−S˜jSjθj−ξjSj+∑j=1nτjϱjβjS˜jV−∑j=1nτjϱj1−I˜jIjϱjIj−∑j=1nτjϱjI˜jIjβjSjV+1−V˜V∑j=1nτjIj−ηV+ψ1V˜H1+ψ2V˜H2−ψ1ϕ1μ1H1−ψ2ϕ2μ2H2+∑j=1nτjϱj1−S˜jSjdjΔSj+∑j=1nτjϱj1−I˜jIjdn+jΔIj+1−V˜Vd2n+1ΔV+ψ1μ1d2n+2ΔH1+ψ2μ2d2n+3ΔH2.
The steady state conditions atΠ1 are given by(53)θj=βjS˜jV˜+ξjS˜j,βjS˜jV˜=ϱjI˜j,ηV˜=∑j=1nτjI˜j.
Using equation (53), equation (52) becomes(54)∂L1∂t=−∑j=1nξjτjϱjSj−S˜j2Sj+∑j=1nτj1−S˜jSjI˜j−∑j=1nτjI˜jI˜jSjVIjS˜jV˜+2∑j=1nτjI˜j−∑j=1nτjI˜jIjV˜I˜jV+ψ1V˜−ϕ1μ1H1+ψ2V˜−ϕ2μ2H2+∑j=1nτjϱj1−S˜jSjdjΔSj+∑j=1nτjϱj1−I˜jIjdn+jΔIj+1−V˜jVjd2n+1ΔV+ψ1μ1d2n+2ΔH1+ψ2μ2d2n+3ΔH2.
Simplify (54), we get(55)∂L1∂t=−∑j=1nξjτjϱjSj−S˜j2Sj+∑j=1nτjI˜j3−S˜jSj−I˜jSjVIjS˜jV˜−IjV˜I˜jV+ψ1V˜−V¯H1+ψ2V˜−V^H2+∑j=1nτjϱj1−S˜jSjdjΔSj+∑j=1nτjϱj1−I˜jIjdn+jΔIj+1−V˜jVjd2n+1ΔV+ψ1μ1d2n+2ΔH1+ψ2μ2d2n+3ΔH2.
Calculating the time derivative ofL˜1t, we get(56)dL˜1dt=−∑j=1nξjτjϱj∫ΛSj−S˜j2Sjdx+∑j=1nτjI˜j∫Λ3−S˜jSj−I˜jSjVIjS˜jV˜−IjV˜I˜jVdx+ψ1V˜−V¯∫ΛH1dx+ψ2V˜−V^∫ΛH2dx+∑j=1nτjϱjdj∫Λ1−S˜jSjΔSjdx+∑j=1nτjϱjdn+j∫Λ1−I˜jIjΔIjdx+d2n+1∫Λ1−V˜jVjΔVdx+ψ1μ1d2n+2∫ΛΔH1dx+ψ2μ2d2n+3∫ΛΔH2dx.
Using equalities (42) to get(57)dL˜1dt=−∑j=1nξjτjϱj∫ΛSj−S˜j2Sjdx+∑j=1nτjI˜j∫Λ3−S˜jSj−I˜jSjVIjS˜jV˜−IjV˜I˜jVdx+ψ1∫ΛV˜−V¯H1dx+ψ2∫ΛV˜−V^H2dx−∑j=1nτjϱjdjS˜j∫ΛΔSj2Sj2dx−∑j=1nτjϱjdn+jI˜j∫ΛΔIj2Ij2dx−d2n+1V˜j∫ΛΔV2Vj2dx.
It follows from the relation between arithmetic and geometric means that(58)3≤S˜jSj+I˜jSjVIjS˜jV˜+IjV˜I˜jV.
Next we prove thatV˜≤V¯ when ℜ1≤1. From (53) we have ∑j=1nτjβjθj/ηϱjξj+βjV˜=1, then(59)0≥ℜ1−1=∑j=1nτjβjθjηϱjξj+βjV¯−∑j=1nτjβjθjηϱjξj+βjV˜=∑j=1nτjβj2θjV˜−V¯ηϱjξj+βjV¯ξj+βjV˜=V˜−V¯χ˜,where χ˜=∑j=1nτjβj2θj/ηϱjξj+βjV¯ξj+βjV˜. Likewise, we can prove that V˜≤V^ when ℜ2≤1. Thus, it can be concluded that dL˜1/dt≤0 for all S1,⋯,Sn,I1,⋯,In,V,H1,H2>0 since ℜ1≤1 and ℜ2≤1. Moreover, dL˜1/dt=0 at S1,⋯,Sn,I1,⋯,In,V,H1,H2=S˜1,⋯,S˜n,I˜1,⋯,I˜n,V˜,0,0. The solutions of model (3) tend to Υ1′=Π1. It follows that Π1 is G.A.S as per LaSalle’s invariance principle [45–47].Theorem 3.
For model (3), suppose that ℜ1>1 and ℜ2≤ℜ1, then Π2 is G.A.S.Proof.
Define a functionL2 as:(60)L2=∑j=1nτjϱjS¯jgSjS¯j+I¯jgIjI¯j+V¯gVV¯+ψ1μ1H¯1gH1H¯1+ψ2μ2H2.
As a result of computing∂L2/∂t along the solutions of model (3), we obtain(61)∂L2∂t=∑j=1nτjϱj1−S¯jSjdjΔSj+θj−ξjSj−βjSjV+1−I¯jIjdn+jΔIj+βjSjV−ϱjIj+1−V¯Vd2n+1ΔV+∑j=1nτjIj−ηV−ψ1VH1−ψ2VH2+ψ1μ11−H¯1H1d2n+2ΔH1+μ1VH1−ϕ1H1+ψ2μ2d2n+3ΔH2+μ2VH2−ϕ2H2.
Equation (61) is simplified to(62)∂L2∂t=∑j=1nτjϱj1−S¯jSjθj−ξjSj+βjS¯jV−βjSjVI¯jIj+ϱjI¯j−V¯V∑j=1nτjIj−ηV+ηV¯+ψ1V¯H1+ψ2V¯H2−ψ1ϕ1μ1H1−ψ1H¯1V+ψ1ϕ1μ1H¯1−ψ2ϕ2μ2H2+∑j=1nτjϱj1−S¯jSjdjΔSj+∑j=1nτjϱj1−I¯jIjdn+jΔIj+1−V¯Vd2n+1ΔV+ψ1μ11−H¯1H1d2n+2ΔH1+ψ2μ2d2n+3ΔH2.
The steady state conditions atΠ2 are given by(63)θj=ξjS¯j+βjS¯jV¯,βjS¯jV¯=ϱjI¯j,ηV¯=∑j=1nτjI¯j−ψ1V¯H¯1,V¯=ϕ1μ1.
When applied the aforementioned conditions, equation (62) becomes(64)∂L2∂t=−∑j=1nτjξjϱjSj−S¯j2Sj+∑j=1nτj1−S¯jSjI¯j−∑j=1nτjI¯jSjVI¯jS¯jV¯Ij+2∑j=1nτjI¯j−∑j=1nτjI¯jIjV¯I¯jV+ψ2V¯−ϕ2μ2H2+∑j=1nτjϱj1−S¯jSjdjΔSj+∑j=1nτjϱj1−I¯jIjdn+jΔIj+1−V¯Vd2n+1ΔV+ψ1μ11−H¯1H1d2n+2ΔH1+ψ2μ2d2n+3ΔH2=−∑j=1nτjξjϱjSj−S¯j2Sj+∑j=1nτjI¯j3−S¯jSj−SjVI¯jS¯jV¯Ij−IjV¯I¯jV+ψ2V¯−ϕ2μ2H2+∑j=1nτjϱj1−S¯jSjdjΔSj+∑j=1nτjϱj1−I¯jIjdn+jΔIj+1−V¯Vd2n+1ΔV+ψ1μ11−H¯1H1d2n+2ΔH1+ψ2μ2d2n+3ΔH2.
CalculatingdL˜2/dt, and then using equalities (42), we get(65)dL˜2dt=−∑j=1nτjξjϱj∫ΛSj−S¯j2Sjdx+∑j=1nτjI¯j∫Λ3−S¯jSj−SjVI¯jS¯jV¯Ij−IjV¯I¯jVdx+ψ2V¯−V^∫ΛH2dx−∑j=1nτjϱjdjS¯j∫Λ▽Sj2Sj2dx−∑j=1nτjϱjdn+jI¯j∫Λ▽Ij2Ij2dx−d2n+1V¯∫Λ▽V2V2dx−ψ1μ1d2n+2H¯1∫Λ▽H12H12dx.
It follows from the relation between arithmetic and geometric means that(66)3≤S¯jSj+SjVI¯jS¯jV¯Ij+IjV¯I¯jV.
Next, we prove thatV¯≤V^ when ℜ2≤ℜ1. We have(67)ℜ2−ℜ1=∑j=1nτjβjθjηϱjξj+βjV^−∑j=1nτjβjθjηϱjξj+βjV¯=∑j=1nτjβj2θjV¯−V^ηϱjξj+βjV^ξj+βjV¯=V¯−V^χ¯,where χ¯=∑j=1nτjβj2θj/ηϱjξj+βjV^ξj+βjV¯. Thus, if ℜ2≤ℜ1, we obtain dL˜2/dt≤0 for all S1,…,Sn,I1,…,In,V,H1,H2>0. Furthermore, dL˜2/dt=0 at S1,…,Sn,I1,…,In,V,H1,H2=S¯1,…,S¯n,I¯1,…,I¯n,V¯,H¯1,0. The solutions of system (3) tend to Υ2′=Π2. It follows that Π2 is G.A.S as per LaSalle’s invariance principle [45–47].Theorem 4.
For model (3), suppose that ℜ2>1 and ℜ1≤ℜ2, then Π3 is G.A.S.Proof.
Consider a Lyapunov functionL3 as:(68)L3=∑j=1nτjϱjS^jgSjS^j+I^jgIjI^j+V^gVV^+ψ1μ1H1+ψ2μ2H^2gH2H^2.
Then, we obtain(69)∂L3∂t=∑j=1nτjϱj1−S^jSjdjΔSj+θj−ξjSj−βjSjV+1−I^jIjdn+jΔIj+βjSjV−ϱjIj+1−V^Vd2n+1ΔV+∑j=1nτjIj−ηV−ψ1VH1−ψ2VH2+ψ1μ1d2n+2ΔH1+μ1VH1−ϕ1H1+ψ2μ21−H^2H2d2n+3ΔH2+μ2VH2−ϕ2H2.
Equation (69) is simplified to(70)∂L3∂t=∑j=1nτjϱj1−S^jSjθj−ξjSj+βjS^jV−βjSjVI^jIj+ϱjI^j−ηV−V^V∑j=1nτjIj+ηV^+ψ1V^H1+ψ2V^H2−ψ1μ1ϕ1H1−ψ2μ2ϕ2H2−ψ2VH^2+ψ2μ2ϕ2H^2+∑j=1nτjϱj1−S^jSjdjΔSj+∑j=1nτjϱj1−I^jIjdn+jΔIj+1−V^Vd2n+1ΔV+ψ1μ1d2n+2ΔH1+ψ2μ21−H^2H2d2n+3ΔH2.
Applying the steady state conditions atΠ3:(71)θj=ξjS^j+βjS^jV^,βjS^jV^=ϱjI^j,ηV^=∑j=1nτjI^j−ψ2V^H^2,V^=ϕ2μ2,
we have(72)∂L3∂t=−∑j=1nτjξjϱjSj−S^j2Sj+∑j=1nτjI^j1−S^jSj−∑j=1nτjI^jSjVI^jS^jV^Ij+2∑j=1nτjI^j−V^V∑j=1nτjIj+ψ1V^−ϕ1μ1H1+∑j=1nτjϱj1−S^jSjdjΔSj+∑j=1nτjϱj1−I^jIjdn+jΔIj+1−V^Vd2n+1ΔV+ψ1μ1d2n+2ΔH1+ψ2μ21−H^2H2d2n+3ΔH2=−∑j=1nτjξjϱjSj−S^j2Sj+∑j=1nτjI^j3−S^jSj−SjVI^jS^jV^Ij−IjV^I^jV+ψ1V^−V¯H1+∑j=1nτjϱj1−S^jSjdjΔSj+∑j=1nτjϱj1−I^jIjdn+jΔIj+1−V^Vd2n+1ΔV+ψ1μ1d2n+2ΔH1+ψ2μ21−H^2H2d2n+3ΔH2.
CalculatingdL˜3/dt, we obtain(73)dL˜3dt=−∑j=1nτjξjϱj∫ΛSj−S^j2Sjdx+∑j=1nτjI^j∫Λ3−S^jSj−SjVI^jS^jV^Ij−IjV^I^jVdx+ψ1V^−V¯∫ΛH1dx+∑j=1nτjϱjdj∫Λ1−S^jSjΔSjdx+∑j=1nτjϱjdn+j∫Λ1−I^jIjΔIjdx+d2n+1∫Λ1−V^VΔVdx+ψ1μ1d2n+2∫ΛΔH1dx+ψ2μ2d2n+3∫Λ1−H^2H2ΔH2dx.
Using equalities (42), we get(74)dL˜3dt=−∑j=1nτjξjϱj∫ΛSj−S^j2Sjdx+∑j=1nτjI^j∫Λ3−S^jSj−SjVI^jS^jV^Ij−IjV^I^jVdx+ψ1V^−V¯∫ΛH1dx−∑j=1nτjϱjdjS^j∫Λ▽Sj2Sj2dx−∑j=1nτjϱjdn+jI^j∫Λ▽Ij2Ij2dx−d2n+1V^∫Λ▽V2V2dx−ψ2μ2d2n+3H^2∫Λ▽H22H22dx,where 3≤S^j/Sj+SjVI^j/S^jV^Ij+IjV^/I^jV, and from equation (67) we have V^≤V¯ if ℜ1≤ℜ2. Therefore, dL˜3/dt≤0 for all S1,…,Sn,I1,…,In,V,H1,H2>0. Additionally, dL˜3/dt=0 at S1,…,Sn,I1,…,In,V,H1,H2=S^1,…,S^n,I^1,…,I^n,V^,0,H^2. The solutions of (3) tend to Υ3′=Π3. It follows that Π3 is G.A.S based on LaSalle’s invariance principle [45–47].
The existence and global stability conditions of the steady statesΠ0,Π1,Π2, and Π3 are outlined in Table 1.Table 1
The conditions of existence and global stability for the steady states of model (3).
Steady stateConditions of existenceConditions of global stabilityΠ0S10,…,Sn0,0,…,0,0︷n+3Noneℜ0≤1Π1S˜1,⋯,S˜n,I˜1,⋯,I˜n,V˜1,0,0ℜ0>1ℜ1≤1<ℜ0 and ℜ2≤1Π2S¯1,⋯,S¯n,I¯1,⋯,I¯n,V¯,Z¯1,0ℜ1>1ℜ1>1 and ℜ2≤ℜ1Π3S^1,⋯,S^n,I^1,⋯,I^n,V^,0,Z^2ℜ2>1ℜ2>1 and ℜ1≤ℜ2


## 6. Numerical Simulations

Numerical simulations are carried out in this section to reinforce the outcome of the Theorems1–4. Our simulation involves considering a special case in which the DENV attacks two types of target cells, i.e. n=2. In this case, model (3) leads to the following:(75)∂S1x,t∂t=d1ΔS1x,t+θ1−β1S1x,tVx,t−ξ1S1x,t,∂S2x,t∂t=d2ΔS2x,t+θ2−β2S2x,tVx,t−ξ2S2x,t,∂I1x,t∂t=d3ΔI1x,t+β1S1x,tVx,t−ϱ1I1x,t,∂I2x,t∂t=d4ΔI2x,t+β2S2x,tVx,t−ϱ2I2x,t,∂Vx,t∂t=d5ΔVx,t+τ1I1x,t+τ2I2x,t−ηVx,t−ψ1H1x,t+ψ2H2x,tVx,t,∂H1x,t∂t=d6ΔH1x,t+μ1Vx,tH1x,t−ϕ1H1x,t,∂H2x,t∂t=d7ΔH2x,t+μ2Vx,tH2x,t−ϕ2H2x,t.The threshold parameters for model (75) are as follows(76)ℜ0=τ1β1θ1ηξ1ϱ1+τ2β2θ2ηξ2ϱ2,ℜ1=μ1τ1β1θ1ηϱ1μ1ξ1+β1ϕ1+μ1τ2β2θ2ηϱ2μ1ξ2+β2ϕ1,ℜ2=μ2τ1β1θ1ηϱ1μ2ξ1+β1ϕ2+μ2τ2β2θ2ηϱ2μ2ξ2+β2ϕ2.We consider the following starting conditions for system (75) based on the work presented in [25, 48, 49]:(77)S1x,0=4001+0.4cos2πx,S2x,0=3001+0.3cos2πx,I1x,0=301+0.5cos2πx,I2x,0=401+0.6cos2πx,Vx,0=41+0.7cos2πx,H1x,0=21+0.7cos2πx,H2x,0=21+0.7cos2πx,x∈0,3.Furthermore, the homogeneous Neumann boundary condition is considered:(78)∂S1∂G⟶=∂S2∂G⟶=∂I1∂G⟶=∂I2∂G⟶=∂V∂G⟶=∂H1∂G⟶=∂H2∂G⟶=0,t>0,x=0,3.To solve system (75) numerically, we choose a domain Λ as Λ=0,3 with a step size 0.02 and step size 0.1 for time. Using MATLAB solver pdepe, we numerically demonstrate the global stability of Π0,Π1,Π2,Π3. During this part, some parameters values are changed while others are assigned fixed estimate parameters (see Table 2). With varied values of parameters β1, β2, μ1, and μ2, we get four scenarios as follows:Case I (Stability of Π0): By choosing β1=0.0001,β2=0.0001, μ1=0.003, and μ2=0.001, we obtain ℜ0=0.5533<1, ℜ1=0.4309<1, and ℜ2=0.3883, and the solutions of system (75) converges to Π0=600,466.6667,0,0,0,0,0. As a result of this, DENV will be cleared (Figure 1). This indicates that Π0 is G.A.S, in harmony with Theorem 1.Case II (Stability of Π1): By choosing β1=0.0004,β2=0.0002, μ1=0.003, and μ2=0.001, we obtain ℜ0=1.7467>1, ℜ1=0.8716483517<1, and ℜ2=0.7067<1. In this case, the affected person has no antibody immune reaction to DENV contamination, and the solutions of system (75) reach to Π1=303.6858,352.1369,19.7542,17.1795,48.7863,0,0 (Figure 2). This indicates that Π1 is G.A.S consistent with Theorem 2.Case III (Stability of Π2): In this situation, we chose β1=0.0004,β2=0.0002, μ1=0.01, and μ2=0.001, and then ℜ0=1.7467>1, ℜ1=1.3261>1, and ℜ2=0.7067. The affected person has only active heterologous antibody immune reaction to DENV infection. The solutions of system (75) reach to Π2=428.5714,411.7647,11.4286,8.2353,20,8.1513,0 (Figure 3). This indicates that Π2 is G.A.S consistent with Theorem 3.Case IV (Stability of Π3): we chose β1=0.0004,β2=0.0002, μ1=0.009, and μ2=0.008. Then, we calculate ℜ0=1.7467>1, ℜ1=1.2926>1, and ℜ2=1.4548>1. In this situation, The solutions of system (75) reach to Π3=480.,430.7692,8.,5.3846,12.50,0,3.7897 (Figure 4), and then the affected person has only active homologous antibody immune reaction to DENV infection. This indicates that Π3 is G.A.S consistent with Theorem 4.Table 2
The values of the parameters of system (75) ([25, 26, 48]).
ParameterValueβ1, β2, μ1, μ2Variedθ112θ214ξ10.02ξ20.03ϱ10.3ϱ20.2τ18τ25η5ψ10.2ψ20.6ϕ10.2ϕ20.1d10.08d20.3d30.04d40.01d50.3d60.3d70.2Figure 1
Solution of system (75) with the initial conditions (77), for Case I. ℜ0≤1, and the steady state Π0=600,466.6667,0,0,0,0,0 is G.A.S. (a, b) Target cells. (c, d) DENV-infected cells. (e) DENV particles. (f) Heterologous antibodies. (g) Homologous antibodies.
(a)(b)(c)(d)(e)(f)(g)Figure 2
Solution of system (75) with the initial conditions (77), for Case II. ℜ1≤1<ℜ0 and ℜ2≤1, and the steady state Π1=303.6858,352.1369,19.7542,17.1795,48.7863,0,0 is G.A.S. (a, b) Target cells. (c, d) DENV-infected cells. (e) DENV particles. (f) Heterologous antibodies. (g) Homologous antibodies.
(a)(b)(c)(d)(e)(f)(g)Figure 3
Solution of system (75) with the initial conditions (77), for Case III. ℜ1>1 and ℜ2≤ℜ1, and the steady state Π2=428.5714,411.7647,11.4286,8.2353,20,8.1513,0 is G.A.S. (a, b) Target cells. (c, d) DENV-infected cells. (e) DENV particles. (f) Heterologous antibodies. (g) Homologous antibodies.
(a)(b)(c)(d)(e)(f)(g)Figure 4
Solution of system (75) with the initial conditions (77), for Case IV. ℜ2>1 and ℜ1≤ℜ2, and the steady state Π3=480.,430.7692,8.,5.3846,12.50,0,3.7897 is G.A.S. (a, b) Target cells. (c, d) DENV-infected cells. (e) DENV particles. (f) Heterologous antibodies. (g) Homologous antibodies.
(a)(b)(c)(d)(e)(f)(g)In the following, we discuss the Local stability of system (75) at its steady states (see [38, 50, 51]).Let0=ω˜1<ω˜2<… be the complete set of eigenvalues of Lablacian operator −Δ on Λ with homogeneous Neumann boundary conditions, and Eω˜m be the eigenspace corresponding to ω˜m, where m=1,2,…. Let Φmj:j=1,2,…,dimEω˜m be an orthonormal basis of Eω˜m, and(79)Amj=kΦmj:k∈ℝ7,A=S1,S2,I1,I2,V,H1,H2∈C1Λ¯7:∂S1∂G⟶=∂S2∂G⟶=∂I1∂G⟶=∂I2∂G⟶=∂V∂G⟶=∂H1∂G⟶=∂H2∂G⟶=0 on∂Λ.Then, we have(80)Am=⊕j=1dimEω˜mAmj,A=⊕m=1∞Am.The linearization of model (75) at steady state Πll=0,1,2,3 is given by(81)∂ℳ∂t=D˜Δℳ+J˜Πlℳ,l=0,1,2,3,where ℳ=S1,S2,I1,I2,V,H1,H2T, D˜=diagd1,d2,d3,d4,d5,d6,d7 is the diagonal matrix made of diffusion constants, and J˜Πl is the Jacobian matrix of the reaction, which is given by(82)J˜=J˜Πl=−β1V−ξ1000−β1S1000−β2V−ξ200−β2S200β1V0−ϱ10β1S1000β2V0−ϱ2β2S20000τ1τ2−η−ψ1H1−ψ2H2−ψ1V−ψ2V0000μ1H1μ1V−ϕ100000μ2H20μ2V−ϕ2.LetLℳ=D˜Δℳ+J˜Πlℳ, then Am is invariant under the operator L for m≥1. Furthermore, Γ is an eigenvalue of L if and only if it is an eigenvalue of the matrix A˜=J˜−D˜ω˜m, for some m≥1. Consequently, there is an eigenvector in Am.To calculating the eigenvalues of the matrixA˜=J˜−D˜ω˜, we will obtain Γi which are the roots of the characteristic equation(83)detΓI−J˜+D˜ω˜i=0,where I is the identity matrix. All eigenvalues Γi of the matrix A˜ must satisfy the condition ReΓi<0,i=1,2,⋯,7, to be locally stable. Using the parameters as given in Table 2, and for all Cases I–IV, we determine the eigenvalues Γi, i=1,2,…,7 as follows: In case I; the characteristic equation of the matrix A˜ at Π0=600,466.6667,0,0,0,0,0 takes the form(84)Γ3+a1Γ2+b1Γ+c10.02+0.08ω˜+Γ0.03+0.3ω˜+Γ0.2+0.3ω˜+Γ0.1+0.2ω˜+Γ=0,where(85)a1=5.5+0.35ω˜,b1=0.411ω˜+1.8467+0.0154ω˜2,c1=0.134+0.05887ω˜+0.0053ω˜2+0.00012ω˜3.Clearly, (84) has four negative real roots: −0.02+0.08ω˜, −0.03+0.3ω˜, −0.2+0.3ω˜, −0.1+0.2ω˜. The other roots of (84) are determined by the equation(86)Γ3+a1Γ2+b1Γ+c1=0.Let(87)f1Γ=Γ3+a1Γ2+b1Γ+c1.Clearlyf10=c1>0, f1′0=b1>0, and limΓ⟶+∞f1Γ=+∞, which implies that the (84) has not a positive real root, therefore Π0 is stable.In case II; the characteristic equation ofA˜ at Π0 is given by the following form(88)Γ3+a2Γ2+b2Γ+c20.02+0.08ω˜+Γ0.03+0.3ω˜+Γ0.2+0.3ω˜+Γ0.1+0.2ω˜+Γ=0,where(89)a2=0.35ω˜Γ2+5.5,b2=0.0154ω˜2+0.411ω˜+0.1733,c2=0.00012ω˜3+0.0053ω˜2+0.0351ω˜−0.224.Define(90)f2Γ=Γ3+a2Γ2+b2Γ+c2.As a resultf20=c2<0 if 0.00012ω˜3+0.0053ω˜2+0.0351ω˜<0.224, f2′0=b2>0, andlimΓ⟶+∞f2Γ=+∞. Therefore, for some ω˜, there is a characteristic root of (88) with positive real part, hence, Π0 is unstable.In a similar manner as above, we present in the Table3 steady states Π0,Π1,Π2 and Π3, along with the real part of all corresponding eigenvalues Γi. We see that, whenever a steady state is G.A.S, all other steady states are unstable.Table 3
Local stability of nonnegative steady statesΠi,i=0,1,2,3 in Cases I–IV of model (75) with ω˜=0.25.
CaseSteady stateRe.Γi,i=1,2,…7StabilityIΠ0=600.,466.6667,0,0,0,0,0−0.1110,−0.2572,−5.2193, −0.04,−0.105,−0.2750,−0.15StableIIΠ0=600.,466.6667,0,0,0,0,00.1706,−0.2277,−5.5304−0.04,−0.105,−0.275,−0.15UnstableΠ1=303.6858,352.1369,19.7543,17.1795,48.7863,0,0−0.1286,−0.1012,−0.0384−0.1081,−0.2391,−5.3377,−0.0384StableIIIΠ0=600.,466.6667,0,0,0,0,00.1706,−0.2277,−5.5304,−0.04−0.105,−0.2750,−0.15UnstableΠ1=303.6858,352.1369,19.7543,17.1795,48.7863,0,0−5.3377,−0.2391,−0.0384−0.0384,−0.1081,0.2129,−0.1012UnstableΠ2=428.5714,411.7647,11.4286,8.2353,20,8.1513,0−6.9265,−0.2326,−0.0657−0.0657,−0.0514,−0.1079,−0.13StableIVΠ0=600.,466.6667,0,0,0,0,0−0.04,−0.105,−5.5303,0.1706−0.2277,−0.275,−0.15UnstableΠ1=303.6858,352.1369,19.7543,17.1795,48.7863,0,0−5.3377,−0.2391,−0.0384−0.0384,−0.1081,0.1641,0.2402UnstableΠ2=415.3846,406.4516,12.3077,9.0323,22.2222,7.3151,0−6.7625,−0.2332,−0.0636−0.0636,−0.0527,−0.1081,0.0278UnstableΠ3=480.,430.7692,8.,5.3846,12.50,0,3.7897−7.5882,−0.2327,−0.04551−0.0455,−0.1063,−0.0456,−0.1625Stable


## 7. Discussion

Infection by dengue is the most significant arboviral disease, causing severe illness and representing an increasing public health threat. In experimental research, mathematical modeling is a critical tool for understanding disease mechanisms. In this study, a dynamical model was used to describe the behavior of secondary DENV infection. The DENV particles’ and cells’ spatial mobility were incorporated into the model. We took into account that several types of target cells are susceptible to DENV infection. By using2n+3 PDEs, Our model describes how target cells, infected cells, DENV particles, heterologous and homologous antibodies interact with each other. In order to demonstrate that the model is biologically relevant, we first established that the key variables are nonnegative and bounded.By using Lyapunov functions and using Lyapunov-LaSalle asymptotic stability theorems, we studied the stability analysis of PDEs system, and we proved that:(a)
If the basic reproduction numberℜ0≤1, the uninfected steady state Π0 is globally asymptotically stable (G.A.S) for all diffusion coefficients dj, it means that the virus is cleared and the infection dies out.(b)
Whenℜ0>1, other steady states exist; namely, Π1, the infected steady state with inactive immune antibody response; Π2, infected steady state with only active heterologous antibodies; Π3, infected steady state with only active homologous antibodies, such as:(i)
If the heterologous antibody immune response activation numberℜ1≤1, and the homologous antibody immune response activation number ℜ2≤1, then Π1 is G.A.S for all dj. The DENV virus persists in the host and the infection becomes chronic and there is no antibody immune reaction to DENV particles. In this case, despite having DENV infection, the patient’s immune system is ineffective.(ii)
Ifℜ1>1 and ℜ2≤ℜ1, then Π2 is G.A.S for all dj. As a result, the DENV infection becomes chronic, and the only immune response to DENV infection is an active heterologous antibody response, and the patient, in this case, suffers from DENV infection with effective heterologous antibody immune response.(iii)
Ifℜ2>1 and ℜ1≤ℜ2, then Π3 is G.A.S for all dj. Thus, DENV infection becomes chronic, with only an active homologous antibody response acting as an immune response. In this case, the patient suffers from DENV infection with an active homologous antibody response.In light of the above theoretical analysis, we conclude that the global dynamics of the model are fully determined by the threshold parametersℜ0,ℜ1,ℜ2. The results of Theorems 1–4 were supported using numerical simulations.As an extension of our model, we can consider delays during the process of infection and production of DENV; thus, DENV dynamics model can be adapted to incorporate delays into the infection equation and/or the DENV production equation.



---

*Source: 2106910-2022-11-21.xml*