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

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) e ﬀ ect, 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 di ﬀ erent scenarios of dengue fever transmission in a population. The model shows a rich dynamical behavior, from ﬁ xed points and periodic oscillations up to chaotic behaviour with complex attractors.


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 crossimmunity (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][4][5][6][7][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][14][15][16][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 hostto-host transmission. In the first approach, the fluctuations  ⋆ 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 R 0 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 R 0 from 1:38 to 7:70, with average 3:2, using data from Thailand. in the mosquito dynamics and climate change are included and assumed to affect the disease transmission [23][24][25][26][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][29][30][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][35][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 3 Computational and Mathematical Methods 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 crossprotection 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 con-sider 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 crossprotection 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.  Computational and Mathematical Methods 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 antibodydependent 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.

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.

Computational and Mathematical Methods
Let NðtÞ 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 (SðtÞ), infected by one serotype i (I i ðtÞ), temporarily immune to all serotypes after infection by serotype i (C i ðtÞ), and recovered from serotype i, and after the temporary cross-immunity protection period, they become susceptible to the other serotype, R i ðtÞ. Secondary infection can occur with a heterologous serotype; that is, I ij ðtÞ 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 RðtÞ 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     [46][47][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 I i ðtÞ 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 C i ðtÞ, due to the crossimmunity protection period, with P i ðtÞ 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 P i ð0Þ = 1 and P i ð∞Þ = 0 and P i ðtÞ is not increasing. Moreover, Thus, the number of cross-protected individuals is given by where C 0 is the number of protected individuals at time t = 0. Of course, C 0 must satisfy lim t⟶∞ C 0 ðtÞ = 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 in R i 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 C 0 = 0, the mathematical model can be written as follows:  Phase space = 3.5

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.
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 Table 1.

Numerical Results
The numerical values for the parameters are shown in the Table 2.

Computational and Mathematical Methods
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),

Stability of the Coexistence Endemic Equilibrium (CEE).
The stability of the CEE will be studied using numerical experiments. The CEE, D 3 , will be in the positive region only if the parameters satisfy the condition R Inv > 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: In this case, the infection rates are the same for both strains. For the simulations, we fixed β i = 180, giving the value for R 0 = 3:46. For the symmetric case scenarios, the   13 Computational and Mathematical Methods 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 for R 0 > 1, R 1 < 1 In this case scenario, with values for the parameters in Table 2, we fixed β 1 = 45 and β 2 = 180, which gives R 1 = 0:87 and R 0 = 3:46. To study the stability of the CEE, it is necessary that the condition R Inv > 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) . Figures 2(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).
The values for the parameters used in the simulations for this case scenario can be found in Table 2. We fixed β 1 = 120 and β 2 = 180, which gives R 1 = 2:31 and R 0 = 3:46, respectively, and thus R Inv > 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). Figures 4(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.

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 Case R 0 = R 1 > 1 Figure 6(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 for R 0 > 1, R 1 < 1 The 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 for R 0 > 1, R 1 > 1 After the Hopf bifurcation, the diagram shows chaotic behavior and complex dynamics for bigger values of the ϕ parameter (see Figure 10 By 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.

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  Phase space = 0.08 15 Computational and Mathematical Methods in the C i class is P i ðtÞ = e −ω i t , with ω i > 0, for i = 1, 2. In this case, the integro-differential equation (IDE) system became an ODE system as follows: 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.

Numerical
Results. The parameters used for the numerical experiments can be found in Table 2, with ω i = 2y −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.     12(b) illustrate the bifurcation diagrams for the symmetric case scenario of β 1 = β 2 = 180, with R 0 = 3:46. The fixed point is found for small values of the parameter ϕ, after the chaotic behaviour is found (Figures 13(a)-13(b)). Figures 14(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 15(a) illustrates the bifurcation diagram for the symmetric case scenario, with β = 104 and R 0 = 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.

Figures 17(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)). Figures 19(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. For 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.

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.

Data Availability
Not applicable.

Disclosure
Vanessa Steindorf's present address is the Basque Center for Applied Mathematics (BCAM), Bilbao, Spain.

Conflicts of Interest
The authors declare that they have no conflicts of interest. 20 Computational and Mathematical Methods