Cross immunity protection and antibody-dependent enhancement in a distributed delay dynamic model

1.   Introduction

Dengue fever is an endemic disease in tropical and subtropical countries. According to the World Health Organization ^[1] in 2016, more than 2.38 million cases were reported in the Americas, where Brazil alone contributed with almost 1.5 million cases. In 2017, a significant reduction in the number of cases in the Americas was reported, even though a recent estimate indicates 390 million cases per year around the world ^[1].

Dengue fever is caused by dengue virus (DENV), which is transmitted to humans by the bite of infected mosquitos. It is known that there are four distinct virus serotypes circulating in the population ^[2]. The disease presents as Dengue classic or Dengue hemorrhagic fever (DHF).

Antibody-dependent enhancement (ADE) was proposed to explain the more frequent occurrence of the severity of Dengue fever and Dengue hemorrhagic fever in secondary infections ^[2,3,4]. Another characteristic of dengue virus is cross-immunity protection so that infection by any of the four serotypes leads to short-term protection for all serotypes. Long-term protection occurs only for the strain with which the individual was infected ^[4].

A number of mathematical models that include ADE were proposed in ^[5], ^[6], ^[7] and ^[8]. Adams et al. ^[5] included in a two-strain model the ADE effect and, described the immunological distance between dengue serotypes through a function with the hypothesis that this function reduces the probability of contracting a secondary infection. Numerical results were presented in this study. In a host-host model, Billings ^[7] demonstrated that ADE causes oscillations. Bianco ^[6] investigated a model with partial temporary cross immunity and showed that weak cross immunity stabilized the system while strong cross immunity caused chaos. While Hu ^[8] considered vector dynamics, revealing that ADE can explain data complexity.

Aguiar et al. ^[9] proposed an Ordinary Differential Equations (ODE) system for two serotypes of the Dengue virus. The total population was stratified by their clinical states. Time series simulations were plotted, and various bifurcation phenomena were observed. Theoretical mathematical analysis of the model in ^[9] was carried out in ^[10] and in ^[11]. Aguiar et al. ^[10], quantified the attractor structure, limit cycle and the chaotic attractor by calculating Lyapunov exponents. Analytic formulations for the equilibrium and analysis of the bifurcation structure were obtained by Kooi et al. in ^[11] for the symmetric case (when the strains were the same).

These aforementioned models were formulated in terms of ordinary differential equations, without incorporating time lags for recovery and incubation. General delay models were introduced by Van den Driessche ^[12], using a special step function for the recovery period. Recently, Nah ^[13] described a mathematical model for malaria transmission with a time delay involving an exposure and a general incubation period.

A mathematical model describing a general multi-strain disease was proposed in ^[14]. The authors proposed a delay diffusive two strain disease model, considering the SIR structure, the constant recruitment rate, and the constant time delay representing the length of the immunity period. The stability of the model was determined by the Basic Reproduction Number. Guan ^[15] also described a Dengue fever model with a time delay which refers to the time of incubation of the virus in an infected population class and in an infected mosquito population class. In a SIR epidemic model, Hattaf ^[16] also included a constant time delay.

Taking into account a simpler SIR model, Huang ^[17] considered an infinite distributed delay on a complex population network. Numerical experiments confirmed that the delay slows down the extinction of the disease when the threshold value is smaller than 1, while the delay accelerates the spread if it is bigger than 1. Distributed delay was also used by Xu ^[18] in a SVEIR model involving a continuous vaccination strategy. Their results showed that distributed delay had no impact on the qualitative behavior and global dynamics.

Since dengue fever constitutes a public health problem, we proposed a mathematical model, namely, an IDE system that can be applied to describe and study the propagation of dengue fever in a population, considering two main characteristics of the disease: ADE and cross-immunity protection. The main purpose was to include and analyze the effect of the general time delay on the model, which is included in order to represent the length of cross immunity protection. In addition, a constant parameter will be added to the model in order to represent the effect of ADE.

In the third section, the model is analyzed. By defining the associated limiting system, we look for equilibrium and show the local stability that is determined by an important threshold value. In Section 4, by constructing a Lyapunov functional, we prove global stability. Numerical experiments are performed in Section 5, where the bifurcation structure and stability of the coexistence equilibrium are numerically studied. Final considerations and conclusions are made in the final section.

2.   Model formulation

In this section we shall propose a model to describe a multi-strain disease, motivated by dengue fever.

Let $N(t)$ be the total population of individuals at time $t$ in a region. We divided the population into disjoint classes according to the individual status, susceptible for all serotypes, infected by serotype $i$, temporarily immune for all serotypes after being infected by serotype $i$, recovered for serotype $i$ but susceptible to the others, and, recovered for all serotypes, represented at time $t$, respectively by $S(t)$, $I_i(t)$, $C_i(t)$, $R_i(t)$ e $R(t)$. Furthermore, we included two more classes for the population, $I_{ij}(t)$, representing the subpopulation reinfected by serotype $j$ after being infected by serotype $i$, with $i, j = 1, 2$.

Assuming that the birth and mortality rates are equal, define $d$ as the human population's constant natural mortality rate. We assumed that the constant rate, $\beta_{i}$ is the transmission rate by serotype $i$, for the first time an individual is infected. Whereas $\alpha_j$ is the transmission rate by different serotypes, by serotype $j$, for the second time an individual is infected.

The mosquito population is not explicitly considered in the model, so the $\beta_i$ and $\alpha_j$ rates will represent the average number of bites and the probability of a susceptible individual being bitten by an infected mosquito by serotype $i$ and the probability of a recovered individual being bitten by an infected mosquito by serotype $j$, respectively.

Individuals in the infectious classes $I_i(t)$, remain in this class with an average time $\frac{1}{\gamma}$ since we assumed that the length of this class is exponentially distributed. Once infected, the individual recovers and moves on to the temporary immunity class ($C_i(t)$. In this class, the individual gains temporary immunity to all serotypes.

Moreover, we assumed for $C_i$ class a general length of immunity. Let $P_i(t)$ be the fraction of individuals that, after recovering from the serotype $i$, remain cross protected by all serotypes, $t$ units after get in the temporary immunity class. It is reasonable to assume that

and, satisfied that $P_i(t)$ is non-increasing and

Thus, the number of cross protected individuals is given by

where $C_0$ is the number of protected individuals at time $t = 0$ that are still immune at time $t$. Of course, $C_0$ must satisfy $\lim_{t \rightarrow \infty } C_0(t) = 0$.

After the cross-protection time the individual remains susceptible to the other strains. In this way, an individual in the $R_i$ class can be infected again at rate $\alpha_{j}$ by a different serotype. Only after being infected by all strains do the individual become immune, recovering at a constant rate $\gamma$, remaining permanently recovered.

Whenever the individual is infected again, by a different serotype, the immune system responds, increasing the infectiousness and enhancing viral replication ^[2,3,19]. This effect is called ADE, when secondary infections are possible in the presence of an intermediate levels of antibodies, increasing viral replication ^[2,3,19].

In the studies in ^{[6,7,8,9,10,11]} the authors assume that this increasing viral load is associated with increasing transmission infections, and they use this assumption to mathematically describe and model the ADE using a constant $\phi$ that represents the increase in the transmission rate, associated with the transmission of the infected individuals during secondary infections.

According to Katzelnick ^[19], ADE occurs in a specific range of antibody concentrations. Low levels of antibodies did not enhance disease, intermediate levels exacerbated disease, and high levels protected against severe disease. The authors verified enhancement in humans and showed that pre-existing antibodies were associated with the severity of the disease, and they also showed that the immune correlate for enhanced severe dengue is distinct from that for protection. Furthermore, Rothman ^[20] affirms that depending on the specific antibody concentration, dengue virus antibodies can inhibit viral infection (neutralization) or enhance infection.

Therefore, ADE takes place or not, increasing or decreasing infection, due to a high or low antibody concentration in people who were infected once (recovered) and again infected. Clearly, this high viral load is not related to the increasing transmission (from a recovered individual to new infections) but rather to the capacity of the immune response system of the primary individual to respond to the secondary infection and then to be protected against infection or enhance the disease.

Thus, the ADE characteristic must be described as a constant rate that can increase or decrease the probability of the recovered individual who has already been infected once to be infected again, assuming that the individual with the primary infection carries a viral load that may neutralize or enhance the infection. Therefore, in this study, we proposed that $\phi$ represents the fraction that decreases ($\phi < 1$) or increases (enhances) ($\phi > 1$) the probability of secondary infections in primary infected individuals. In other words, the epidemiological consequence ADE will be described through the constant coefficient $\phi$, where it is assumed that previous exposure to one serotype results in an increase in the susceptibility for reinfection.

Thus, based on the assumptions and in the studies by Wang et al. ^[21], Cooke et al. ^[22], Van den Driessche et al. ^[23], and taking the derivative of equation (2.3), assuming $C_0 = 0$, the model can be described as follows:

The total population dynamics is determined by

and, since ${\frac{dN}{dt} = 0}$, the total population remains constant in time.

Denote $\frac{S}{N} = S$, $\frac{I_{ij}}{N} = I_{ij}$, $\frac{C_i}{N} = C_i$ and $\frac{R_i}{N} = R_i$ representing, for each class, the fractions of the population, and $N = N^*$. Thus, the sum of the total population satisfies $S+ I_1 + I_2 + C_1 + C_2 + R_1 + R_2 + I_{12} + I_{21}+ R = 1$. In addition, the dynamics of the recovered and cross-immunity classes are decoupled. Therefore, the original system can be studied by analyzing the following subsystem:

This subsystem will be qualitatively analyzed in Section 3. In the next subsection, we are going to discuss a particular case of this model.

2.1.   A particular case: Exponential immunity

In this subsection, we are going to discuss a particular case of the model (2.4). If we assume the length of immunity is exponentially distributed, which means that the fraction of individuals temporarily immune remain in the $C_i$ class is $P_i(t) = e^{- \omega_i t}$, with $\omega_i > 0$, for $i = 1, 2$, then the IDE system (2.4) becomes an ODE system. Thereby, the system already normalized with the variables representing the fractions of the populations can be described as follows:

Qualitative analysis was conducted for this particular model through classic theory for ODE systems. Details are found in the supplementary material. We summarize the results here, but first, we shall define a threshold value that is very important for the stability and existence of equilibrium.

2.1.1.   Basic Reproduction number and Invasion Reproduction number

The Basic Reproduction number $\mathcal{R}_0$ is defined by many authors, as in ^[24], as the expected number of secondary infections produced by one case in a susceptible population and also as a measure of the potential for disease spread in a population. Mathematically, the Basic Reproduction Number is a threshold for the stability of the Disease Free equilibrium ^[24].

We shall define $\mathcal{R}_1 = \frac{\beta_1}{d + \gamma}$ as the Basic Reproduction Number for strain one and, $\mathcal{R}_2 = \frac{\beta_2}{d + \gamma}$ as the Basic Reproduction Number for strain two.

It is usual to define an overall reproduction number for a multi-strain model, with different strains. Thus, it will be defined as:

This definition, according to Martcheva ^[25], is similar to that given in ^[24], where the Basic Reproduction number is defined mathematically as the spectral radius (the maximum of the modulus of the eigenvalues) of the Next Generation matrix.

When working with a multiple-strain model, it is usual to find another important threshold value. This value will determine whether the Coexistence Endemic equilibrium (CEE) will be in the biological positive region and also the stability of the Boundary Equilibrium (BE). This threshold is called the Invasion Reproduction number (of strain one at the equilibrium of strain two) ^[25], and we shall denote it by $\mathcal{R}_{Inv}$. Martcheva ^[25] defined epidemiologically as a number of secondary infections that one individual infected with strain one will produce in a population, in which strain two is at equilibrium during its infectious lifetime.

2.2.   Discussion of qualitative result

For the ODE model (2.6) we proved the existence of four equilibria (proofs and details can be found in the supplementary material). The Disease Free equilibrium (DFE) is locally asymptotically stable if $\mathcal{R}_0 < 1$, and thus the disease will die out. However, the DFE is unstable if $\mathcal{R}_0 > 1$.

We found two BEs, each equilibrium corresponding to only one infection. The infection with the lowest $\mathcal{R}_i$ value is always unstable whereas the other is stable if $\mathcal{R}_{Inv} < 1$. This means that the strain with the biggest force of infection will persist, and the other will die out. In addition, the BE is unstable if $\mathcal{R}_{Inv} > 1$ and thus there is a fourth equilibrium, the CEE, in the positive invariant region. In this scenario, the two strains will coexist.

This ODE model is very similar to the models proposed in ^[9], ^[10], ^[11]. However, the assumption for the ADE effect is different, leading to different dynamics. The analytical result for the CEE and the stability was obtained for the symmetric case in ^[10] when all the parameters were assumed to be equal for different strains. The authors numerically showed the transcritical bifurcation for the asymmetric case in ^[11] as well as a robust diagram showing different bifurcations for the epidemic model, showing the complexity of the referred model. Nonetheless, there is no analytical form of the CEE as we showed for this model in this paper.

We could prove that for certain values of the parameter that represents the enhancement, there is indeed a CEE within the region. This is a very important result because, according to VinodKumar ^[26], different serotypes have been co-circulating in the same area with one of them being dominant during an outbreak.

We also numerically proved that the equilibrium can be stable only for small parameter values that represent the ADE. Additionally, we showed that bifurcation occurs and, for most of the values of the parameter $\phi$ that represent the ADE, equilibrium exists and it is unstable. Thus, the coexistence of the two serotypes is possible but not established, and the solution oscillates for a critical parameter value.

3.   General Case: Time delay system

In the previous section (2.2), we discussed a particular case, considering the time of cross-immunity protection as exponentially distributed. We shall now analyze the model assuming the cross-immunity protection time as any continuous function with certain additional properties.

Consider the system (2.5) together with the assumptions (2.1) and (2.2). Following the idea in ^[27], we will examine the system (2.5) as a perturbation of the limiting system:

According to Miller ^[28], the limiting system ensures that the initial system proposed has equilibrium. In adittion, according to Hethcote ^[29], using Miller's theorems ^[30], it is possible to show that the equilibrium of the system coincides with that of their limiting equation.

For this limiting system (3.1) it is necessary to define the Banach space with memory, as defined in ^[31], ^[32], ^[33] and ^[34], in order to have a well-posed system and to have solutions defined in $(-\infty, 0]$. For qualitative theory, it will be useful to consider the following space phase.

Denote $\Delta_1$, $\Delta_2$ positives constants such that $\Delta_i < d$, satisfying, for $i = 1, 2$, $\int_{0}^{\infty} P_i(u)e^{-du}e^{\Delta_i u} du < \infty$. Define the Banach space, for $i = 1, 2$,

where $||\Psi||_{e} = \mbox{sup}_{s \leq 0} | \Psi(s)| e^{\Delta_i s}$. We consider $X = \mathbb{R} \times X_{\Delta_1} \times X_{\Delta_2} \times \mathbb{R}^7$ as the phase space for the limiting system (3.1).

Denote $I_{i_t}$, for $i = 1, 2$, the solution $I_i(t)$ at time $t$, that is $I_{i_t}(s) = I_i(t + s)$, $s \leq 0$. Let $\Lambda_i = \{ \Psi \in X_{\Delta_i} : \Psi(s) \geq 0, s \in (-\infty, 0] \ \}$, $i = 1, 2$. Then, for initial conditions, $S(0) = s_0 \in \mathbb{R}_+$, $I_{1_0} = \Psi_1 \in \Lambda_1$, $I_{2_0} = \Psi_2 \in \Lambda_2$, $C_i(0) = c_i \in \mathbb{R}_+$, $R_i(0) = r_i \in \mathbb{R}_+$, $I_{ji}(0) = \theta_i \in \mathbb{R}_+$, $i = 1, 2$, the solutions of the limiting systems in $X$ remain non-negative and, $I_{i_t} \in X_{\Delta_i}$, for all $t$, for $i = 1, 2$. Moreover, $\Omega _{X} = \{ (S, I_1(.), I_2(.), C_1, C_2, R_1, R_2, I_{12}, I_{21}, R) \in \mathbb{R}_+ \times \Lambda_1 \times \Lambda_2 \times \mathbb{R}_+^{7} : S+ I_1(0)+ I_2(0)+ C_1+ C_2+ R_1+ R_2+ I_{12}+ I_{21}+ R \leq 1 \}$ is positively invariant for system (3.1).

The trivial equilibrium $D_0 = (1, 0, 0, 0, 0, 0, 0, 0, 0, 0)$ is always in the invariant set $\Omega_X$. Now using the assumption

we get $0 < - \int_0^\infty P_i'(u)e^{-d(u)}du < \infty$. Thus, we will define

Furthermore, we also assume that

The assumption (3.4) together with (3.2) leads to

Since $0 < M_i < 1$, in the case of the extinction of one of the strains, the BE of the system (3.1) are

and,

and it will be in the $\Omega_X$ region, as long as the parameters satisfy ${\frac{{\beta_1}}{d + \gamma}} > 1$ and, ${\frac{\beta_2}{d + \gamma}} > 1$, respectively.

In the case of coexistence of the two strains, the CEE is given by

and, $S^*$ is the root of the cubic polynomial $O(S) = b_3S^{3} + b_2S^{2} + b_1 S + b_0$ where

Moreover, $S^*$ has to satisfy that $S^* < \frac{d+ \gamma}{{\beta_i}}$, $i = 1, 2$ in order to have the equilibrium in $\Omega_X$. This discussion proves the following theorems about the equilibria of the system (3.1).

Theorem 1. If $\frac{{ \beta_1}}{d + \gamma} > 1$ then the system of equations (3.1), always has a BE, $D_1$, in $\Omega_X$.And, if $\frac{{\beta_2}}{d + \gamma} > 1$ then the system of equations (3.1), always has a BE, $D_2$, in $\Omega_X$.

Theorem 2. Without loss of generality, assume $\beta_2 > \beta_1$. If max $\{\frac{{\beta_1}}{d + \gamma}, \frac{\beta_2}{d + \gamma} \} > 1$ and,

then, the system (3.1) admits an equilibrium of the coexistence with two strains (a CEE).

Proof. The independent term, $b_0$, of the polynomial $O(S)$ is always positive. Since the equilibrium is given by (3.8) with $S^*$ being a root of the polynomial $O$, and this equilibrium will be in the region $\Omega$ if $S^* < \frac{d + \gamma}{{\beta_i}}$, for $i = 1, 2$ let us define

Then, if $\frac{{\beta_2}}{d + \gamma} > 1$ and $\mathcal{R}_{Inv} > 1$ then

This shows that we have a root $S^*$ of the polynomial $O$, such that $0 < S^* < S_{min}$ for $i = 1, 2$. Therefore, for this $S^*$, we have a positive equilibrium in $\Omega_X$ with the coexistence of the two strains.

3.1.   Stability analysis

In Section 3 we calculated the equilibria of the limiting system in order to know the equilibria of the delay system. In this section we shall introduce results that connect the local stability of the limiting system with the local stability of the delay system.

For the purpose of finding the stability of the solutions of the system (2.5) we follow the study conducted by Brauer ^[27], denoting the limiting system (3.1) as the unperturbed system.

So, we regard

as a perturbation, where $H$ is a vector that contains functions, $Y$ is the matrix containing the variables of the populations, $P'$ is the vector containing the functions $P_i'(s) e^{-d(s)}$. This perturbation function tends to zero as $t$ goes to infinity. Adding the perturbation (3.10) to the limiting system (3.1) we have the initial system (2.5).

The stability of the equilibria of the limiting system is a consequence of stability of the zero solution of the linearized system. Moreover, the asymptotic stability of the zero solution of the linear system $X'(t) = A X(t) + \int_0^{\infty} B(s)X(t-s)ds$ is equivalent to finding no solutions in the right half of the plane $Re \lambda \geq 0$ of $det (\lambda I - A - \hat B(\lambda)) = 0$, where $I$ is the identity matrix, $A$ is a matrix, and $\hat B (\lambda)$ denotes the Laplace transform of $B$ ^[35].

Since the assumptions (3.5) hold, and the perturbation function of the system is integrable, we have the necessary assumptions to use Theorem 2 in ^[27]. First, we need the following results of the stability of the limiting system (3.1), mathematically defining $\mathcal{R}_0 = \mbox{max} \{ \mathcal{R}_1, \mathcal{R}_2\} = \mbox{max} \{ \frac{\beta_1}{d + \gamma}, \frac{\beta_2}{d + \gamma} \}$.

Theorem 3. If $\mathcal{R}_0 = max \{ \mathcal{R}_1, \mathcal{R}_2 \} < 1$ then the DFE ($D_0$), of the system (3.1) is locally asymptotically stable. And, $D_0$ will be unstable if $\mathcal{R}_0 > 1$.

Proof. Solving the characteristic equation

where $I$ is the identity matrix $9 \times 9$, $H_0$ is the matrix

and,

with ${\hat G_i(\lambda) = \int_0^{\infty} e^{-\lambda z}g_i(z) dz = \int_0^{\infty} e^{-\lambda z}P_{i}'(z)e^{-d z}} dz$, gives us the following eigenvalues of the system at DFE:

If $\mathcal{R}_0 < 1 \ $ then all the eigenvalues are negative. If $\mathcal{R}_0 > 1$ then at least one eigenvalue will be positive. This guarantees the proof.

Remark: Without loss of generality, we are assuming that the infection has different forces, and therefore we assume for now on that $\beta_2 > \beta_1$.

Theorem 4. The BE, $D_1$, of the system (3.1) is always unstable in the $\Omega_X$ region, and the BE, $D_2$, is locally stable, in $\Omega_X$, if $\mathcal{R}_{Inv} < 1$ and is unstable in $\Omega_X$ if $\mathcal{R}_{Inv} > 1$.

Proof. Solving the characteristic equation

where $I$ is the identity matrix $9 \times 9$, $H_1$ is the matrix

where $w_1 = \frac{{\alpha_2}\phi}{{\beta_1}}\gamma M_1(\mathcal{R}_1 -1)$, $w_2 = \frac{{\alpha_1}\phi}{{\beta_1}}d (\mathcal{R}_1 -1)$, and $\hat G(\lambda)$ is the matrix whose elements are the Laplace transform of the $g_i(z) = e^{-d z}P'_{i}(z)$ of the associated linear system gives us the following eigenvalues of the system at $D_1$, given by (3.6):

Since $\beta_2 > \beta_1$, the eigenvalue $\lambda_9$ will be always positive. On the other hand, solving the characteristic equation of the associated linear system in $D_2$, given by (3.7), gives us the following eigenvalues:

Therefore, all the eigenvalues have a negative real part, except for one eigenvalue that changes its sign. If $\mathcal{R}_{Inv} < 1$ the eigenvalue $\lambda_9$ is negative. On the other hand, if $\mathcal{R}_{Inv} > 1$ thus, $\lambda_9$ is positive. This proves the result.

Remark: Since $M_2 < 1$ then, if ${\frac{{\alpha_1} \phi}{{\beta_1}} \leq 1}$, we have that $\mathcal{R}_{Inv} < \mathcal{R}_1$. Biologically speaking, this results means that there is a range of values for $\beta_1$ which strain one cannot invade the population if strain two is endemic. In this way, infection force two may protect the population from infection force one. After this range of values, the strains coexist.

3.1.1.   Stability of the solutions of the Time Delay System

In Section 3 we showed that the equilibria of the unperturbed system (3.1) and the theorems proved the uniformly asymptotic stability of the zero solution of the linear limiting system, and therefore, the stability of the equilibria of the limiting system. Thus, using Theorem 2 ^[27] we have the following results on the stability of solutions of the initial system (2.5) and, therefore for (2.4).

Corollary 1. If $\mathcal{R}_0 < 1$, thus the DFE, $D_0$, of the system (2.5) is locally asymptotically stable. And $D_0$ is unstable if $\mathcal{R}_0 > 1$.

Corollary 2. The BE, $D_1$, of the system (2.5) is always unstable in $\Omega$. And, the BE, $D_2$, is locally stable in $\Omega$ if $\mathcal{R}_{Inv} < 1$ and is unstable in $\Omega$ if $\mathcal{R}_{Inv} > 1$.

Remark: The analysis of the local stability of the CEE using this theory was not successful since we have to deal with a characteristic transcendental equation, with an infinite number of roots. Thus we will study the stability of the CEE numerically in Section 5.

4.   Global stability

In this section we shall investigate the global dynamics of the system (2.5) by constructing a suitable Lyapunov functional.

4.1.   Global stability of the DFE

Theorem 5. If $\mathcal{R}_0 \leq 1$, then the DFE, $D_0$, of the system (2.5) is globally attractive in $\Omega$.

Proof. First, denote the positive function, for $i = 1, 2$,

Consider $h(z) = z - 1 - ln (z), z \in \mathbb{R}_+^*$. Then, $h \geq 0$.

Let $L$ be the function

so,

is well-posed, is a Lyapunov functional, with $L \geq 0$, where the equality being true if and only if $S = 1$ and, $I_i = 0$, $R_i = 0$, $I_{ij} = 0$ for $i, j = 1, 2$, since $h\left(S\right)\geq 0$, for any $S > 0$.

Now, differentiating $L$ along the solution of the system (2.5) we have

Since $S + \frac{1}{S} \geq 2$, we have that the term in (4.3) is always non-positive in $\Omega$. Also, since $\mathcal{R}_0 \leq 1$, we have that the term in (4.4) is always non-positive in $\Omega$. $R_1$ and $R_2$ are in $\Omega$ and, $\Pi_1$ and $\Pi_2$ are positive functions of $t$, thus the term in (4.5) is always non-positive. The last term in (4.6) is non-positive because $P_i(t)$ is decreasing, and therefore, $P_i'$ is negative. Therefore, $L' \leq 0$ in $\Omega$ and, the equality is true if and only if each term of the equation is zero. From (4.3) the equality is true if and only if $S = 1$ and, with $\mathcal{R}_0 < 1$, from (4.4) to (4.6), we conclude that $I_1 = I_2 = R_1 = R_2 = C_1 = C_2 = I_{12} = I_{21} = 0$.

In the case that $\mathcal{R}_0 = 1$, since $S = 1$, from the first equation of the system we have that $S' = -\beta_i (I_i + I_{ji}) < 0$. But this is a contradiction with the fact that, since $S = 1$, it implies that $S' = 0$. Thus, $I_i + I_{ji} = 0$. Therefore, for instance, if $\mathcal{R}_0 = 1$, we still have that $I_2 = 0$ and $I_{12} = 0$. In this way, defining $E = \{ (S, I_1, I_2, C_1, C_2, R_1, R_2, I_{12}, I_{21}, R) \in \Omega \ \; \ L'(t) = 0 \}$, thus, the singleton ${D_0}$ is the largest invariant set in $E$. By the Invariance Principle for IDE ^[37] we have that the DFE, $D_0$, of the system (2.5) is globally asymptotically stable in $\Omega$.

4.2.   Global stability of the Boundary Equilibrium

We proved, in the previous section (3.1.1), that for $\beta_2 > \beta_1$ the BE, $D_2$, is locally asymptotically stable, when $\mathcal{R}_2 > 1$ and $\mathcal{R}_{Inv} < 1$. Now we are able to show that under these same conditions the trajectories with initial conditions in $\Omega - \{(S, I_1, I_2, C_1, C_2, R_1, R_2, I_{12}, I_{21}, R) \ \in \Omega \; I_2 = 0 \}$ approach this equilibrium. We need, first, the following result.

Lemma 1. If $\mathcal{R}_2 > 1$ and $\mathcal{R}_{Inv} < 1$ then the trajectories of system (2.5) that start in $\Omega$ approach the invariant set $\Omega_2 = \{(S, I_1, I_2, C_1, C_2, R_1, R_2, I_{12}, I_{21}, R) \in \Omega \ \; I_1 = C_1 = R_1 = I_{12} = I_{21} = R = 0\}$.

Proof. Consider the Lyapunov functional

with $\Pi_1$ being the same function definite in (4.1).

Thus, differentiating $L_1$ along the solution of the system (2.5) we have

Therefore, $L_{1}' \leq 0$ in $\Omega$ and, the equality is true if, and only if, $R_1 = 0$, $I_{12} = 0$, $\Pi_1 = 0$ and $- \int_t^{\infty} \gamma I_1(t)P_1'(t-s)e^{-d(t-s)ds} = 0$. In addition, from the equation of the system and, from the equality above, we have directly that $I_1 = 0$ and, $C_1$ goes to zero when t goes to infinity. Also, since $I_1 = 0$ we have that

Then, $S = 0$ or $I_{21} = 0$. If $S = 0$ then, for the first equation of the system, we have $S'(t) = d > 0$, which is a contradiction, since $S = 0$ implies $S' = 0$. Then $I_{21} = 0$. Therefore, we have shown that the maximal invariant set contained in the set of all points in $\Omega$, where $L_1^{'} = 0$, is $\Omega_2$. This shows the lemma.

This lemma shows that under certain conditions is sufficient to study the dynamics of the delay system (2.5) only on the projection of $\Omega_2$, this is, in the set $\Omega_{2_p} = \{(S, I_2, C_2, R_2); \ \ S+I_2+C_2+R_2 \leq 1)$ in $\mathbb{R}^4$. In this set, the initial system is reduced to

This system has two equilibria $D_{0_p} = (1, 0, 0, 0)$ and $D_{2_p} = (\frac{d + \gamma}{ {\beta_2}}, \frac{d}{{\beta_2}}(\mathcal{R}_2 - 1), C^*, R^*)$ corresponding to the projections of $D_0$ and $D_2$, respectively. The following theorem shows that all solutions of the system (4.8) with initial conditions in $\Omega_{2_p} - \{(S, I, C, R) \in \mathbb{R}^4; I = 0 \}$ approach the equilibrium $D_{2_p}$ when $\mathcal{R}_2 > 1$.

Theorem 6. If $\mathcal{R}_2 > 1$ then the Endemic equilibrium $D_{2_p}$ of the system (4.8) is globally stable in $\Omega_{2_p}$.

Proof. Proof is straightforward from the classical SIR Model. Since the first two equations do not depend on C and R, when $\mathcal{R}_0 > 1$, the solutions converge to the endemic equilibrium. And, once $I$ converges to $I^*$, $C$ converges to $C^*$.

The local asymptotic stability of $D_2$, the lemma and the theorem above, demonstrate the global stability of $D_2$ in $\Omega$ under the conditions $\mathcal{R}_2 > 1$ and $\mathcal{R}_{Inv} < 1$.

5.   Numerical results

In this paper we shall use a numerical approach to obtain information on the local stability of the CEE. The numerical values for the parameters are shown in the Table 1.

It is also necessary to choose a function that can represent the immunity period. This function was chosen to describe the biological feature of the temporary cross-protection period as described in ^[2]. For this initial numerical analysis, we choose to work with the cubic polynomial,

in order to not work with infinity time, that describes the decreasing and loss of cross protection, considering the average time being $\frac{1}{2}$ year ^[2]. In addition, we assume that after one year cross immunity is no longer effective and after this period the individual is susceptible again to other strains ^[2].

5.1.   Stability of the Coexistence Endemic equilibrium

As we saw in Section 3, the CEE $D_3$ exists only if $\mathcal{R}_0 > 1$ and, if $\mathcal{R}_{Inv} > 1$. In this case, the BE loses stability, and the CEE, with the coexistence of the two strains, rises and remains in the invariant region.

Figures 1a and 1b show the parameter region for the stability of the BE and, for the DFE and the parameter region for the existence of the CEE. Note that, as the value $\phi$ increases, the parameter region of the stability for the BE decreases, forcing the CEE to coexist within the region. Additionally, there is a threshold for the value of $\phi$, in each case, that satisfies $\mathcal{R}_{Inv} > 1$.

The characteristic equation of the system is a transcendental equation, typically having infinitely many roots. Thus, we shall make a numerical analysis in order to study the stability of the CEE, analyzing the roots of the characteristic equation, numerically, for some values of the parameter $\phi$, using it as a bifurcation parameter.

We separated in two cases:

Case (i): $\mathcal{R}_0 > 1$, $\mathcal{R}_1 < 1$

In this case, with values for the parameters in table (1) chosen $\beta_1 = 45$ and $\beta_2 = 180$, which gives $\mathcal{R}_1 = 0.87$, and $\mathcal{R}_0 = 3.46$. The $\mathcal{R}_{Inv}$ is bigger than one for $\phi = 1.23$. Then, for all value of $\phi > 1.23$ we have the existence of the CEE, and the corresponding eigenvalues, as show in Figures 2a to 2f.

Figures 2a to 2f show that the characteristic equation has a pair of conjugated complex roots, that change the sign of the real part as $\phi$ increases. Thus, a Hopf bifurcation occurs when the parameter $\phi \approx 2.19$.

Case (ii): $\mathcal{R}_0 > 1$, $\mathcal{R}_1 > 1$

In this case, with values for the parameters in Table 1 chosen $\beta_1 = 120$ and $\beta_2 = 180$, which gives $\mathcal{R}_1 = 2.31$, and $\mathcal{R}_0 = 3.46$. The $\mathcal{R}_{Inv}$ is bigger than one for $\phi = 0.205$. Then, for all value of $\phi > 0.205$ we have the existence of the CEE, and the corresponding eigenvalues, as show in Figures 3a to 3f.

Figures 3a to 3f show that the characteristic equation has a pair of conjugated complex roots, that change the sign of the real part as $\phi$ increases. Thus, a Hopf bifurcation occurs when the parameter $\phi \approx 0.244$.

Furthermore, after the Hopf Bifurcation, the CEE will be unstable, leading to a complex dynamic where sometimes the strain one resists resulting in the other becoming almost extinct, and sometimes the the opposite. Models that carry these mechanisms that can cause the coexistence of pathogens, and therefore this switch between the strains has an important role since, according to VinodKumar ^[26], different serotypes have been co-circulating in the same area, with one of them being dominant during an outbreak.

5.1.1.   Bifurcation Structure

As seen in the figures in Section 5.1, the CEE change the stability as the parameter $\phi$ changes. As $\phi$ increase from a small value to the critical value $\phi_c$, the steady state changes from a stable focus to an unstable steady state. Therefore, Hopf bifurcation occurs, and thus we conclude that closed periodic orbit will be found in a small neighborhood of $\phi_c$.

The Hopf bifurcation occurs at $\phi_c = 2.19$ (Figure 4a) and at $\phi_c = 0.244$ (Figure 4b) and, therefore, 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 occurred.

This change of stability, and thus, the Hopf bifurcation, is local. Therefore, the Hopf bifurcation does not specify what happens when the parameter is further beyond the vicinity of its critical bifurcation value ^[45,46]. Because of this, solutions will be plotted, in the next section, for different values of $\phi$, to show the asymptotic behavior also for parameter values further from the bifurcation value.

5.1.2.   Solutions of the system

Case (i): $\mathcal{R}_0 > 1$, $\mathcal{R}_1 < 1$

Figures 5a to 5e show the solutions of the system for the case $\mathcal{R}_1 < 1$ and $\mathcal{R}_0 > 1$. Figure (5a) for $\phi = 0.5$ and $\mathcal{R}_{Inv} < 1$. Thus, as $\mathcal{R}_0 > 1$ the solution converges to the value of BE.

Figures 5b and 5c for $\phi = 1.4$ and $\phi = 2.1$, respectively, where $\mathcal{R}_{Inv} > 1$. Since $\phi$ is smaller than the critical value of HB, the solution thus converges to the CEE, that is, the CEE is stable.

Figure 5d for $\phi = 2.21$ and $\mathcal{R}_{Inv} > 1$. Since the value of $\phi$ is closed to the HB critical value, periodic solutions are seen. Figure 5e for $\phi = 3.5$ where the parameter is bigger than the HB critical value. Thus, for long term behavior the solutions seen not to converge to periodic orbit, showing complex dynamics.

Case (ii): $\mathcal{R}_0 > 1$, $\mathcal{R}_1 > 1$

Figures 6a to 6e show the solutions of the system for the case $\mathcal{R}_1 > 1$ and $\mathcal{R}_0 > 1$. In Figure 6a for $\phi = 0.1$, the $\mathcal{R}_{Inv} < 1$. Thus, as $\mathcal{R}_0 > 1$ the solution converges to the value of BE. In Figure 6b for $\phi = 0.23$, smaller than the HB critical value, thus the solution converges to the CEE, that is, the CEE is stable.

Figure 6c shows the solution for $\mathcal{R}_{Inv} > 1$ and $\phi = 0.26$, closed to the HB critical value, thus periodic solutions are seen. In Figure 6d for $\phi = 0.8$, the $\mathcal{R}_{Inv} > 1$, then the parameter is bigger than the HB critical value. Thus, for long term behavior the solutions seem to converge to periodic orbit. Figure 6e for $\phi = 2$, bigger than the HB critical value. Also, in this case solutions converge to periodic orbits.

Although we conclude numerically that the solutions of the system, for values of the parameter far from the bifurcation critical value, reach an equilibrium or to a periodic orbit, in some cases the solutions seem not to converge, showing complex dynamics.

Therefore, biologically speaking, since the value for the ADE parameter is unknown, it is hard to give the smallest interval that is a representative value for ADE. This means that it is hard to predict the next episode of the disease since we may have different scenarios: coexistence of infections, periodic outbreaks, or even complex dynamics, depending on the parameter interval value for ADE.

6.   Discussion and conclusions

In this study, we have developed a mathematical model that can be applied to a multi strain infectious disease such as dengue fever, which is endemic in more than 100 countries ^[1] and remains a major public health problem.

Over the years, some mathematical models for infectious diseases, especially for dengue fever, have been developed, most of which are in the format of ODEs. However, the propagation of epidemics is not instantaneous, and it is more appropriate to model epidemics incorporating continuously distributed delay and not discrete delay or constant time.

We proposed a model with a time delay in temporary immunity, allowing general periods of immunity. The time delay was used to model phenomena so that an individual may not be immediately susceptible after recovering. In addition, a constant parameter describes the ADE effect.

A special form for the temporary immunity function was considered in Section 2. By considering the temporary immunity time as exponentially distributed, we are able to reduce the IDE system to an ODE system. The qualitative analysis of the ODE system thus gave us a visual picture of the dynamic behavior.

In Section 3, using results from IDE systems, we were able to prove the existence of four equilibria and their stability, proving that the disease will die out if $\mathcal{R}_0 < 1$. If $\mathcal{R}_{Inv} < 1$, one strain will die out, and the other will persist. And, in the last scenario, the two strains will coexist if $\mathcal{R}_{Inv} > 1$.

Lyapunov functionals can be an effective tool to prove the global stability of IDE systems. Nevertheless, we were able to prove the global stability of DFE and of the BE when $\mathcal{R}_{Inv} < 1$. This result means, biologically, that there is a range of values for the infection force so that strain one cannot invade the population if strain two is endemic. In this way, infection two may protect the population from infection one.

When comparing the particular case (exponential distributed function) with the general case, the qualitative behavior of the system was not altered by the distributed delay. However, the invasion number depends on the average cross-protection time, which depends on the choice of the function, and this can affect whether the infection can coexist.

The ADE effect was studied by considering it as a bifurcation parameter. Thus, the possibility of Hopf bifurcation was examined through a numerical analysis. This kind of bifurcation is local, and therefore, only for values close to the bifurcation critical point, can we conclude that a limit cycle exists, and that the solutions will show periodic behavior. This means that, in the scenario where the strains coexist, they can either coexist in an equilibrium or, for a specific value, they can show periodic behavior, according to the value of the ADE parameter.

For solutions of the system, for values further from the bifurcation critical value, numerical results show that the solutions either reach an equilibrium or go into a periodic orbit. However, it will be necessary to further analyze other kinds of bifurcations.

Since ADE is still an unknown factor, we can only suggest scenarios such as periodic solutions or coexistence of strains, as well as complex dynamics, depend on its value. However, a possible outbreak will be hard to predict.

Lastly, mathematical models focus on understanding the spread of infectious diseases as well as suggesting interventions to control or even predict the consequences of their propagation. With this model, we tried to understand the dynamics of different strains at the population level and explain why it is so hard to predict dengue fever. We concluded that there are certain mechanisms and intrinsic characteristics of dengue fever, such as ADE and cross protection, which may hinder the prediction of the next outbreaks of the disease.

Acknowledgments

This research was supported by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES) - Finance Code 001 and LIAM - Laboratory for Industrial and Applied Mathematics, Department of Mathematics and Statistics, York University-CA. The first author would like to thank Fernando Valdes Ravelo for recommendations in the Matlab Program, which helped to improve compilation and results.

Conflict of interest

The authors declare that they have no conflict of interest.

Supplementary material

A particular case: Exponential immunity

Here we are going to analyze qualitatively the ODE system, a particular case of the model studied. If we assume the length of immunity being exponentially distributed, then the fraction of temporarily immune individuals who remain in the $C_i$ class is $P^i(t) = e^{- \omega_i t}$, with $\omega_i > 0$, for $i = 1, 2$, then the IDE system became an ODE system. Thereby, the system already normalized, where the variables represent the fractions of the populations, can be described as follows:

Clearly, the system (6.1) always has a Disease-free equilibrium (DFE), namely,

In the case of the extinction of one of the strains, we are able to find the Boundary equilibrium (BE) of the system:

And, the BE

For biological reasons, we are looking for steady states that belong to $\Omega$, a positively invariant region, where $\Omega = \{ (S, I_1, I_2, C_1, C_2, R_1, R_2, I_{12}, I_{21}, R) \in \mathbb{R}^{10}_{+} \mbox{ such that} \ S + I_1 + I_2 + C_1 + C_2 + R_1 + R_2 + I_{12} + I_{21} + R \leq 1 \}$. In this way, the BE, $E_i$, $i = 1, 2$ is in $\Omega$ as long as ${\frac{ \beta_i}{d + \gamma} > 1}$.

In the case of the coexistence of the two strains, we are able to find the Coexistence Endemic equilibrium (CEE) that is given by

and, $S^*$ is the root of the cubic polynomial $O(S) = b_3S^{3} + b_2S^{2} + b_1 S + b_0$ where

Moreover, $S^*$ has to satisfy that $S^* < \frac{d+ \gamma}{\beta_i}$, $i = 1, 2$, in order to have the CEE, $E_3$, in the $\Omega$ region. Otherwise, if $S^*$ does not satisfy this inequality, the variables that represent the recovered and infected populations will be negative.

These discussions support and demonstrate the following theorems.

Theorem 7. If $\frac{ \beta_1}{d + \gamma} > 1$ then the system of equations (6.1), always has a BE, $E_1$, given by (6.2) in $\Omega$.And, if $\frac{\beta_2}{d + \gamma} > 1$ then the system of equations (6.1), always has a BE, $E_2$, given by (6.3) in $\Omega$.

Theorem 8. Without loss of generality, suppose that $\beta_2 > \beta_1$. If max$\{\frac{\beta_1}{d + \gamma}, \ \frac{\beta_2}{d + \gamma}\} > 1$ and,

then, the system (6.1) admits a CEE in $\Omega$ with the coexistence of the two strains.

Proof. The independent term, $b_0$, of the polynomial $O(S)$ is always positive. Since the equilibrium is given by (6.4) with $S^*$ being a root of the polynomial $O$, and this equilibrium will be in the region $\Omega$ if $S^* < \frac{d + \gamma}{\beta_i}$, for $i = 1, 2$ let us define

Then, if $\frac{{\beta_2}}{d + \gamma} > 1$ and $\mathcal{R}_{Inv} > 1$, the polynomial $O$ at $S_{min}$ is

This shows that we have a root $S^*$ of the polynomial $O$, such that, $0 < S^* < {S_{min}}$. Therefore, for this $S^*$ we have a positive equilibrium of the system (6.1), in $\Omega$, with the coexistence of the two strains.

Basic Reproduction number

We defined the Basic Reproduction number in the subsection 2.1.1 of this study. We defined

as the Basic Reproduction n for strain one. And,

as the Basic Reproduction number for strain two.

In addition, the system's overall reproduction number was defined as $\mathcal{R}_0 = \mbox{max} \{ \mathcal{R}_1, \mathcal{R}_2\}$.

Another important threshold value, the Invasion Reproduction number, was defined in that subsection, and it is given by

Remark: For our purpose and without loss of generality, we are assuming that the infection has different forces and, therefore, we assume for now on that $\beta_2 > \beta_1$. Thus, $\mathcal{R}_0 = \mathcal{R}_2$.

Stability Analysis

The local stability of the equilibria will be determined by the classical method of determining the stability of the steady states of some ODE systems, by analysis of the eigenvalues of the Jacobian matrix of the system at each equilibrium.

Theorem 9. If $\mathcal{R}_0$ then the DFE, $E_0$, of the system (6.1) is locally asymptotically stable. And, $E_0$ is unstable if $\mathcal{R}_0 > 1$.

Proof. The eigenvalues of Jacobian matrix of the system evaluated at the DFE ($E_0$) are

Clearly if $\mathcal{R}_0 < 1$ thus all the eigenvalues will be negative. If $\mathcal{R}_0 > 1$ thus at least one eigenvalue will be positive. Which proves the theorem.

Theorem 10. The BE, $E_1$, of the system (6.1) is always unstable in $\Omega$ region. And the BE, $E_2$, is stable in $\Omega$ if $\mathcal{R}_{Inv} < 1$, and unstable in $\Omega$, if $\mathcal{R}_{Inv} > 1$.

Proof. The eigenvalues of the Jacobian matrix of the system evaluated at $E_1$ are

Since $\beta_2 > \beta_1$, the eigenvalue $\lambda_9$ will be always positive. In addition, the eigenvalues of the Jacobian matrix of the system evaluated at $E_2$ are

If $\mathcal{R}_{Inv} < 1$ thus, the eigenvalue $\lambda_{9} < 0$. Furthermore, the other eigenvalue has a negative real part. Therefore, the BE, $E_2$, will be stable.

If $\mathcal{R}_{Inv} > 1$ thus, the eigenvalue $\lambda_{9} > 0$. This proves the theorem.

Remark: We can rewrite the threshold $\mathcal{R}_{Inv}$ as

Thus, if ${\frac{\alpha_1 \phi}{\beta_1} \leq 1}$, we have that $\mathcal{R}_{Inv} < \mathcal{R}_1$. This result means that there is a range of values for $\beta_1$ for which strain one cannot invade the population if the strain two is endemic. In this way, strain two may protect the population from strain one.

The analysis of the local stability of the CEE using the classical theory was not successful. We will have to deal with a characteristic polynomial of a $9 \times 9$ matrix and, with the fact that it was not possible to describe the value of $S^*$ in terms of the parameters. Thus, in order to study the stability, we will analyze it numerically, in Section 6.

Numerical analysis

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

Stability of the Coexistence Endemic equilibrium

In this section we are going to explore the stability of the CEE numerically. As we saw in section (2.1) the CEE, $E_3$, exists only if $\mathcal{R}_{0} > 1$ and $\mathcal{R}_{Inv} > 1$. In this case, the BE loses stability, and the CEE, with the coexistence of the two strains, rises and remains in the invariant region.

Figures 7a and 7b show the parameter region for the stability of the BE and, for the DFE and, the parameter region for the existence of the CEE. Note that, as the value $\phi$ increases, the parameter region of the stability for the BE decreases, forcing the CEE to coexist within the region.

To show the stability of the CEE numerically, we need to guarantee that the values of the parameters satisfy $\mathcal{R}_{Inv} > 1$. $\phi$, the parameter used to describe the ADE effect, is a parameter that its value is unknown. As we saw in Figures 7a and 7b there is a threshold for the value of $\phi$, in each case, that satisfies $\mathcal{R}_{Inv} > 1$. Starting from this threshold value, the CEE will be in the positive region and, therefore, we can look for the eigenvalues of the Jacobian Matrix at CEE. Note that the CEE can be in the invariant region even if the strain one is not established. Thus, we separated in two cases:

Case (i): $\mathcal{R}_0 > 1$, $\mathcal{R}_1 < 1$

In this case, with values for the parameters in Table 2, chosen $\beta_1 = 45$ and $\beta_2 = 180$, which gives $\mathcal{R}_1 = 0.87$, and $\mathcal{R}_0 = 3.46$. The $\mathcal{R}_{Inv}$ is bigger than one for $\phi$ bigger than $1.23$. Then, for all values of $\phi > 1.23$ we have the CEE and, the corresponding eigenvalues, as shown in Figures 8a to 8f.

Figures 8a to 8f give the stability of the CEE, $E_3$, for each value of $\phi$. We can see that the matrix has two conjugated complex roots that change the sign of the real part as $\phi$ increases. Thus, a Hopf bifurcation occurs when the parameters $\phi$ $\approx 2.52$ (Figure 8d).

Case (ii): $\mathcal{R}_0 > 1$, $\mathcal{R}_1 > 1$

In this case, with the values for the parameters in Table 2 chosen $\beta_1 = 120$ and $\beta_2 = 180$, which gives $\mathcal{R}_1 = 2.31$, and $\mathcal{R}_0 = 3.46$. The $\mathcal{R}_{Inv}$ is bigger than one for $\phi$ bigger than $0.21$. Then, for all value of $\phi > 0.21$ we have the CEE and, the corresponding eigenvalues, as shows in Figures 9a to 9f.

Figures 9a to 9f give the stability of the CEE, $E_3$ for each value of $\phi$. We can see that the matrix has two conjugated complex roots, which change the sign of the real part as $\phi$ increases. Thus, a Hopf bifurcation occurs when parameters $\phi \approx 0.244$ (Figure 9c).

Bifurcation Structure

As seen in the figures in Section 6, as $\phi$ increase from small value trough critical value, $\phi_c$, the steady state changes from a stable focus to an unstable steady state. Therefore, Hopf bifurcation occurs and thus, we conclude that closed periodic orbit will be found in a small neighbourhood of $\phi_c$.

In order to see the limit cycle around the equilibrium, in a small vicinity of the critical value, the bifurcation diagrams are shown in Figures 10a and 10b.

The Hopf bifurcation occurs at $\phi_c = 2.52$ (Figure 10a) and at $\phi_c = 0.244$ (Figure 10b) and, therefore, the solutions exhibit a small amplitude limit cycle around the CEE. A stable limit cycle arises close to the critical bifurcation point and moves away from the unstable equilibrium. Thus, it is possible to conclude that a supercritical Hopf bifurcation occurred.

This change in the stability, and thus, the Hopf bifurcation is local. Therefore, the Hopf bifurcation does not specify what happens when the parameter is further beyond the vicinity of its critical bifurcation value ^[45,46,47]. Because of this, solutions will be plotted, in the next section, for different values of $\phi$, to also show the asymptotic behaviour for parameter values further from the bifurcation value.

Solutions of the system: a particular case

In this section, we shall numerically explore the solutions of the system in order to obtain information about the dynamics and also for values of $\phi$ further from the bifurcation value.

Case (i): $\mathcal{R}_0 > 1$, $\mathcal{R}_1 < 1$

Figures 11a to 11e show the solutions of the system for the case $\mathcal{R}_1 < 1$ and $\mathcal{R}_0 > 1$. Figure 11a for $\phi = 0.5$ and $\mathcal{R}_{Inv} < 1$. Thus, as $\mathcal{R}_0 > 1$ the solution converges to the value of BE.

Figures 11b and 11c the $\mathcal{R}_{Inv} > 1$, for $\phi = 1.4$ and $\phi = 2.1$ respectively. Since the value of $\phi$ is smaller than the value of $\phi_c$, thus the solutions converge to the CEE, this is, the CEE is stable. Figure 11d the $\mathcal{R}_{Inv} > 1$, for $\phi = 2.55$ close to $\phi_c$. Thus, periodic solutions are seen.

Figure 11e the $\mathcal{R}_{Inv} > 1$, for $\phi = 3.5$ bigger than the value of $\phi_c$. For long term behavior the solutions seem not to converge to periodic orbit, showing complex dynamic.

Case (ii): $\mathcal{R}_0 > 1$, $\mathcal{R}_1 > 1$

Figures 12a to 12e show the solutions of the system for the case $\mathcal{R}_1 > 1$ and $\mathcal{R}_0 > 1$. Figure 12a for $\phi = 0.1$ the $\mathcal{R}_{Inv} < 1$. Thus, as $\mathcal{R}_0 > 1$ the solution converges to the value of BE.

Figure 12b the $\mathcal{R}_{Inv} > 1$, for $\phi = 0.23$, smaller than the value of $\phi_c$, thus the solution converges to the CEE, this is, the CEE is stable. Figure 12c the $\mathcal{R}_{Inv} > 1$, for $\phi = 0.26$ close to $\phi_c$. Thus, periodic solutions are seen.

Figure 12d the $\mathcal{R}_{Inv} > 1$, for $\phi = 0.8$ bigger than the value of $\phi_c$. For long term behaviour the solutions seen to converge to periodic orbit. Figure 12e the $\mathcal{R}_{Inv} > 1$, for $\phi = 2$ bigger than the value of $\phi_c$. Also, in this case, solutions converge to periodic orbits.

Although we conclude numerically that the solutions of the system reach an equilibrium or a periodic orbit, for values further from the bifurcation value, in some cases solutions do not converge, showing complex dynamics. Thus, analysis of other kinds of bifurcations would be necessary.

Therefore, biologically speaking, since the value of the ADE parameter is unknown, it is hard to give the smallest interval that is a representative value for ADE. This means that it is hard to predict the next episode of the disease, since we can have different scenarios: coexistence of infections, periodic outbreaks, or even complex dynamics, depending on the parameter interval value for ADE.

Definition Lyapunov Functional

Let $\Omega \subset D$ be a compact set, $L : D \rightarrow \mathbb{R}$ be a continuously differentiable function such that:

$\cdot$ $L(x) = 0$ if and only if $x = 0$;

$\cdot$ $V(x) > 0$ $\forall x \in \Omega$, $x \neq 0$;

$\cdot$ $\dot{L}(x) = \nabla L \cdot x' \leq 0$, $\forall x \in \Omega$.

Then, $L$ is called a Lyapunov functional ^[36,48,49].