# Effect of General Cross-Immunity Protection and Antibody-Dependent Enhancement in Dengue Dynamics

Authors: Vanessa Steindorf; Sergio Oliva; Jianhong Wu; Maíra Aguiar
Journal: Computational and Mathematical Methods
 (2022)

---

## Abstract

A mathematical model to describe the dynamic of a multiserotype infectious disease at the population level is studied. Applied to dengue fever epidemiology, we analyse a mathematical model with time delay to describe the cross-immunity protection period, including a key parameter for the antibody-dependent enhancement (ADE) effect, the well-known features of dengue fever infection. Numerical experiments are performed to show the stability of the coexistence equilibrium, which is completely determined by the basic reproduction number and by the invasion reproduction number, as well as the bifurcation structures for different scenarios of dengue fever transmission in a population. The model shows a rich dynamical behavior, from fixed points and periodic oscillations up to chaotic behaviour with complex attractors.


---

## Body

## 1. Introduction

Dengue fever (DF) is a mosquito-borne disease, a huge public health problem affecting specially the tropical and subtropical countries. With approximately 3 billion people at risk of acquiring the infection, it is estimated that 390 million dengue infections occur every year, of which 96 million manifest symptoms with any level of disease severity [1].DF is caused by four antigenically related but distinct serotypes (DENV-1 to DENV-4) [2, 3]. While a primary natural dengue infection is often asymptomatic or mild, the clinical response on exposure to a second serotype is complex and may depend on factors such as patient age, dengue serotype, sequence of infection, and the interval between infection by one serotype and exposure to a second serotype. Infection by one serotype confers life-long immunity to that serotype and a period of temporary cross-immunity (TCI) protection to other serotypes. However, individuals undergoing a secondary dengue infection with a heterologous serotype have a higher risk of developing the severe form of the disease. There is good evidence that sequential infection increases the risk of developing severe disease, due to a process described as antibody-dependent enhancement (ADE), where the preexisting antibodies to previous dengue infection do not neutralize but rather enhance the new infection [3–8]. Treatment of uncomplicated dengue cases is only supportive, and severe dengue cases require proactive treatment of hemorrhagic symptoms.Up to date, two tetravalent dengue vaccines have completed phase 3 clinical trial: Dengvaxia (developed by Sanofi Pasteur) that is now licensed in several countries [9, 10] and the DENVax vaccine, developed by Takeda Pharmaceutical Company [11, 12]. While Dengvaxia has resulted in serious adverse events in seronegative individuals compared with age-matched seronegative controls [13–17], long-term surveillance consisting of prudent and careful observation of DENVax vaccine recipients is required [18, 19], since negative vaccine efficacy was estimated for vaccinated seronegative individuals who were infected with serotype 3, as opposed to an intermediate efficacy for seropositive [20, 21].Dengue fever epidemiological dynamics shows large fluctuations in disease incidence, and several mathematical models describing the transmission of dengue viruses have been proposed to explain the irregular behavior of dengue epidemics. A careful review of deterministic dengue modeling was published in 2011 [22], where two main modeling approaches were considered, the vector-host and the host-to-host transmission. In the first approach, the fluctuations in the mosquito dynamics and climate change are included and assumed to affect the disease transmission [23–27], while in the second approach, the effect of seasonality (mimicking the vectorial dynamics) appears to be essential to explain the intra-annual fluctuations in disease cases [28–31].Multistrain dengue dynamics are generally modelled with extended SIR- (susceptible-infected-recovered-) type models, and the combination of biological aspects such as TCI and ADE effects has been studied by several authors involving four strains [23, 24, 32, 33], but always limiting the effect of ADE to increase the contribution of secondary cases to the force of infection. Aguiar et al. [28, 29, 34–36] have investigated a minimalistic two-infection dengue model, an extension of a model initially suggested and preliminary analysed in Ferguson et al. [23], where deterministic chaos was found in wider parameter regions, not needing to restrict the infectivity on secondary infection to one or another region in parameter space. Consideration of temporary cross-immunity protection together with ADE gave bifurcations up to chaotic attractors in much wider and also unexpected parameter regions, not predicted by previous models [37, 38].Theoretical mathematical analysis of this model was made by Aguiar et al. [28, 35] and by Kooi et al. [38]. A detailed time series simulations were investigated, and various bifurcation phenomena were observed. Considering symmetry of strains, the authors have shown the analytic formulations for the equilibrium, as well as the analysis of the bifurcation structures observed in the system. Furthermore, using numerical simulations to obtain the bifurcation points, the analytic steady states and the stability for the disease-free equilibrium and for the boundary equilibrium were shown. Aguiar et al. [29] have also shown that the combination of TCI and ADE is the most important feature to drive the complex dynamics in the system, more than the detailed number of dengue serotypes to be added in the model. These models, however, are formulated in terms of Ordinary Differential Equations (ODE), without incorporating time lags in any of the state class for the human population, neither for recovery nor for the incubation period, for example. Consideration of time delay as a mechanism to model the dynamic of dengue fever would be eventually useful to better describe the disease immunological response where a short cross-protection period to related serotypes is observed after the individual recovers from a primary natural infection.Simple delay models are frequently used to describe biological models and infectious disease dynamics. Delay differential equations are useful to describe, for example, the population or disease behaviour which intrinsically depends on a certain period of past time. Models with general time delay (relapse period) for recovery individuals were used by Van den Driessche in [39]. A mathematical model describing a general multistrain disease has also been recently described by Chen et al. in [40]. The authors proposed a delay diffusive two-strain disease model considering a SIR structure, constant recruitment rate, and constant time delay representing the length of the immunity period. A model for dengue fever that describes a time delay between infectious and infected hosts was described by Sakdanupahp and Moore in [41]. This model is relatively simple, with four equations, but brings two different constant time delays, one for the infectious host to the infected host and another one for the vector population. However, this study did not consider different serotypes. Steady states were obtained, stability was proven, and numerical analysis and comparison with empirical data were performed.The model described by Cai et al. in [25] preserves the SIR structure, considers multistrain structure, and includes incubation period time and vector population dynamic. The authors focused on the evolution of pathogen, and according to the authors, the model could be used to describe the dynamics of dengue fever or malaria. Guan et al. [42] also proposed a dengue fever model with time delay. The time delay included in the model refers to a time incubation of the virus in an infected subpopulation and in the infected mosquito population. In this work, numerical simulations and analytical results were obtained assuming four cases for the constant time delay.Using a simpler SIR model, Huang et al. [43] considered an infinite distributed delay on complex population network in order to give some insights on biological and social networks. Numerical experiments confirmed that, in this case, with the delay slowing down the extinction of the disease when the threshold value is smaller than 1, while the delay accelerates the spreading of the disease if the threshold value is bigger than 1. Distributed delay was also used by Xu et al. [44] in a SVEIR type model, including vaccination and general incidence function. Results showed that the distributed delay has no impact on the qualitative behaviour and global dynamics of the system.Infinite distributed delay was used by Steindorf et al. [45] to describe a general time delay representing the cross-protection period on dengue modelling. Allowing a general period of immunity, the time delay was used to model the phenomenon of temporary cross-immunity protection, assuming that the individual recovered from a primary infection may not be immediately susceptible to a heterologous dengue serotype infection. The model includes a key parameter ϕ representing ADE, with results showing that the disease fades out if the basic reproduction number is smaller than 1. Otherwise, depending on ϕ values, different scenarios are observed, including coexistence of serotypes, prevalence of one serotype, periodic outbreaks, and also complex dynamics for a certain parameter regions where the disease transmission is difficult to be predicted.In this work, we analyse a mathematical model to describe the propagation of multistrain infectious diseases. We propose a system of integro-differential equations (IDE), motivated by dengue fever epidemiology and its well-studied antibody-dependent enhancement (ADE) phenomenology and temporary cross-immunity period. Aimed at understanding the effect of the general time delay on the model, which describes the length of the temporary cross-immunity protection period and assuming a constant parameterϕ representing the ADE effect, numerical experiments are performed, showing the bifurcation structures for different scenarios of dengue fever transmission in a population.


## 2. Mathematical Model

In this section, we present the modeling framework developed to describe dengue multistrain dynamics, considering the temporary cross-immunity period and antibody-dependent enhancement effects. For the sake of simplicity, we evaluate a two-serotype host-to-host dynamical model.LetNt be the total host population size in a region, divided in classes according to the disease-related status. Individuals are assumed to be susceptible to all serotypes (St), infected by one serotype i (Iit), temporarily immune to all serotypes after infection by serotype i (Cit), and recovered from serotype i, and after the temporary cross-immunity protection period, they become susceptible to the other serotype, Rit. Secondary infection can occur with a heterologous serotype; that is, Iijt represent infected individuals by serotype j after being recovered from serotype i, with i,j=1,2 and i≠j. Finally, recovered individuals for all serotypes Rt at time t. Let d be the natural mortality.For the sake of simplicity, the mosquito population is not considered explicitly in the model. That is supported by the findings published in [31], showing that the addition of seasonality into these models can explain well the vector dynamics, mimicking its abundance and decrease overtime, without the need of adding many more equations into the system. In this paper, we focus our analysis into the multistrain aspect of dengue epidemiology. We analyze first the non-seasonal model with already shows a very rich dynamical behaviour. Ongoing refinements consider the use of seasonal forcing to describe the vector behaviour in endemic countries. Explicit vector dynamics may be needed to evaluate disease control measures such as vector control and vaccination [46–48].Letβi be the transmission rate of a primary infection with serotype i and αi be the transmission rate of a secondary infection with serotype i. Individuals in the infectious classes Iit remain in this class with average time 1/γ, assuming that the length in this class is exponentially distributed. After the infectious period, individual recovers and remain temporarily immune to all serotypes in the class Cit, due to the cross-immunity protection period, with Pit being the fraction of individuals recovered from an infection with serotype i at time t, remaining cross-immune protected to all serotypes for a certain amount of time. Here, we assume that Pi0=1 and Pi∞=0 and Pit is not increasing. Moreover,
(1)∫0∞Pisds=1ωi<∞.Thus, the number of cross-protected individuals is given by(2)Cit=C0t+∫0tγIsPit−se−dt−sds,where C0 is the number of protected individuals at time t=0. Of course, C0 must satisfy limt⟶∞C0t=0.After the cross-protection period, individual becomes susceptible again, now to a different serotype, but life-long protected to the serotype causing the first infection. Thus, an individual inRi class can be infected with a rate αj with a heterologous dengue serotype. After being infected by the two strains, individual became recovered from all strains, remaining lifelong immune.The ADE phenomenon explains the enhanced risk of disease severity in secondary infections caused by a heterologous dengue serotype, due to increasing viral replication mediated by preexisting antibodies that were produced during a primary infection [2, 7]. According to Katzelnick et al. [7], ADE occurs at a specific range of antibody concentrations. Furthermore, Rothman et al. [8] affirms that depending on the specific antibody concentration, dengue virus antibodies can inhibit viral infection by neutralization or enhance the infection, with increased viral replication, by facilitating the entry of the complex antibody-heterologous virus into target cells. Hence, the ADE effect enhances infection during a secondary infection caused by a heterologous serotype due to high concentrations of preexisting antibodies on individuals who already recovered from a primary infection.Here, the epidemiological effect of ADE is described by a constant parameterϕ. We assume that this increasing (or decreasing) of viral load during the secondary infection is not related to the increased transmissibility but, rather to the capacity of the immune system to respond to the secondary infection, protecting or enhancing the disease. Thus, in this study, the ADE effect is parametrized with a factor ϕ, increasing (ϕ>1) or decreasing (ϕ<1) the probability of a recovered individual to acquire a secondary infection with a heterologous serotype, assuming that the individual recovered from a primary infection carries antibody load levels to either neutralize or to enhance the infection.Based on the assumption described above and in the studies by Van den Driessche in [49, 50] and taking the derivative of equation (2), assuming C0=0, the mathematical model can be written as follows:
(3)dStdt=dNt−dSt−β1StNtI1t−β2StNtI2t−β2StNtI12t−β1StNtI21t,dI1tdt=−dI1t+β1StNtI1t+β1StNtI21t−γI1t,dI2tdt=−dI2t+β2StNtI2t+β2StNtI12t−γI2t,dC1tdt=γI1t−dC1t+∫0tγI1sP1′t−se−dt−sds,dC2tdt=γI2t−dC2t+∫0tγI2sP2′t−se−dt−sds,dR1tdt=−dR1t−α2ϕR1tNtI12t−α2ϕR1tNtI2t−∫0tγI1sP1′t−se−dt−sds,dR2tdt=−dR2t−α1ϕR2tNtI21t−α1ϕR2tNtI1t−∫0tγI2sP2′t−se−dt−sds,dI12tdt=−dI12t−γI12t+α2ϕR1tNtI2t+α2ϕR1tNtI12t,dI21tdt=−dI21t−γI21t+α1ϕR2tNtI1t+α1ϕR2tNtI21t,dRtdt=−dRt+γI12t+γI21t.


## 3. Results

In Steindorf et al. [45], the model is analysed following Miller’s [51] results that the limiting system ensures that the system proposed has an equilibrium and the equilibrium of the system coincides with those of their limiting equation, followed by defining an appropriated Banach space in order to have a well-posed system and to have solutions well defined.The qualitative analysis of this model is available in Steindorf et al. [45], where the analytical four equilibria of the limiting system in the invariant region are presented. The local stability analysis for the disease-free equilibria and for the boundary endemic equilibrium was proven. Using the theory proposed by Brauer [52], the stability of the solutions of the system was proven and was completely determined by the basic reproduction number and by the invasion reproduction number defined mathematically as a threshold value for stability. In this work, Steindorf et al. [45] used Lyapunov functions to prove the global dynamics.It was proven that if the basic reproduction numberR0=maxR1,R2=maxβ1/d+γ,β2/d+γ≤1, the disease-free equilibrium is globally asymptotically stable. If R0>1 and the invasion reproduction number RInv=R1/R2+R2−1α1ϕγM2/β2d+γ≤1, where M2=∫0∞P′2se−dsds, the boundary equilibrium related to the the highest infection rate, βi, is globally asymptotically stable, while the other boundary equilibrium is always unstable. Finally, if R0>1 and RInv>1, the coexistence endemic equilibrium exists and it is positive in the invariant region.In order to show the stability of the coexistence equilibrium, numerical experiments were performed. Although the existence of a coexistence endemic equilibrium could be shown and proven analytically, the theoretical analysis of the stability using the linearized theory, for this case, was not successful, since we will have to deal with a transcendental characteristic equation. Thus, analysis of the local stability of this equilibrium is studied through numerical approach.The main results are summarized in Table1.Table 1
Stability of the system: for asymmetric parameters.
Basic numberInvasion numberDisease-freeBoundaryD1BoundaryD2CoexistenceD3R1<R2RInvR0≤1RInv<1GASNoNoNoR0>1RInv<1UnstableNo (ifR1<1)GASNoUnstableUnstable (ifR1>1)GASNoR0>1RInv>1UnstableUnstableUnstableLAS forϕ<ϕcUnstableUnstableUnstableHB (stable limit cycle)UnstableUnstableUnstableUnstableϕ>ϕc


## 4. Numerical Results

The numerical values for the parameters are shown in the Table2.Table 2
Numerical values of the parameters used for the simulations.
ParameterMeaningValueReferencedMortality rate0.015y−1[53]γRecovery rate52y−1[2, 54]βiInfection rate (individuals susceptible)40−200⋆αiReinfection rate (individuals recovered fromj)40−200⋆ϕADE factor0−5[7]⋆These values were calculated to give reasonable basic reproduction numbers for dengue. For instance, using data for Brazil, Massad et al. [55, 56] estimated the range of R0 to be between 1.38 and 7.86. Reiner et al. [57] estimated the range of the basic reproduction number to vary from 0.76 to over 5, using data from Peru, while Fergusson et al. [58] estimated the R0 from 1.38 to 7.70, with average 3.2, using data from Thailand.We consider a function to represent the cross-immunity period, assuming that after one year, the individual is no longer protected to a new infection [2]. For this initial numerical analysis, we choose to work with the cubic polynomial (in order not to work with infinity time),
(4)Ps=2s3−3s2+1,0≤s≤1,0,s≥1,with average time being 1/2 year [2].


### 4.1. Stability of the Coexistence Endemic Equilibrium (CEE)

The stability of the CEE will be studied using numerical experiments. The CEE,D3, will be in the positive region only if the parameters satisfy the condition RInv>1. If it satisfied, then also the boundary equilibrium (BE) will lose the stability. In order to study the stability of the CEE, we are going to analyse numerically the roots of the characteristic equation. To do that, we choose some values of the parameter ϕ, using it as a bifurcation parameter.We evaluate three different case scenarios:(i)
Symmetric CaseR0=R1>1In this case, the infection rates are the same for both strains. For the simulations, we fixedβi=180, giving the value for R0=3.46. For the symmetric case scenarios, the coexistence equilibrium will be always positive in the invariant region. This fact was also proven analytically in [45]. The equilibrium will be stable for values of ϕ<0.032; that is, when the Hopf bifurcation occurs, this is, for ϕc=0.032, where a purely imaginary part of eigenvalues appears (Figure 1(a)).
(ii)
Asymmetric Case forR0>1, R1<1Figure 1
Symmetric case (R0=R1>1). In a vicinity of ϕc=0.032, the maximum and minimum values for the susceptible population are plotted.In this case scenario, with values for the parameters in Table2, we fixed β1=45 and β2=180, which gives R1=0.87 and R0=3.46. To study the stability of the CEE, it is necessary that the condition RInv>1 it satisfied. This occurs for ϕ>1.23. Then, for all values of ϕ>1.23, we have the existence of the CEE and the correspondent eigenvalues, as shown in Figures 2(a)–2(f) .Figure 2
Asymmetric case (R0>1, R1<1). Figures show the eigenvalues of the coexistence endemic equilibrium (CEE) in the complex plane, for each value of ϕ. (d) A purely imaginary eigenvalue appears for ϕ≈2.19.
(a)
ϕ=1.24(b)
ϕ=1.4(c)
ϕ=2(d)
ϕ=2.19(e)
ϕ=2.4(f)
ϕ=2.52Figures2(a)–2(f) show that the characteristic equation has a pair of conjugated complex roots that change the sign of the real part as the parameter ϕ increases. As ϕ increases from small values through critical value, ϕc, the steady state changes from a stable focus to an unstable steady state. Therefore, the Hopf bifurcation occurs when the parameter ϕ≈2.19. Thus, closed periodic orbit will be found in a small neighbourhood of ϕc. This can also be seen in the bifurcation diagram shown in Figure 3(a).
(iii)
Asymmetric Case forR0>1, R1>1Figure 3
Asymmetric case (R0>1, R1<1). In the vicinity of ϕc=2.19, the maximum and minimum values for the susceptible population are plotted.The values for the parameters used in the simulations for this case scenario can be found in Table2. We fixed β1=120 and β2=180, which gives R1=2.31 and R0=3.46, respectively, and thus RInv>1 for ϕ>0.205. Therefore, for all values of ϕ>0.205, we have the CEE at the positive region, and we plotted the correspondents’ eigenvalues, as shown in Figures 4(a)–4(f).Figure 4
Asymmetric case (R0>1, R1>1). Figures show the eigenvalues of the coexistence endemic equilibrium (CEE) in the complex plane, for each value of ϕ. (c) A purely imaginary eigenvalue appears for ϕ≈0.244.
(a)
ϕ=0.209(b)
ϕ=0.235(c)
ϕ=0.244(d)
ϕ=0.8(e)
ϕ=1(f)
ϕ=2Figures4(a)–4(f) show that the characteristic equation has a pair of conjugated complex roots that change the sign of the real part of the eigenvalues as the parameter ϕ increases. Thus, the Hopf bifurcation occurs when the parameter ϕ≈0.244. This can also be seen in Figure 5(a), and for this case scenario, also the supercritical Hopf bifurcation occurs. Furthermore, after the Hopf bifurcation, the CEE will be unstable.Figure 5
Asymmetric case (R0>1, R1>1). Parameter ϕ varies in a vicinity of ϕc=0.244, while the maximum and minimum values for the susceptible population are plotted.


### 4.2. Bifurcation Structure

In the previous section, we have shown, using numerical experiments, that for a range of the parameterϕ, the equilibrium is stable for the three case scenarios. A Hopf bifurcation occurs for a critical value ϕc, followed by periodic oscillations. We evaluate bifurcation diagrams to study the dynamics for the different ϕ values being farther than ϕc.
(i)
Symmetric CaseR0=R1>1Figure6(a) shows a bifurcation diagram where the maximum and minimum values for the infected population are plotted, while the value for the parameter ϕ is varying. After the Hopf bifurcation, a chaotic behavior is observed. For bigger values, a periodic behavior is found. Figures 7(a)–7(f) illustrate the space phase for different values of the parameter ϕ, with complicated attractor and limit cycles identified.
(ii)
Asymmetric Case forR0>1, R1<1Figure 6
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) are plotted.Figure 7
Symmetric case (R0=R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.01(b)
ϕ=0.032(c)
ϕ=0.1(d)
ϕ=0.5(e)
ϕ=1.5(f)
ϕ=2.5The Hopf bifurcation occurs atϕc=2.19 (see Figure 3(a)). The solutions exhibit a small amplitude limit cycle around the endemic equilibrium. A stable limit cycle arises close to the critical bifurcation point and goes away from the unstable equilibrium. Thus, it is possible to conclude that a supercritical Hopf bifurcation has occurred. After the Hopf bifurcation, the diagram shows chaotic behavior for bigger ϕ values (see Figure 8(a)). Figures 9(a)–9(f) illustrate the space phase for different values of the parameter ϕ, with a limit cycle occuring at ϕ=2.19. For bigger values of ϕ, chaotic dynamics is observed.
(iii)
Asymmetric Case forR0>1, R1>1Figure 8
Asymmetric case (R0>1,R1<1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.Figure 9
Asymmetric case (R0>1,R1<1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=2.1(b)
ϕ=2.19(c)
ϕ=2.21(d)
ϕ=2.25(e)
ϕ=2.5(f)
ϕ=3.5After the Hopf bifurcation, the diagram shows chaotic behavior and complex dynamics for bigger values of theϕ parameter (see Figure 10(a)). Figures 11(a)–11(f) illustrate the space phase for different values of the parameter ϕ, showing chaotic dynamics and complicated attractors for different values of ϕ.Figure 10
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.Figure 11
Asymmetric case (R0>1, R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum for the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.22(b)
ϕ=0.261(c)
ϕ=0.55(d)
ϕ=0.8(e)
ϕ=0.99(f)
ϕ=2By varying the ADE parameter, the system shows a rich dynamical behavior for all the three case scenarios. For small values ofϕ, we have a fixed point; hence, the coexistence endemic equilibrium is stable. The Hopf bifurcation occurs for a specific range of ϕ values, from where disease propagation is possible and periodic solutions are observed. Complicated attractors and chaotic behavior are also observed.


## 4.1. Stability of the Coexistence Endemic Equilibrium (CEE)

The stability of the CEE will be studied using numerical experiments. The CEE,D3, will be in the positive region only if the parameters satisfy the condition RInv>1. If it satisfied, then also the boundary equilibrium (BE) will lose the stability. In order to study the stability of the CEE, we are going to analyse numerically the roots of the characteristic equation. To do that, we choose some values of the parameter ϕ, using it as a bifurcation parameter.We evaluate three different case scenarios:(i)
Symmetric CaseR0=R1>1In this case, the infection rates are the same for both strains. For the simulations, we fixedβi=180, giving the value for R0=3.46. For the symmetric case scenarios, the coexistence equilibrium will be always positive in the invariant region. This fact was also proven analytically in [45]. The equilibrium will be stable for values of ϕ<0.032; that is, when the Hopf bifurcation occurs, this is, for ϕc=0.032, where a purely imaginary part of eigenvalues appears (Figure 1(a)).
(ii)
Asymmetric Case forR0>1, R1<1Figure 1
Symmetric case (R0=R1>1). In a vicinity of ϕc=0.032, the maximum and minimum values for the susceptible population are plotted.In this case scenario, with values for the parameters in Table2, we fixed β1=45 and β2=180, which gives R1=0.87 and R0=3.46. To study the stability of the CEE, it is necessary that the condition RInv>1 it satisfied. This occurs for ϕ>1.23. Then, for all values of ϕ>1.23, we have the existence of the CEE and the correspondent eigenvalues, as shown in Figures 2(a)–2(f) .Figure 2
Asymmetric case (R0>1, R1<1). Figures show the eigenvalues of the coexistence endemic equilibrium (CEE) in the complex plane, for each value of ϕ. (d) A purely imaginary eigenvalue appears for ϕ≈2.19.
(a)
ϕ=1.24(b)
ϕ=1.4(c)
ϕ=2(d)
ϕ=2.19(e)
ϕ=2.4(f)
ϕ=2.52Figures2(a)–2(f) show that the characteristic equation has a pair of conjugated complex roots that change the sign of the real part as the parameter ϕ increases. As ϕ increases from small values through critical value, ϕc, the steady state changes from a stable focus to an unstable steady state. Therefore, the Hopf bifurcation occurs when the parameter ϕ≈2.19. Thus, closed periodic orbit will be found in a small neighbourhood of ϕc. This can also be seen in the bifurcation diagram shown in Figure 3(a).
(iii)
Asymmetric Case forR0>1, R1>1Figure 3
Asymmetric case (R0>1, R1<1). In the vicinity of ϕc=2.19, the maximum and minimum values for the susceptible population are plotted.The values for the parameters used in the simulations for this case scenario can be found in Table2. We fixed β1=120 and β2=180, which gives R1=2.31 and R0=3.46, respectively, and thus RInv>1 for ϕ>0.205. Therefore, for all values of ϕ>0.205, we have the CEE at the positive region, and we plotted the correspondents’ eigenvalues, as shown in Figures 4(a)–4(f).Figure 4
Asymmetric case (R0>1, R1>1). Figures show the eigenvalues of the coexistence endemic equilibrium (CEE) in the complex plane, for each value of ϕ. (c) A purely imaginary eigenvalue appears for ϕ≈0.244.
(a)
ϕ=0.209(b)
ϕ=0.235(c)
ϕ=0.244(d)
ϕ=0.8(e)
ϕ=1(f)
ϕ=2Figures4(a)–4(f) show that the characteristic equation has a pair of conjugated complex roots that change the sign of the real part of the eigenvalues as the parameter ϕ increases. Thus, the Hopf bifurcation occurs when the parameter ϕ≈0.244. This can also be seen in Figure 5(a), and for this case scenario, also the supercritical Hopf bifurcation occurs. Furthermore, after the Hopf bifurcation, the CEE will be unstable.Figure 5
Asymmetric case (R0>1, R1>1). Parameter ϕ varies in a vicinity of ϕc=0.244, while the maximum and minimum values for the susceptible population are plotted.


## 4.2. Bifurcation Structure

In the previous section, we have shown, using numerical experiments, that for a range of the parameterϕ, the equilibrium is stable for the three case scenarios. A Hopf bifurcation occurs for a critical value ϕc, followed by periodic oscillations. We evaluate bifurcation diagrams to study the dynamics for the different ϕ values being farther than ϕc.
(i)
Symmetric CaseR0=R1>1Figure6(a) shows a bifurcation diagram where the maximum and minimum values for the infected population are plotted, while the value for the parameter ϕ is varying. After the Hopf bifurcation, a chaotic behavior is observed. For bigger values, a periodic behavior is found. Figures 7(a)–7(f) illustrate the space phase for different values of the parameter ϕ, with complicated attractor and limit cycles identified.
(ii)
Asymmetric Case forR0>1, R1<1Figure 6
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) are plotted.Figure 7
Symmetric case (R0=R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.01(b)
ϕ=0.032(c)
ϕ=0.1(d)
ϕ=0.5(e)
ϕ=1.5(f)
ϕ=2.5The Hopf bifurcation occurs atϕc=2.19 (see Figure 3(a)). The solutions exhibit a small amplitude limit cycle around the endemic equilibrium. A stable limit cycle arises close to the critical bifurcation point and goes away from the unstable equilibrium. Thus, it is possible to conclude that a supercritical Hopf bifurcation has occurred. After the Hopf bifurcation, the diagram shows chaotic behavior for bigger ϕ values (see Figure 8(a)). Figures 9(a)–9(f) illustrate the space phase for different values of the parameter ϕ, with a limit cycle occuring at ϕ=2.19. For bigger values of ϕ, chaotic dynamics is observed.
(iii)
Asymmetric Case forR0>1, R1>1Figure 8
Asymmetric case (R0>1,R1<1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.Figure 9
Asymmetric case (R0>1,R1<1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=2.1(b)
ϕ=2.19(c)
ϕ=2.21(d)
ϕ=2.25(e)
ϕ=2.5(f)
ϕ=3.5After the Hopf bifurcation, the diagram shows chaotic behavior and complex dynamics for bigger values of theϕ parameter (see Figure 10(a)). Figures 11(a)–11(f) illustrate the space phase for different values of the parameter ϕ, showing chaotic dynamics and complicated attractors for different values of ϕ.Figure 10
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.Figure 11
Asymmetric case (R0>1, R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum for the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.22(b)
ϕ=0.261(c)
ϕ=0.55(d)
ϕ=0.8(e)
ϕ=0.99(f)
ϕ=2By varying the ADE parameter, the system shows a rich dynamical behavior for all the three case scenarios. For small values ofϕ, we have a fixed point; hence, the coexistence endemic equilibrium is stable. The Hopf bifurcation occurs for a specific range of ϕ values, from where disease propagation is possible and periodic solutions are observed. Complicated attractors and chaotic behavior are also observed.


## 5. Particular Case: Exponential Immunity

We extend our analysis to investigate a particular case of the proposed model. Here, we assume the length of immunity being exponentially distributed, that is, considering that the fraction of temporarily immune individuals remaining in theCi class is Pit=e−ωit, with ωi>0, for i=1,2. In this case, the integro-differential equation (IDE) system became an ODE system as follows:
(5)dStdt=d−dSt−β1StI1t+I21t−β2StI2t+I12t,dI1tdt=−d+γI1t+β1StI1t+I21t,dI2tdt=−d+γI2t+β2StI2t+I12t,dC1tdt=−dC1t+γI1t−ω1C1t,dC2tdt=−dC2t+γI2t−ω2C2t,dR1tdt=−dR1t−α2ϕR1tI12t+I2t+ω1C1,dR2tdt=−dR2t−α1ϕR2tI21t+I1t+ω2C2,dI12tdt=−d+γI12t+α2ϕR1tI2t+I12t,dI21tdt=−d+γI21t+α1ϕR2tI1t+I21t.This qualitative analysis of this particular case was only initially evaluated in Steindorf et al. [45]. Here, a detailed analysis of the dynamical behaviour is presented.


### 5.1. Numerical Results

The parameters used for the numerical experiments can be found in Table2, with ωi=2y−1.


#### 5.1.1. Bifurcation Structure

To analyse the stability of the coexistence endemic equilibrium (CEE), we evaluate the bifurcation structures shown in a diagram where the parameterϕ is the varying parameter. The stability for the CEE using ϕ as a bifurcation parameter is also shown numerically in [45]. In detail, as ϕ increases from small values towards the critical value, ϕc, the steady state changes from a stable focus to an unstable steady state. Therefore, the Hopf bifurcation occurs and a closed periodic orbit appears in a small neighbourhood of ϕc. In order to see the limit cycle around the equilibrium, in a small vicinity of the critical value, bifurcation diagrams are shown.
(i)
Symmetric CaseR0=R1>1Figures12(a) and 12(b) illustrate the bifurcation diagrams for the symmetric case scenario of β1=β2=180, with R0=3.46. The fixed point is found for small values of the parameter ϕ, after the chaotic behaviour is found (Figures 13(a)–13(b)).Figure 12
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the maximum and minimum values for the susceptible population are plotted.
(a)
Symmetric case. Hopf bifurcation atϕc=0.032(b)
Maximum and minimum values for the susceptible population forϕ∈0.1,1Figure 13
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.
(a)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,4(b)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,1Figures14(a)–14(f) show the space phase for different values of the parameter ϕ for the symmetric case scenario, where chaotic behaviour for ϕ<1 and limit cycles for larger values are observed.Figure 14
Symmetric case (R0=R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.01(b)
ϕ=0.032(c)
ϕ=0.2(d)
ϕ=0.6(e)
ϕ=2(f)
ϕ=3.5Figure15(a) illustrates the bifurcation diagram for the symmetric case scenario, with β=104 and R0=1.99. In this case, the Hopf bifurcation is found for the critical value of ϕc=0.08. Figures 16(a)–16(f) show the space phase for different values of the parameter ϕ. A limit cycle is found for ϕ=0.08, and chaotic behaviour for value of ϕ<1 is observed. For ϕ=3.5, more complicated limit cycles are observed, in this case with β=104, while for β=180, a limit cycle is seen.
(ii)
Asymmetric Case ofR0>1, R1>1Figure 15
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.Figure 16
Symmetric case (R0=R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.01(b)
ϕ=0.08(c)
ϕ=0.2(d)
ϕ=0.5(e)
ϕ=2(f)
ϕ=3.5Figures17(a)–17(b) illustrate the bifurcation diagrams for the case where β1=120 and β2=180. The fixed point is found for small values of parameter ϕ, after chaotic behaviour is observed (Figures 18(a) and 18(b)).Figure 17
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the maximum and minimum values for the susceptible population are plotted.
(a)
Case:R1=2.31. Hopf bifurcation at ϕc=0.244(b)
Maximum and minimum population forϕ∈0.26,1Figure 18
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the log of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.
(a)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,4(b)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,1Figures19(a)–19(f) show the space phase for different values of the parameter ϕ for the asymmetric case of β1=120 and β2=180. Chaotic behaviour for ϕ<1, limit cycle for ϕ=0.244, and complicated attractors after these values are shown.
(iii)
Asymmetric Case forR0>1, R1<1Figure 19
Asymmetric case (R0>1, R1>1). Space phase for different ϕ values. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.244(b)
ϕ=0.26(c)
ϕ=0.55(d)
ϕ=0.8(e)
ϕ=0.99(f)
ϕ=2Figures20(a)–20(a) illustrate the bifurcation diagrams for the case scenario where β1=45 and β2=180. Hopf bifurcation occurs for ϕ=2.52 (Figures 21(a) and 21(b)).Figure 20
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the maximum and minimum values for the susceptible population are plotted.Figure 21
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the log of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.
(a)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,4(b)
Hopf bifurcation atϕ=2.5Figures22(a)–22(f) show the space phase for different values of the parameter ϕ for the asymmetric case scenario of β1=45 and β2=180. Chaotic behaviour for ϕ>1, limit cycle for ϕ=2.52, and complex attractors after these values are shown.Figure 22
Asymmetric case (R0>1, R1>1). Space phase for different ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=2(b)
ϕ=2.1(c)
ϕ=2.3(d)
ϕ=2.52(e)
ϕ=2.8(f)
ϕ=3.5For the three case scenarios, the equilibrium is stable for a small range of theϕ parameter. While for the symmetric case scenario, the Hopf bifurcation occurs at ϕ=0.032, for the asymmetric case scenarios, the bifurcation points are ϕ=0.244 and ϕ=2.52, respectively, for scenario ii and scenario iii. Periodic solutions are found for critical values of ϕ, as it is shown in Figures 16, 19, and 22. For the symmetric case scenario, chaotic behaviour is observed for ϕ<1. For the asymmetric case scenarios, complex attractors are observed for wide parameter values below 1 and complicated limit cycles above 1, for case scenario ii, while chaotic behaviour is observed only for values above ϕ=2.52 in case scenario iii.


## 5.1. Numerical Results

The parameters used for the numerical experiments can be found in Table2, with ωi=2y−1.


### 5.1.1. Bifurcation Structure

To analyse the stability of the coexistence endemic equilibrium (CEE), we evaluate the bifurcation structures shown in a diagram where the parameterϕ is the varying parameter. The stability for the CEE using ϕ as a bifurcation parameter is also shown numerically in [45]. In detail, as ϕ increases from small values towards the critical value, ϕc, the steady state changes from a stable focus to an unstable steady state. Therefore, the Hopf bifurcation occurs and a closed periodic orbit appears in a small neighbourhood of ϕc. In order to see the limit cycle around the equilibrium, in a small vicinity of the critical value, bifurcation diagrams are shown.
(i)
Symmetric CaseR0=R1>1Figures12(a) and 12(b) illustrate the bifurcation diagrams for the symmetric case scenario of β1=β2=180, with R0=3.46. The fixed point is found for small values of the parameter ϕ, after the chaotic behaviour is found (Figures 13(a)–13(b)).Figure 12
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the maximum and minimum values for the susceptible population are plotted.
(a)
Symmetric case. Hopf bifurcation atϕc=0.032(b)
Maximum and minimum values for the susceptible population forϕ∈0.1,1Figure 13
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.
(a)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,4(b)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,1Figures14(a)–14(f) show the space phase for different values of the parameter ϕ for the symmetric case scenario, where chaotic behaviour for ϕ<1 and limit cycles for larger values are observed.Figure 14
Symmetric case (R0=R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.01(b)
ϕ=0.032(c)
ϕ=0.2(d)
ϕ=0.6(e)
ϕ=2(f)
ϕ=3.5Figure15(a) illustrates the bifurcation diagram for the symmetric case scenario, with β=104 and R0=1.99. In this case, the Hopf bifurcation is found for the critical value of ϕc=0.08. Figures 16(a)–16(f) show the space phase for different values of the parameter ϕ. A limit cycle is found for ϕ=0.08, and chaotic behaviour for value of ϕ<1 is observed. For ϕ=3.5, more complicated limit cycles are observed, in this case with β=104, while for β=180, a limit cycle is seen.
(ii)
Asymmetric Case ofR0>1, R1>1Figure 15
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.Figure 16
Symmetric case (R0=R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.01(b)
ϕ=0.08(c)
ϕ=0.2(d)
ϕ=0.5(e)
ϕ=2(f)
ϕ=3.5Figures17(a)–17(b) illustrate the bifurcation diagrams for the case where β1=120 and β2=180. The fixed point is found for small values of parameter ϕ, after chaotic behaviour is observed (Figures 18(a) and 18(b)).Figure 17
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the maximum and minimum values for the susceptible population are plotted.
(a)
Case:R1=2.31. Hopf bifurcation at ϕc=0.244(b)
Maximum and minimum population forϕ∈0.26,1Figure 18
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the log of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.
(a)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,4(b)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,1Figures19(a)–19(f) show the space phase for different values of the parameter ϕ for the asymmetric case of β1=120 and β2=180. Chaotic behaviour for ϕ<1, limit cycle for ϕ=0.244, and complicated attractors after these values are shown.
(iii)
Asymmetric Case forR0>1, R1<1Figure 19
Asymmetric case (R0>1, R1>1). Space phase for different ϕ values. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.244(b)
ϕ=0.26(c)
ϕ=0.55(d)
ϕ=0.8(e)
ϕ=0.99(f)
ϕ=2Figures20(a)–20(a) illustrate the bifurcation diagrams for the case scenario where β1=45 and β2=180. Hopf bifurcation occurs for ϕ=2.52 (Figures 21(a) and 21(b)).Figure 20
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the maximum and minimum values for the susceptible population are plotted.Figure 21
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the log of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.
(a)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,4(b)
Hopf bifurcation atϕ=2.5Figures22(a)–22(f) show the space phase for different values of the parameter ϕ for the asymmetric case scenario of β1=45 and β2=180. Chaotic behaviour for ϕ>1, limit cycle for ϕ=2.52, and complex attractors after these values are shown.Figure 22
Asymmetric case (R0>1, R1>1). Space phase for different ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=2(b)
ϕ=2.1(c)
ϕ=2.3(d)
ϕ=2.52(e)
ϕ=2.8(f)
ϕ=3.5For the three case scenarios, the equilibrium is stable for a small range of theϕ parameter. While for the symmetric case scenario, the Hopf bifurcation occurs at ϕ=0.032, for the asymmetric case scenarios, the bifurcation points are ϕ=0.244 and ϕ=2.52, respectively, for scenario ii and scenario iii. Periodic solutions are found for critical values of ϕ, as it is shown in Figures 16, 19, and 22. For the symmetric case scenario, chaotic behaviour is observed for ϕ<1. For the asymmetric case scenarios, complex attractors are observed for wide parameter values below 1 and complicated limit cycles above 1, for case scenario ii, while chaotic behaviour is observed only for values above ϕ=2.52 in case scenario iii.


## 5.1.1. Bifurcation Structure

To analyse the stability of the coexistence endemic equilibrium (CEE), we evaluate the bifurcation structures shown in a diagram where the parameterϕ is the varying parameter. The stability for the CEE using ϕ as a bifurcation parameter is also shown numerically in [45]. In detail, as ϕ increases from small values towards the critical value, ϕc, the steady state changes from a stable focus to an unstable steady state. Therefore, the Hopf bifurcation occurs and a closed periodic orbit appears in a small neighbourhood of ϕc. In order to see the limit cycle around the equilibrium, in a small vicinity of the critical value, bifurcation diagrams are shown.
(i)
Symmetric CaseR0=R1>1Figures12(a) and 12(b) illustrate the bifurcation diagrams for the symmetric case scenario of β1=β2=180, with R0=3.46. The fixed point is found for small values of the parameter ϕ, after the chaotic behaviour is found (Figures 13(a)–13(b)).Figure 12
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the maximum and minimum values for the susceptible population are plotted.
(a)
Symmetric case. Hopf bifurcation atϕc=0.032(b)
Maximum and minimum values for the susceptible population forϕ∈0.1,1Figure 13
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.
(a)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,4(b)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,1Figures14(a)–14(f) show the space phase for different values of the parameter ϕ for the symmetric case scenario, where chaotic behaviour for ϕ<1 and limit cycles for larger values are observed.Figure 14
Symmetric case (R0=R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.01(b)
ϕ=0.032(c)
ϕ=0.2(d)
ϕ=0.6(e)
ϕ=2(f)
ϕ=3.5Figure15(a) illustrates the bifurcation diagram for the symmetric case scenario, with β=104 and R0=1.99. In this case, the Hopf bifurcation is found for the critical value of ϕc=0.08. Figures 16(a)–16(f) show the space phase for different values of the parameter ϕ. A limit cycle is found for ϕ=0.08, and chaotic behaviour for value of ϕ<1 is observed. For ϕ=3.5, more complicated limit cycles are observed, in this case with β=104, while for β=180, a limit cycle is seen.
(ii)
Asymmetric Case ofR0>1, R1>1Figure 15
Symmetric case (R0=R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the logarithm of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.Figure 16
Symmetric case (R0=R1>1). Space phase for different values of parameter ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.01(b)
ϕ=0.08(c)
ϕ=0.2(d)
ϕ=0.5(e)
ϕ=2(f)
ϕ=3.5Figures17(a)–17(b) illustrate the bifurcation diagrams for the case where β1=120 and β2=180. The fixed point is found for small values of parameter ϕ, after chaotic behaviour is observed (Figures 18(a) and 18(b)).Figure 17
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the maximum and minimum values for the susceptible population are plotted.
(a)
Case:R1=2.31. Hopf bifurcation at ϕc=0.244(b)
Maximum and minimum population forϕ∈0.26,1Figure 18
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the log of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.
(a)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,4(b)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,1Figures19(a)–19(f) show the space phase for different values of the parameter ϕ for the asymmetric case of β1=120 and β2=180. Chaotic behaviour for ϕ<1, limit cycle for ϕ=0.244, and complicated attractors after these values are shown.
(iii)
Asymmetric Case forR0>1, R1<1Figure 19
Asymmetric case (R0>1, R1>1). Space phase for different ϕ values. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=0.244(b)
ϕ=0.26(c)
ϕ=0.55(d)
ϕ=0.8(e)
ϕ=0.99(f)
ϕ=2Figures20(a)–20(a) illustrate the bifurcation diagrams for the case scenario where β1=45 and β2=180. Hopf bifurcation occurs for ϕ=2.52 (Figures 21(a) and 21(b)).Figure 20
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the maximum and minimum values for the susceptible population are plotted.Figure 21
Asymmetric case (R0>1, R1>1). In the horizontal axis, the parameter ϕ varies, while in the vertical axis, the log of the maximum and minimum values for the infected population (logI1+I2+I12+I21t) is plotted.
(a)
Maximum and minimum values for the infected populationlogI1+I2+I12+I21t for ϕ=0,4(b)
Hopf bifurcation atϕ=2.5Figures22(a)–22(f) show the space phase for different values of the parameter ϕ for the asymmetric case scenario of β1=45 and β2=180. Chaotic behaviour for ϕ>1, limit cycle for ϕ=2.52, and complex attractors after these values are shown.Figure 22
Asymmetric case (R0>1, R1>1). Space phase for different ϕ. The susceptible population is plotted in the horizontal axis, while the logarithm of the sum of the infected population (logI1+I2+I12+I21t) is plotted in the vertical axis.
(a)
ϕ=2(b)
ϕ=2.1(c)
ϕ=2.3(d)
ϕ=2.52(e)
ϕ=2.8(f)
ϕ=3.5For the three case scenarios, the equilibrium is stable for a small range of theϕ parameter. While for the symmetric case scenario, the Hopf bifurcation occurs at ϕ=0.032, for the asymmetric case scenarios, the bifurcation points are ϕ=0.244 and ϕ=2.52, respectively, for scenario ii and scenario iii. Periodic solutions are found for critical values of ϕ, as it is shown in Figures 16, 19, and 22. For the symmetric case scenario, chaotic behaviour is observed for ϕ<1. For the asymmetric case scenarios, complex attractors are observed for wide parameter values below 1 and complicated limit cycles above 1, for case scenario ii, while chaotic behaviour is observed only for values above ϕ=2.52 in case scenario iii.


## 6. Conclusions

While the qualitative analysis of the particular case considering exponential immunity was initially evaluated in [45], here we have performed a detailed analysis for the ADE parameter, gaining insights into the dynamical behavior of disease propagation in a population.In this study, we have shown that the model has an equilibrium of the coexistence of serotypes. The mechanisms that can lead to this finding are the cross-immunity protection period and the ADE effect. These mechanisms have an important role in dengue modelling, since the disease exhibits different serotypes cocirculating in a given region. It was also shown that, if the ADE effect is small enough, there is a range ofϕ values where if one serotype is endemic; the other serotype will not be able invade the population. In that case, the presence of one serotype will eventually protect the population from the invasion of the other serotype, decreasing the risk of a secondary infection to occur and hence decreasing severe disease in that population. Otherwise, the serotypes will coexist.If the key parameter describing the ADE reaches the threshold value, then a supercritical Hopf bifurcation occurs, leading to periodic solutions. If the parameter is big enough, then the endemic equilibrium is unstable and a chaotic behaviour occurs. Within this parameter region where the complex dynamics scenario is observed, the disease propagation can only be predicted for a limited period of time.Results showing complex attractors are of major importance to describe dengue fever epidemiology dynamics, since the available empirical data dynamics show irregular behaviour resembling chaotic dynamics, as previously described by Aguiar et al. [26]. We note that this study is an ongoing work, with many aspects still to be investigated. As an example, Lyapunov exponent calculations will certainly complement the bifurcations diagram analysis and phase space plots, able to show the exact parameter value where the chaotic behaviour starts and gets eventually stabilized.Comparing the particular case of the exponential distributed function for the immunity period, with the general delay model, we observed that the qualitative behaviour of the system was not altered by the choice of the function. However, the invasion reproduction number depends on the average cross-protection time, which can affect whether the infection could coexist or not. Moreover, we observe slight differences in the numerical experiments for the bifurcations at the critical values, as well as for the range of periodic solutions and complex attractors.



---

*Source: 2074325-2022-02-13.xml*