# Multi-serotype nested immuno-epidemiological model for dengue hemorrhagic fever involving backward bifurcation and serotype invasion

Mostafa Adimy, Charlotte Dugourd-Camus, Ruben Taieb

► To cite this version:

Mostafa Adimy, Charlotte Dugourd-Camus, Ruben Taieb. Multi-serotype nested immuno-epidemiological model for dengue hemorrhagic fever involving backward bifurcation and serotype invasion. 2025. hal-05289175

HAL Id: hal-05289175

https://inria.hal.science/hal-05289175v1

Preprint submitted on 30 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Multi-serotype nested immuno-epidemiological model for
dengue hemorrhagic fever involving backward bifurcation
and serotype invasion

Mostafa Adimy$^a$, Charlotte Dugourd-Camus$^{a,*}$, Ruben Taieb$^{b,a}$

$^a$Inria, CNRS, Centrale de Lyon, INSA Lyon, Universite Claude Bernard Lyon 1,
Universite Jean Monnet, ICJ UMR5208, Villeurbanne, France

$^b$Université Le Havre Normandie, Normandie Univ., LMAH UR 3821, Le Havre, France

Abstract

Reinfection with the same dengue serotype is generally benign, as individuals develop protective immunity. On the other hand, in the case of reinfection with a different serotype, pre-existing antibodies can increase the risk of developing Dengue Hemorrhagic Fever (DHF), by inducing Antibody-Dependent Enhancement (ADE). To model this dynamic, we introduce a multi-scale immuno-epidemiological system. The immunological part is described by a system of ODEs representing the interaction between two antibodies (from previous and current infection) and the virus. The epidemiological part is represented by an infection-age structured SIRS system (for both the primary and secondary infections) and a recovery-age structured equation (for the first infection). A detailed mathematical analysis of the equilibrium points of the multi-scale re-infection model, including disease-free, mono-endemic and bi-endemic states, is performed. We establish necessary and sufficient conditions for the existence of backward bifurcations and derive an expression for the invasion reproduction number, which shows that the second serotype can invade the population after a mono-endemic first serotype. We also investigate the dependence of the basic and invasion reproduction numbers on the immunological parameters of the first and second infections. This gives us a better understanding of the relationship between DHF and ADE during secondary infection.

Keywords: Epidemiological Model, Immunological Model, Multi-scale
Modelling, Age-Structured PDEs, Backward bifurcation, Invasion, Dengue
Hemorrhagic Fever, Antibody-Dependent Enhancement

<footnote>*Corresponding author.</footnote>

<footnote>Email addresses: mostafa.adimy@inria.fr (Mostafa Adimy),
charlotte.dugourd@inria.fr (Charlotte Dugourd-Camus), ruben.taieb@inria.fr (Ruben
Taieb)</footnote>
# 1. Introduction

Dengue is a mosquito-borne disease caused by four viral serotypes (DENV-1, DENV-2, DENV-3, and DENV-4), none of which permanently dominate the others. Primary infection with one serotype confers lifelong immunity to that serotype but not to the others. However, a secondary infection with a different serotype can increase the risk of severe disease manifestations, such as Dengue Hemorrhagic Fever (DHF). Identifying the circulating serotype is crucial for predicting severe cases. DHF are largely driven by Antibody-Dependent Enhancement (ADE), a phenomenon in which non-neutralizing antibodies from a previous infection facilitate rather than block subsequent infections [19].

Numerous within-host models have been developed to study the interactions between dengue and the immune system, as well as to elucidate the relationship between DHF and ADE. For instance, [7] uses patient data to model the role of antibodies in controlling dengue infection, offering insights for the development of vaccines and antivirals. Similarly, [16] and [20] examine how cross-reactive antibodies amplify viral replication and discuss partial success in attenuating ADE by modifying antibody structures. The study [26] provides empirical evidence of 1-3 years of cross-protection following dengue infection, which has implications for vaccine development. Additionally, [27] shows that DENV-1 infections generally result in higher viral loads than DENV-2 or DENV-3, with severe symptoms correlating to early viral levels and faster clearance in previously infected individuals. Studies [2, 6] simulate ADE by modeling the competition between neutralizing and enhancing antibodies, thus predicting the dynamics and severity of the disease. The study [1] focuses on ADE, where antibodies of an earlier serotype increase the risk of infection, thus contributing to DHF. It uses an epidemic model to examine the impact of ADE on dengue susceptibility, transmission, and mortality. It finds that while increased transmission and mortality can cause periodicity or chaos, only increased transmission allows closely related serotypes to coexist. In [8], a model of dengue dynamics shows how variations in immune responses can explain differences in viral dynamics between primary and secondary infections, supporting ADE with enhanced viral entry during secondary infections. Using a system of ODEs, the paper [5] investigates the dynamics induced by the ADE phenomenon in multi-serotype diseases. To increase secondary infectivity in the presence of immunity resulting from previous infection with a different serotype, the viral growth rate is increased, giving rise to the ADE phenomenon. This phenomenon is associated with the emergence of oscillations via a Hopf bifurcation. The paper derives approximations of the ADE parameter required to induce these oscillations and analyzes the associated bifurcations. It then examines the stability of these dynamics by introducing stochastic perturbations. It also presents a preliminary analysis of the effect of a single serotype vaccination. Our recent study [10] extends the approaches from [2, 6] by developing a model to analyze viral activity in the presence of both neutralizing and enhancing antibodies, helping to identify their roles in disease progression and ADE during secondary infections.

The nested modeling approach, first introduced by Gilchrist and Sasaki in
2002 [14], was inspired by the work of Halstead [19], who introduced the concept of ADE in dengue infections. Halstead's study laid the foundation for the development of nested immuno-epidemiological models, which integrate immunological and epidemiological scales into a unified framework. This approach improves our understanding of the dynamics of an epidemic by explicitly considering the host immune status and its influence on the disease spread. Several studies have employed this multi-scale approach. For example, [17] investigates the complex interactions between dengue virus serotypes and host immunity, particularly the role of cross-reactive antibodies in enhancing infection severity and contributing to DHF. A similar approach is used in [18], which develops a nested immuno-epidemiological model for HIV transmission on complex networks, combining ordinary differential equations (ODEs) with partial differential equations (PDEs). The study [24] links a within-host model of Mycobacterium Avium ssp. Paratuberculosis (MAP) with an epidemiological framework to better understand Johne's disease in dairy cattle.

Nested immuno-epidemiological models are useful tools for studying the interaction between host immunity, pathogen dynamics and epidemic spread, particularly in multi-serotype (and multi-strain) diseases. These models include immunological factors such as immune decay and cross-protection, as well as epidemiological factors such as infection and recovery rates, leading to a more precise analysis of disease transmission. These models are particularly adapted to the study of reinfections, in particular to understand the relationship between DHF and ADE.

In this paper, we investigate a multi-scale immuno-epidemiological system for dengue reinfection. The immunological level is described by a system of ODEs representing the interaction between two antibodies, one from a previous infection and the other from the current infection, and the virus. The epidemiological component is represented by an SIRS system structured by time-since-last-infection (both primary and secondary) and time-since-last-recovery from the primary infection. The transition rates between these classes may depend on the immunological state of the host, and the initial viral load of an infected individual may be influenced by the viral load of the individual who transmits the disease. We sometimes consider that the viral load at the beginning of the infection is independent of that transmitted by the infected individual. Indeed, as we are focusing on dengue, which is a mosquito-borne disease, this hypothesis is consistent. On the other hand, as dengue reinfections beyond two are infrequent, we mainly consider the first and second infections. Nevertheless, the proposed model can be extended to the case of transmission of the viral load according to a relatively general probability law and also to multiple consecutive or simultaneous infections, and thus it could be extended to other diseases involving multi-serotypes (and multi-strains).

We analyze the existence and stability of the disease-free, mono-endemic and bi-endemic equilibrium of the multi-scale reinfection model. We establish a necessary and sufficient condition for the existence of a backward bifurcation. Then, we obtain an expression for the invasion reproduction number, which determines whether the second serotype can invade a population already endemic
to the first serotype. We also study how the basic and invasion reproduction numbers depend on the immunological parameters associated with the first and second infection. In particular, we examine the role of the immune response - such as cross-protection and ADE - in the transmission dynamics of two successive infections. Finally, we complete our theoretical analysis with numerical simulations to highlight the potential of the proposed model. In particular, we use the backward bifurcation and the invasion reproduction number to explain the relationship between ADE phenomenon and the outbreak of the DHF epidemic during secondary infections.

Models similar to the one proposed in this paper have been studied by Gandolfi et al. [13] and Gulbudak et al. [17]. Our model extends the study of [13] by accounting for reinfection and multi-serotypes (or multi-strains). Additionally, we complement the epidemiological analysis in [17] by explicitly identifying and proving a necessary and sufficient condition for backward bifurcation, and obtaining an explicit expression for the invasion reproduction number, which provides a condition for the two serotypes to coexist in the population. All results in this paper are proven within the age-structured immuno-epidemiological framework. Furthermore, unlike [17], we take into account the difference in within-host dynamics between primary and secondary infections, which is an important consideration for dengue, where secondary infections significantly increase the risk of DHF due to ADE phenomenon. The study [17] also assumes that an individual develops DHF when their viral load exceeds a certain threshold, which is a quantification of DHF cases that we did not use. Finally, unlike in [17], we only implicitly model the dynamics of vector populations, considering the large differences in the time scale between the dynamics of the host and that of the vector.

## 2. Within-host model

In this section, we focus on the immunological scale and present two within-host models. First, we consider a primary infection, followed by a secondary infection in a separate model. A more detailed study of secondary infection modeling, particularly in the case of dengue, can be found in [10]. For the immune response, we consider only those antibodies that are neutralizing during the primary infection. These antibodies may either remain neutralizing (more strongly than during the first infection) or become enhancing, as observed in dengue, during the secondary infection, which could involve a similar or different viral serotype.

### 2.1. Immunological model of first infection

Let us consider a primary infection by a virus. The within-host model used is as follows:

$$ 
\left\{ 
\begin{aligned}
\frac{dv}{d\tau} &= r(a)v\left(1 - \frac{v}{K}\right) - \gamma_v v &:= F(v, a), \\
\frac{da}{d\tau} &= \eta(v) - \gamma_a a &:= H(v, a),
\end{aligned}
\right. 
\qquad (1) $$
where $\tau \mapsto v(\tau)$ represents the virus and $\tau \mapsto a(\tau)$ represents the antibody concentration during the time of infection $\tau$. The virus is neutralized by the antibody due to the growth rate $r$ which is a positive and decreasing function. An example is

$$r(a) = c \frac{1 + \sum_{i=1}^{n} \xi_i C_n^i a^i}{1 + \sum_{i=1}^{n} C_n^i a^i},$$

where $c > 0$, $\xi_i = 1$, for $1 \le i \le n-1$ and $0 \le \xi_n < 1$, and $C_n^i$ the binomial coefficient [10]. $n$ denotes the number of antigenic epitopes that can be recognized by the antibody. Values such as $n = 2$ or $n = 3$ are sometimes introduced to account for cross-reactivity or competition between antibodies [10]. However, this parameter does not affect the qualitative dynamics of the immuno-epidemiological system. In all our simulations, we fixed $n = 3$. For additional details regarding this parameter, see [10] and the references therein. The carrying capacity of the virus is $K$. Antibody production is stimulated by the virus through the function $v \mapsto \eta(v)$, which is an increasing and positive function. The mortality rates of the virus and antibodies are given by $\gamma_v$ and $\gamma_a$, respectively. Finally, the parameters of System (1) depend on the specific virus serotype under consideration. System (1) always has a trivial equilibrium,

$$E_{01} := (0, \bar{a}) = \left(0, \frac{\eta(0)}{\gamma_a}\right).$$

The immunological reproduction number is

$$R_{01} := \frac{1}{\gamma_v} r(\bar{a}) = \frac{1}{\gamma_v} r\left(\frac{\eta(0)}{\gamma_a}\right).$$

**Theorem 2.1.** • If $R_{01} < 1$, then $E_{01}$ is the only equilibrium and it is locally asymptotically stable.

• If $R_{01} > 1$, then $E_{01}$ is unstable and there is a unique coexistence equilibrium $E_1^* := (a^*, v^*)$ that is locally asymptotically stable.

*Proof.* We have $H(v, \bar{a}) = \eta(v) - \eta(0) > 0$, for all $v > 0$. The Jacobian reads

$$J(v, a) = \begin{pmatrix} -\gamma_a & \eta'(v) \\ v\left(1 - \frac{v}{K}\right) r'(a) & r(a) - 2v\frac{r(a)}{K} - \gamma_v \end{pmatrix}.$$

At the trivial equilibrium $E_{01}$, this reads

$$J(0, \bar{a}) = \begin{pmatrix} -\gamma_a & \eta'(0) \\ 0 & r(\bar{a}) - \gamma_v \end{pmatrix}.$$

The eigenvalues are $-\gamma_a < 0$ and $r(\bar{a}) - \gamma_v = \gamma_v (R_{01} - 1)$. Thus, $E_{01}$ is locally asymptotically stable if $R_{01} < 1$ and it is unstable if $R_{01} > 1$. The existence of an equilibrium $E_1^* := (v^*, a^*)$ with $v^* > 0$, is equivalent to solve the equation

$$L_1(v^*) := \frac{1}{\gamma_v} r\left(\frac{1}{\gamma_a}\eta(v^*)\right) \left(1 - \frac{v^*}{K}\right) = 1, \text{ where } a^* = \frac{1}{\gamma_a}\eta(v^*) \text{ and } v^* > 0.$$
If a solution $v^*$ exists, it satisfies $0 < v^* < K$. Thus, we have $L_1 : [0, K] \to [0, +\infty)$. Since $r$ is decreasing and $\eta$ is increasing, $L_1$ is decreasing, with $L_1(0) = R_{01}$ and $L_1(K) = 0$. Therefore, a unique $v^* \in (0, K)$ satisfying $L_1(v^*) = 1$ exists if and only if $R_{01} > 1$. At this equilibrium $E_1^* = (v^*, a^*)$, the Jacobian reads

$$J(v^*, a^*) = \begin{pmatrix} -\gamma_a & \eta'(v^*) \\ \gamma_v v^* \frac{r'(a^*)}{r(a^*)} & -v^* \frac{r(a^*)}{K} \end{pmatrix}.$$

The corresponding characteristic equation is

$$\lambda^2 - \mathrm{Tr}(J(v^*, a^*))\lambda + \det(J(v^*, a^*)) = 0.$$

According to the Routh-Hurwitz criteria, $E_1^*$ is locally asymptotically stable if and only if $\mathrm{Tr}(J(v^*, a^*)) < 0$ and $\det(J(v^*, a^*)) > 0$. In fact, since $\eta'(v^*) > 0$ and $r'(a^*) < 0$, we have

$$\mathrm{Tr}(J(v^*, a^*)) = -\gamma_a - v^* \frac{r(a^*)}{K} < 0$$

$$\text{and } \det(J(v^*, a^*)) = \left( \gamma_a \frac{r(a^*)}{K} - \eta'(v^*) \gamma_v \frac{r'(a^*)}{r(a^*)} \right) v^* > 0.$$

Thus, $E_1^*$ is locally asymptotically stable. $\square$

**Theorem 2.2.** • The set $[0, K] \times [\bar{a}, \eta(K)/\gamma_a]$ is positively invariant and globally attractive for System (1). In addition, the two equilibrium points $E_{01}$ and $E_1^*$ belong to $[0, K] \times [\bar{a}, \eta(K)/\gamma_a]$.

• If $R_{01} < 1$, the trivial equilibrium $E_{01}$ is globally asymptotically stable.

*Proof.* Let $(v, a)$ be a solution. If there exists $\tau \ge 0$ such that $v(\tau) = K$ (respectively, $v(\tau) = 0$), then $v'(\tau) < 0$ (respectively, $v'(\tau) \ge 0$). Then, $v$ cannot cross 0 or $K$ from within $[0, K]$. We conclude that $[0, K]$ is invariant for $v$. Using the same reasoning, we obtain that $[\bar{a}, \eta(K)/\gamma_a]$ is invariant for $a$. Furthermore, it is straightforward to see that $E_{01} \in [0, K] \times [\bar{a}, \eta(K)/\gamma_a]$. Moreover, when $R_{01} > 1$, the coexistence equilibrium $E_1^*$ belongs to $[0, K] \times [\bar{a}, \eta(K)/\gamma_a]$. Now, to prove the global attractiveness of $[0, K] \times [\bar{a}, \eta(K)/\gamma_a]$, we start by noting that if an initial condition $v(0) > K$, then the solution $v$ is decreasing as long as $v$ does not enter the interval $[0, K]$. With the same observation for $a$, we prove that the set $[0, K] \times [\bar{a}, \eta(K)/\gamma_a]$ is globally attractive for System (1). To prove the global asymptotic stability of $E_{01}$, we use the set $[0, K] \times [\bar{a}, \eta(K)/\gamma_a]$ as trapping region for the solution when $R_{01} < 1$. As this set contains only the equilibrium $E_{01}$ and no closed trajectories (because the flow has negative divergence), we apply the Poincaré-Bendixson theorem to conclude that $E_{01}$ is globally asymptotically stable. $\square$

Figure 1 shows the trajectory $(v(\tau), a(\tau))$ of System (1) for different initial conditions $(v(0), a(0))$ in the case $R_{01} < 1$. We observe that if the initial level of
<figure><img src="image_2.png" /><figcaption>Figure 1: The phase portrait in the case of a first infection for different initial conditions and $R_{01} < 1$, is plotted. The red star corresponds to the disease-free equilibrium $E_{01} := (0, \eta(0)/\gamma_a)$. The following parameters were used: $r(a) = 7.3 \times (1+2a+0.8a^2)/(1+2a+a^2)$, $\eta(v) = 0.2(1+v)$, $\gamma_v = 24 \log(2)/2.7$, $\gamma_a = \log(2)/40$ and $K = 10^7$. Note that the parameter values are arbitrary.</figcaption></figure>

antibodies is too low, the virus population increases. This increase stimulates antibody production, which, after reaching a certain threshold, leads to a reduction in the viral load. As the virus quantity decreases, antibody production also declines, eventually leading the system to the trivial equilibrium $E_{01} := (0, \bar{a})$ (represented by the red star in Figure 1).

## 2.2. Immunological model of second infection

Consider now a secondary-serotype viral infection following a primary-serotype infection. Two antibody species are involved, the antibody a associated with the first serotype and the antibody b associated with the second serotype. In this case, the within-host model is

$$ 
\left\{ 
\begin{aligned}
\frac{dv}{d\tau} &= r(a, b)v \left(1 - \frac{v}{K}\right) - \gamma_v v \quad := F(v, a, b), \\
\frac{da}{d\tau} &= \eta(v) - \gamma_a a \quad := H(v, a), \\
\frac{db}{d\tau} &= \bar{\eta}(v) - \gamma_b b \quad := \bar{H}(v, b),
\end{aligned}
\right. 
\quad (2) $$

where r is the growth rate of the virus. It depends on the two antibodies as follows

$$ r(a, b) = c \frac{1 + \sum_{i=1}^{n_1} \xi_i C_{n_1}^i a^i}{1 + \sum_{i=1}^{n_1} C_{n_1}^i a^i} \times \frac{1 + \sum_{j=1}^{n_2} \xi_j C_{n_2}^j b^j}{1 + \sum_{j=1}^{n_2} C_{n_2}^j b^j}, $$
with $c > 0$, $\xi_i \ge 1$ for $1 \le i \le n_1 - 1$ if the first antibodies enhance the second infection, and $\xi_i \le 1$ otherwise. Also, $\xi_j = 1$ for $1 \le j \le n_2 - 1$ and $0 \le \xi_{n_1}, \xi_{n_2} < 1$, with $C_{n_1}^i, C_{n_2}^j$ being the binomial coefficients. It corresponds to an independent binding described in [10]. One could also consider competitive and synergistic bindings. The system parameters depend on the two viral serotypes under consideration. System (2) always has a trivial equilibrium,

$$E_{02} := (0, \bar{a}, \bar{b}) = \left(0, \frac{\eta(0)}{\gamma_a}, \frac{\bar{\eta}(0)}{\gamma_b}\right).$$

We define the immunological basic reproduction number for secondary infection as

$$R_{02} := \frac{1}{\gamma_v} r(\bar{a}, \bar{b}) = \frac{1}{\gamma_v} r\left(\frac{\eta(0)}{\gamma_a}, \frac{\bar{\eta}(0)}{\gamma_b}\right). \quad (3)$$

Since the two functions $\eta$ and $\bar{\eta}$ are assumed to be increasing, we have the properties

$$H(v, \bar{a}) = \eta(v) - \eta(0) > 0 \quad \text{and} \quad \bar{H}(v, \bar{b}) = \bar{\eta}(v) - \bar{\eta}(0) > 0, \quad \text{for all } v > 0.$$

Figure 2 illustrates the phase portrait $(v(\tau), a(\tau), b(\tau))$ of System (2) for different initial conditions in the cases $R_{02} < 1$ and $R_{02} > 1$.

**Theorem 2.3.** • If $R_{02} < 1$, then $E_{02}$ is locally asymptotically stable.

• If $R_{02} > 1$, then $E_{02}$ is unstable and there exists at least one coexistence equilibrium $E_2^* = (v^*, a^*, b^*)$.

*Proof.* The Jacobian of the system reads

$$J(v, a, b) = \begin{pmatrix} r(a, b) \left(1 - 2\frac{v}{K}\right) - \gamma_v & \frac{\partial r}{\partial a}(a, b)v \left(1 - \frac{v}{K}\right) & \frac{\partial r}{\partial b}(a, b)v \left(1 - \frac{v}{K}\right) \\ \frac{\eta'(v)}{\bar{\eta}'(v)} & -\gamma_a & 0 \\ \frac{\eta'(v)}{\bar{\eta}'(v)} & 0 & -\gamma_b \end{pmatrix}.$$

At the trivial equilibrium, it becomes

$$J(0, \bar{a}, \bar{b}) = \begin{pmatrix} \gamma_v (R_{02} - 1) & 0 & 0 \\ \eta'(0) & -\gamma_a & 0 \\ \bar{\eta}'(0) & 0 & -\gamma_b \end{pmatrix}.$$

Then, the trivial equilibrium $E_{02}$ is locally asymptotically stable if $R_{02} < 1$ and it is unstable if $R_{02} > 1$. The set of coexistence equilibrium points is in one-to-one correspondence with the roots of

$$L_2(v) := r\left(\frac{\eta(v)}{\gamma_a}, \frac{\bar{\eta}(v)}{\gamma_b}\right)\left(1 - \frac{v}{K}\right) - \gamma_v = 0. \quad (4)$$

We have

$$L_2(0) := \gamma_v(R_{02} - 1) \quad \text{and} \quad L_2(K) = -\gamma_v < 0.$$

Then, if $R_{02} > 1$, there is at least one coexistence equilibrium, with $v^* \in (0, K)$. □
It is clear that depending on the function $r : (a, b) \mapsto r(a, b)$, we may have more than one coexistence equilibrium and even coexistence equilibrium points in the case $R_{02} < 1$. Now, we put

$$ \eta = (\eta, \bar{\eta}), \quad \Gamma = (\gamma_a, \gamma_b), \quad \Gamma^{\sigma} = (\gamma_b, \gamma_a), \quad \Gamma\eta = (\gamma_a\eta, \gamma_b\bar{\eta}) \quad \text{and} \quad \Gamma^{\sigma}\eta = (\gamma_b\eta, \gamma_a\bar{\eta}), $$

and we note

$$ \nabla r(a^*, b^*) := \left( \frac{\partial r}{\partial a}(a^*, b^*), \frac{\partial r}{\partial b}(a^*, b^*) \right)^T. $$

**Proposition 2.4.** Let $E_2^* = (v^*, a^*, b^*)$ be a coexistence equilibrium. Suppose that

$$ \eta'(v^*) \cdot \nabla r(a^*, b^*) \le 0, \quad \Gamma\eta'(v^*) \cdot \nabla r(a^*, b^*) \le 0 \quad \text{and} \quad \Gamma^\sigma\eta'(v^*) \cdot \nabla r(a^*, b^*) \le 0. \tag{5} $$

Then, $E_2^*$ is locally asymptotically stable.

**Remark 2.5.** If we assume that $\gamma_a = \gamma_b$, then Condition (5) becomes

$$ \eta'(v^*) \cdot \nabla r(a^*, b^*) \le 0. $$

*Proof.* The Jacobian at the coexistence equilibrium $E_2^* = (v^*, a^*, b^*)$ is

$$ J(E_2^*) = \begin{pmatrix} -r(a^*, b^*) \frac{v^*}{K} & \frac{\partial r}{\partial a}(a^*, b^*) v^* \left(1 - \frac{v^*}{K}\right) & \frac{\partial r}{\partial b}(a^*, b^*) v^* \left(1 - \frac{v^*}{K}\right) \\ \eta'(v^*) & -\gamma_a & 0 \\ \bar{\eta}'(v^*) & 0 & -\gamma_b \end{pmatrix}. $$

Then, the characteristic polynomial $\chi$ reads

$$ \begin{aligned} \chi(X) ={}& (X + \gamma_a)(X + \gamma_b) \left( X + r(a^*, b^*) \frac{v^*}{K} \right) \\ & - (X + \gamma_b)\eta'(v^*) \frac{\partial r}{\partial a}(a^*, b^*) v^* \left( 1 - \frac{v^*}{K} \right) \\ & - (X + \gamma_a)\bar{\eta}'(v^*) \frac{\partial r}{\partial b}(a^*, b^*) v^* \left( 1 - \frac{v^*}{K} \right). \end{aligned} $$

It has the form

$$ \chi(X) = p_3 X^3 + p_2 X^2 + p_1 X + p_0, $$

with

$$ p_3 = 1 > 0, $$

$$ p_2 = r(a^*, b^*) \frac{v^*}{K} + \gamma_a + \gamma_b > 0, $$

$$ p_1 = (\gamma_b + \gamma_a) r(a^*, b^*) \frac{v^*}{K} + \gamma_a \gamma_b - v^* \left(1 - \frac{v^*}{K}\right) \eta'(v^*) \cdot \nabla r(a^*, b^*), $$

$$ p_0 = \gamma_a \gamma_b r(a^*, b^*) \frac{v^*}{K} - v^* \left(1 - \frac{v^*}{K}\right) \Gamma^\sigma \eta'(v^*) \cdot \nabla r(a^*, b^*). $$
Furthermore, we have

$$p_2p_1 - p_0 = r(a^*, b^*) \left( (\gamma_a + \gamma_b)r(a^*, b^*) \frac{v^*}{K} - v^* \left(1 - \frac{v^*}{K}\right) \eta'(v^*) \cdot \nabla r(a^*, b^*) \right) \\ + (\gamma_a + \gamma_b) \left( (\gamma_a + \gamma_b)r(a^*, b^*) \frac{v^*}{K} + \gamma_a \gamma_b \right) \\ - v^* \left(1 - \frac{v^*}{K}\right) \Gamma \eta'(v^*) \cdot \nabla r(a^*, b^*).$$

Thanks to Condition (5), we conclude that

$$p_1 > 0, \quad p_2 > 0 \quad \text{and} \quad p_2p_1 - p_0 > 0.$$

By the Routh-Hurwitz criterion, $E_2^*$ is locally asymptotically stable. □

**Remark 2.6.** • If we suppose $\gamma := \gamma_a = \gamma_b$, we obtain

$$L'_2(v^*) = \frac{1}{\gamma} \eta'(v^*) \cdot \nabla r(a^*, b^*) \left(1 - \frac{v^*}{K}\right) - \frac{r(a^*, b^*)}{K}.$$

This means that

$$p_2p_1 - p_0 = r(a^*, b^*)v^*\gamma \left(\frac{r(a^*, b^*)}{K} - L'_2(v^*)\right) \\ + \gamma^2 v^* \left(3\frac{r(a^*, b^*)}{K} - L'_2(v^*)\right) + 2\gamma^3.$$

We conclude that if

$$L'_2(v^*) < \frac{r(a^*, b^*)}{K},$$

then Condition (5) is satisfied. Specifically, we know that there exists at least one solution $v^* > 0$ such that $L'_2(v^*) < 0$, since $L_2(0) > 0$ and $L_2(K) < 0$. Therefore, there exists at least one $E_2^*$ locally asymptotically stable.

* The set $[0, K] \times [0, \eta(K)/\gamma_a] \times [0, \bar{\eta}(K)/\gamma_b]$ is globally attractive and positively invariant for System (2). In addition, all the equilibrium points belong to this set. We also expect global asymptotic stability of $E_{02}$ when $R_{02} < 1$. But we do not intend to prove it.

* A third infection could be considered, and so on. Nevertheless, we limit our analysis to two infections.

### 2.3. Between two infections

As in the paper by Gulbudak and Browne [17], we assume that, between two infections, pre-existing antibodies a decrease according to a nonnegative function w. We consider that individuals who recover from their first infection
<figure><img src="image_3.png" /><figcaption>Figure 2: The phase portrait in the case of a second infection (with one enhancing and one neutralizing antibody) for different initial conditions, where $R_{02} < 1$ (left) and $R_{02} > 1$ (right). The red star corresponds to the steady state reached (left: disease-free; right: coexistence). We have taken $r(a, b) = 6.8 \times (1+2 \times 1.5a+0.95a^2)/(1+2a+a^2) \times (1+2b+0.85b^2)/(1+2b+b^2)$ (left), $r(a, b) = 7.2 \times (1+2 \times 1.5a+0.95a^2)/(1+2a+a^2) \times (1+2b+0.85b^2)/(1+2b+b^2)$ (right), $\eta(v) = \bar{\eta}(v) = 0.2(1+v)$, $\gamma_v = 24 \log(2)/2.7$, $\gamma_a = \gamma_b = \log(2)/40$ and $K = 10^7$. Note that these values are arbitrary.</figcaption></figure>

remain immune to the corresponding serotype. Therefore, the antibody level decreases only to a certain threshold $\bar{a} > 0$. The reached antibody level can be critical during a second infection by a serotype different from the first, potentially inducing the ADE phenomenon (and potentially DHF). If the antibody level is too low, it may fail to induce a cross-neutralizing reaction. Hence, the antibodies a satisfy the following equation between two infections

$$ \frac{da}{d\tau} = -w(a) \quad \text{with} \quad \lim_{a \to \bar{a}} w(a) = 0. \qquad (6) $$

Figure 3 shows the trajectory of the within-host model during two successive infections, with two neutralizing antibodies (second infection with homologous serotype) or a mix of neutralizing and enhancing antibodies (second infection with heterologous serotype) during the second infection. The ADE phenomenon occurs during the second heterologous infection, leading to a higher viral load, which may result in DHF.

### 3. Epidemiological model

We consider a first viral infection with a serotype we call m. We denote $(v(\tau, v_0), a(\tau, v_0))$ the solution at time $\tau$ of System (1), with $F = F_m$ and $H = H_m$, where F and H are defined in (1), associated with the initial condition $(v_0, \bar{a})$, $\bar{a} > 0$ fixed. We can therefore write

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial v}{\partial \tau}(\tau, v_0) = F_m(v, a), \\ \\ \displaystyle \frac{\partial a}{\partial \tau}(\tau, v_0) = H_m(v, a), \end{array} \right. \quad \text{with the initial condition} \quad \left\{ \begin{array}{l} v(0, v_0) = v_0, \\ a(0, v_0) = \bar{a}. \end{array} \right. \qquad (7) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Example of the trajectory of the within-host model during two successive infections, with two neutralizing antibodies (top) or neutralizing/enhancing antibodies (bottom) during the second infection. We have taken $\gamma_v = 24 \log(2)/2.7$, $\gamma_a = \gamma_b = \bar{w} = \log(2)/40$, $K = 10^7$, $\bar{a} = 5$, $w(a) = -\bar{w}(a - \bar{a})$, $a > \bar{a}$, $r(a) = 7.3 \times (1 + 2a + 0.8a^2)/(1 + 2a + a^2)$, $r(a, b) = 7.2 \times (1 + 2a + 0.95a^2)/(1 + 2a + a^2) \times (1 + 2b + 0.85b^2)/(1 + 2b + b^2)$ (top), $r(a, b) = 7.2 \times (1 + 2 \times 1.5a + 0.95a^2)/(1 + 2a + a^2) \times (1 + 2b + 0.85b^2)/(1 + 2b + b^2)$ (bottom), $\eta(v) = \bar{\eta}(v) = 0.2(1 + v)$.</figcaption></figure>

We choose $\bar{a}$ such that $F_m(0, \bar{a}) = H_m(0, \bar{a}) = 0$. Let $\tilde{a}(\tau, a_0)$ be the solution at time $\tau$ of System (6), with $w = w_m$, associated with initial condition $a_0 \ge \bar{a}$:

$$ \frac{\partial \tilde{a}}{\partial \tau}(\tau, a_0) = -w_m(\tilde{a}(\tau, a_0)), \quad \tilde{a}(0, a_0) = a_0. \qquad (8) $$

We choose the function $w_m$ such that the solution $\tau \mapsto \tilde{a}(\tau, a_0)$ keeps the interval $[\bar{a}, +\infty)$ invariant. For simplicity, we choose $w_m(a) = \bar{w}_m(a - \bar{a})$, $a > \bar{a}$, where $\bar{w}_m$ is a positive constant. The amount of antibodies after a certain recovery time $\tau$ is then obtained directly by

$$ \tilde{a}(\tau, a_0) = e^{-\bar{w}_m \tau} (a_0 - \bar{a}) + \bar{a}. \qquad (9) $$

Similarly, we consider a secondary viral infection with a serotype $n$. We denote by $(v(\tau, v_0, a_0), a(\tau, v_0, a_0), b(\tau, v_0, a_0))$ the solution at time $\tau$ of System
(2), with $F = F_{mn}$, $H = H_{mn}$ and $\bar{H} = \bar{H}_{mn}$, where $F$, $H$ and $\bar{H}$ are defined in (2), associated to the initial condition $(v_0, a_0, \bar{b})$, with $a_0 \ge \bar{a}$ and $\bar{b} > 0$ fixed:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial v}{\partial \tau}(\tau, v_0, a_0) = F_{mn}(v, a, b), \\[1.5em] \displaystyle \frac{\partial a}{\partial \tau}(\tau, v_0, a_0) = H_{mn}(v, a), \quad \text{with the initial condition} \\[1.5em] \displaystyle \frac{\partial b}{\partial \tau}(\tau, v_0, a_0) = \bar{H}_{mn}(v, b), \end{array} \right. \qquad (10) $$

Consider now the epidemiological scale of dengue, with its four different serotypes. For the sake of simplicity, we assume that each individual is exposed to no more than two serotypes in the course of a lifetime. Nevertheless, the model we are developing here can be extended to the case of more than two infections. We denote $S(t)$ the susceptible individuals (never infected) and $i_m(t, \tau, v_0)$ the individuals infected with serotype

$$ m \in J := \{1, \dots, 4\}. $$

The individuals $i_m(t, \tau, v_0)$ have been infected for a time $\tau$ and have an initial viral load $v_0$. They can recover and become again susceptible. They keep a certain immunity from the serotype $m$. The antibodies they developed during the first infection with the serotype $m$ will begin to decrease once they have recovered. We denote $S_m(t, \tau, a_0)$ the recovered individuals from the serotype $m$. They are susceptible to serotype $n \in J$, with respect to the time of recovery $\tau$ and initial primary antibody level $a_0$. Then, $i_{mn}(t, \tau, v_0, a_0)$ represents the number of individuals reinfected with serotype $n$, since a time $\tau$, after a primary infection with serotype $m$. Their initial viral load and primary initial antibody level at $\tau = 0$ are $v_0$ and $a_0$, respectively. Finally, $R(t)$ is the number of individuals recovered for both serotypes. We assume that all new individuals entering the population are susceptible, $\Lambda$ is the constant rate of new arrivals. The natural mortality rate is given by $\delta$. Infections can induce additional mortality, therefore we introduce $\delta_m$, $\delta_{mn} \ge \delta$, with $m, n \in J$, the modified mortality rates for individuals in $i_m$ and $i_{mn}$ compartments. The recovery rates $\tilde{\mu}_m$, $\tilde{\mu}_{mn}$ depend on the serotype of the primary and secondary infection, as well as the time since the infection, the initial viral load and antibody level. Susceptible individuals $S(t)$ are infected by any serotype and any primary or secondary infected individual with a force of infection of serotype $m$,

$$ \Gamma_m(t) := \lambda_m(t) + \sum_{l=1}^{4} \lambda_{lm}(t), \quad m \in J, $$

where the different forces of infection, associated to infectious individuals $i_m$ and $i_{mn}$, are given by the integrals

$$ \lambda_m(t) = \int_0^{+\infty} \int_0^{+\infty} \tilde{\beta}_m(\tau, v_0) i_m(t, \tau, v_0) dv_0 d\tau, \qquad (11) $$
and

$$ \lambda_{mn}(t) = \int_{0}^{+\infty} \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} \tilde{\beta}_{mn}(\tau, v_0, a_0) i_{mn}(t, \tau, v_0, a_0) da_0 dv_0 d\tau. \quad (12) $$

$\tilde{\beta}_m$ and $\tilde{\beta}_{mn}$ are the rates of infection depending on the immunological variables,

$$ \tilde{\beta}_m(\tau, v_0) = \beta_m(v(\tau, v_0), a(\tau, v_0)) $$

$$ \text{and} \quad \tilde{\beta}_{mn}(\tau, v_0, a_0) = \beta_{mn}(v(\tau, v_0, a_0), a(\tau, v_0, a_0), b(\tau, v_0, a_0)), $$

with $(v(\tau, v_0), a(\tau, v_0))$ solution of (7) and $(v(\tau, v_0, a_0), a(\tau, v_0, a_0), b(\tau, v_0, a_0))$ solution of (10). Individuals $S_m(t)$, recovered from the serotype $m$, become susceptible to any serotype $n \in J$, and are thus infected with a force

$$ \Gamma_{mn}(t, \tau, a_0) := \tilde{\alpha}_{mn}(\tau, a_0) \left( \lambda_n(t) + \sum_{l=1}^{4} \lambda_{ln}(t) \right), \quad m, n \in J, $$

where the parameter $\tilde{\alpha}_{mn}$ depends on the immunological state of individuals that recovered from the serotype $m$ and are susceptible to serotype $n$,

$$ \tilde{\alpha}_{mn}(\tau, a_0) := \alpha_{mn}(\tilde{a}(\tau, a_0)), \quad m, n \in J, $$

where $\tilde{a}(\tau, a_0)$ is the solution of (8). The parameter $\tilde{\alpha}_{mn}$ can be used to regulate reinfections. A value greater than 1 indicates that successive infections with serotype $m$ followed by serotype $n$ are favored, whereas a value less than 1 suggests the opposite. In most cases, we take $\tilde{\alpha}_{mn} < 1$ if $m = n$ because second infections with the same serotype are infrequent. We could take $\tilde{\alpha}_{mn} = 0$ if $m = n$ meaning that we are neglecting such secondary infections. To simplify, we consider $\tilde{\alpha}_{mn} = 1$ if $m \neq n$. The epidemiological model is illustrated in Figure 4 and described by System (13). There are actually 26 equations, but only 5 are written down because we have the three middle equations depending on $m, n \in J := \{1, \dots, 4\}$:

$$ \left\{ \begin{aligned} S'(t) &= \Lambda - \delta S(t) - S(t) \sum_{l=1}^{4} \Gamma_l(t), \\ \frac{\partial i_m(t, \tau, v_0)}{\partial t} + \frac{\partial i_m(t, \tau, v_0)}{\partial \tau} &= -(\delta_m + \tilde{\mu}_m(\tau, v_0)) i_m(t, \tau, v_0), \\ \frac{\partial S_m(t, \tau, a_0)}{\partial t} + \frac{\partial S_m(t, \tau, a_0)}{\partial \tau} &= -\left(\delta + \sum_{s=1}^{4} \Gamma_{ms}(t, \tau, a_0)\right) S_m(t, \tau, a_0), \\ \frac{\partial i_{mn}(t, \tau, v_0, a_0)}{\partial t} + \frac{\partial i_{mn}(t, \tau, v_0, a_0)}{\partial \tau} &= -(\delta_{mn} + \tilde{\mu}_{mn}(\tau, v_0, a_0)) i_{mn}(t, \tau, v_0, a_0), \\ R'(t) &= -\delta R(t) + \sum_{l,s=1}^{4} \int_{0}^{+\infty} \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} \tilde{\mu}_{ls}(\tau, v_0, a_0) i_{ls}(t, \tau, v_0, a_0) da_0 dv_0 d\tau. \end{aligned} \right. \quad (13) $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Diagram of the epidemiological model with the four dengue serotypes and a maximum of two successive infections.</figcaption></figure>

The initial conditions are

$$ 
\begin{array}{l}
S(0) = S_0, \quad i_m(0, \tau, v_0) = i_{0m}(\tau, v_0), \quad S_m(0, \tau, a_0) = S_{0m}(\tau, a_0), \\
i_{mn}(0, \tau, v_0, a_0) = i_{0mn}(\tau, v_0, a_0), \quad R(0) = R_0.
\end{array} 
\quad (14) $$

The recovery rates depend on immunological variables, with

$$ 
\begin{array}{l}
\tilde{\mu}_m(\tau, v_0) := \mu_m(v(\tau, v_0), a(\tau, v_0)) \\
\text{and} \\
\tilde{\mu}_{mn}(\tau, v_0, a_0) := \mu_{mn}(v(\tau, v_0, a_0), a(\tau, v_0, a_0), b(\tau, v_0, a_0)),
\end{array}
$$

where $(v(\tau, v_0), a(\tau, v_0))$ and $(v(\tau, v_0, a_0), a(\tau, v_0, a_0), b(\tau, v_0, a_0))$ are the solutions of (7) and (10), respectively. We assume that all the functions $\tilde{\beta}_m, \tilde{\beta}_{mn}, \tilde{\mu}_m, \tilde{\mu}_{mn}$ and $\tilde{\alpha}_{mn}$ are essentially bounded on their respective domains, and that the initial functions $i_{0m}, S_{0m}$ and $i_{0mn}$ are integrable on their respective domains. As the solutions of (7), (8) and (10) are bounded, it suffices to assume that the functions $\beta_m, \beta_{mn}, \mu_m, \mu_{mn}$ and $\alpha_{mn}$ are continuous in order to ensure the essential boundedness of their corresponding functions. The probabilities of an infected individual having a viral load $v_0$ or antibodies $a_0$ from the previous infection $m$ at the beginning of the current infection $n$ are as follows

$$ 
\left\{ 
\begin{array}{l}
\tilde{p}_m(\tau, v'_0, v_0) = p_m(v(\tau, v'_0), v_0), \\
\tilde{p}_{mn}(\tau, v'_0, v_0) = p_{mn}(v(\tau, v'_0), v_0), \\
\tilde{q}_m(\tau, v_0, a_0) = q_m(a(\tau, v_0), a_0), \\
\tilde{r}_m(\tau, a'_0, a_0) = r_m(\tilde{a}(\tau, a'_0), a_0),
\end{array}
\right. 
$$

with

$$ \int_{0}^{+\infty} p_m(v, v_0)dv_0 = \int_{0}^{+\infty} p_{mn}(v, v_0)dv_0 = 1, $$
and

$$ \int_{\tilde{a}}^{+\infty} q_m(a, a_0) da_0 = \int_{\tilde{a}}^{+\infty} r_m(\tilde{a}, a_0) da_0 = 1. $$

In fact, $p_m(v, v_0)$ corresponds to the probability that an individual infected with serotype $m$, with a viral load $v$, infects susceptible individuals and transmits the viral load $v_0$. This indirectly takes into account the fact that the disease is transmitted by mosquitoes. Similarly, $q_m(a, a_0)$ is the probability that an infected individual (serotype $m$), with antibody level $a$, recovers with antibody level $a_0$. The susceptible individuals to serotype $n$, previously infected by serotype $m$, with a pre-existing antibody level $\tilde{a}$, are infected (by serotype $n$), with an initial pre-existing antibody level $a_0 > \tilde{a}$, with a probability $r_m(\tilde{a}, a_0)$. This also takes into account the decrease in antibodies during the recovery from serotype $m$ through the use of Expression (9). The probability that infected individuals (with serotype $n$ and previously with serotype $m$) with a viral load $v$, infect susceptible individuals and transmit the viral load $v_0$ is $p_{mn}(v, v_0)$. Newly infected individuals with serotype $m$ and an initial viral load $v_0$ are

$$ i_m(t, 0, v_0) = S(t) \left( \int_0^{+\infty} \int_0^{+\infty} \tilde{p}_m(\tau, v_0, v_0) \tilde{\beta}_m(\tau, v_0') i_m(t, \tau, v_0') dv_0' d\tau + \sum_{l=1}^4 \int_0^{+\infty} \int_0^{+\infty} \int_{\tilde{a}}^{+\infty} \tilde{p}_{lm}(\tau, v_0', v_0) \tilde{\beta}_{lm}(\tau, v_0', a_0) i_{lm}(t, \tau, v_0', a_0) da_0 dv_0' d\tau \right). \quad (15) $$

Newly recovered individuals from serotype $m$, with an initial antibody level $a_0$, are

$$ S_m(t, 0, a_0) = \int_0^{+\infty} \int_0^{+\infty} \tilde{q}_m(\tau, v_0, a_0) \tilde{\mu}_m(\tau, v_0) i_m(t, \tau, v_0) dv_0 d\tau. \quad (16) $$

Newly reinfected individuals with serotype $n$ (after a primary infection with $m$), with initial viral load $v_0$ and primary antibody level $a_0$, are

$$ i_{mn}(t, 0, v_0, a_0) = \left( \int_0^{+\infty} \int_{\tilde{a}}^{+\infty} \tilde{r}_m(\tau, a_0', a_0) \tilde{\alpha}_{mn}(\tau, a_0') S_m(t, \tau, a_0') d\tau da_0' \right) \times \left( \int_0^{+\infty} \int_0^{+\infty} \tilde{p}_n(\tau, v_0', v_0) \tilde{\beta}_n(\tau, v_0') i_n(t, \tau, v_0') dv_0' d\tau + \sum_{l=1}^4 \int_0^{+\infty} \int_0^{+\infty} \int_{\tilde{a}}^{+\infty} \tilde{p}_{ln}(\tau, v_0', v_0) \tilde{\beta}_{ln}(\tau, v_0', a_0') i_{ln}(t, \tau, v_0', a_0') da_0' dv_0' d\tau \right). \quad (17) $$

**Proposition 3.1.** The nested immuno-epidemiological system formed by the equations (7), (8), (10), (13), (14), (15), (16), and (17) is well-posed. Further-
more, the total number of individuals in the population at time t,

$$
\begin{aligned}
N(t) = S(t) &+ \sum_{m=1}^{4} \left[ \int_{0}^{+\infty} \int_{0}^{+\infty} i_m(t, \tau, v_0) dv_0 d\tau \right. \\
& \qquad \left. + \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} S_m(t, \tau, a_0) da_0 d\tau \right. \\
& \qquad \left. + \sum_{n=1}^{4} \int_{0}^{+\infty} \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} i_{mn}(t, \tau, v_0, a_0) da_0 dv_0 d\tau \right] + R(t),
\end{aligned}
$$

satisfies

$$
\limsup_{t \to +\infty} N(t) \le \frac{\Lambda}{\delta} \quad \text{with} \quad \delta = \min \{\delta, \delta_m, \delta_{mn} \mid m, n \in J\}.
$$

Proof. As established by the conditions above, the nested immuno-epidemiological system is well-posed and represents an extension of the age-structured Kermack-McKendrick model (see, for instance, [9, 21, 22, 23]). On the other hand, the derivative of N(t), the total population at time t, is given by

$$
\begin{aligned}
N'(t) = S'(t) &+ \sum_{m=1}^{4} \left[ \int_{0}^{+\infty} \int_{0}^{+\infty} \frac{\partial}{\partial t} i_{m}(t, \tau, v_{0}) dv_{0} d\tau \right. \\
& \qquad \left. + \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} \frac{\partial}{\partial t} S_{m}(t, \tau, a_{0}) da_{0} d\tau \right. \\
& \qquad \left. + \sum_{n=1}^{4} \int_{0}^{+\infty} \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} \frac{\partial}{\partial t} i_{mn}(t, \tau, v_{0}, a_{0}) da_{0} dv_{0} d\tau \right] + R'(t).
\end{aligned}
$$

Using System (13) and bringing together mortality terms, we get

$$
\begin{aligned}
N'(t) ={}& \Lambda - \left( \delta S(t) + \delta R(t) \right. \\
& \qquad \left. + \sum_{m=1}^{4} \left[ \delta_m \int_0^{+\infty} \int_0^{+\infty} i_m(t, \tau, v_0) + \delta \int_0^{+\infty} \int_{\bar{a}}^{+\infty} S_m(t, \tau, a_0) \right. \right. \\
& \qquad \left. \left. + \sum_{n=1}^{4} \delta_{mn} \int_0^{+\infty} \int_0^{+\infty} \int_{\bar{a}}^{+\infty} i_{mn}(t, \tau, v_0, a_0) \right] \right) \\
& - S(t) \sum_{m=1}^{4} \left( \lambda_m(t) + \sum_{n=1}^{4} \lambda_{mn}(t) \right) \\
& - \sum_{m=1}^{4} \left[ \int_0^{+\infty} \int_0^{+\infty} \tilde{\mu}_m(\tau, v_0) i_m(t, \tau, v_0) dv_0 d\tau \right. \\
& \qquad \left. + \int_0^{+\infty} \int_{\bar{a}}^{+\infty} S_m(t, \tau, a_0) \sum_{n=1}^{4} \tilde{\alpha}_{mn}(\tau, a_0) \left( \lambda_n(t) + \sum_{l=1}^{4} \lambda_{nl}(t) \right) da_0 d\tau \right. \\
& \qquad \left. + T_{im} + T_{S_m} + \sum_{n=1}^{4} T_{imn} \right],
\end{aligned}
$$
where

$$T_{im} := \int_{0}^{+\infty} \int_{0}^{+\infty} \frac{\partial i_m(t, \tau, v_0)}{\partial \tau} dv_0 d\tau, \quad T_{Sm} := \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} \frac{\partial S_m(t, \tau, a_0)}{\partial \tau} da_0 d\tau$$

and

$$T_{imn} := \int_{0}^{+\infty} \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} \frac{\partial i_{mn}(t, \tau, v_0, a_0)}{\partial \tau} da_0 dv_0 d\tau.$$

From the boundary condition (16), we obtain

$$T_{Sm} = - \int_{\bar{a}}^{+\infty} S_m(t, 0, a_0) da_0 = - \int_{0}^{+\infty} \int_{0}^{+\infty} \tilde{\mu}_m(\tau, v_0) i_m(t, \tau, v_0) dv_0 d\tau.$$

By switching the integrals we get

$$T_{im} = \int_{0}^{+\infty} \int_{0}^{+\infty} \frac{\partial i_m(t, \tau, v_0)}{\partial \tau} d\tau dv_0 = - \int_{0}^{+\infty} i_m(t, 0, v_0) dv_0.$$

Using the boundary condition (15), we obtain

$$T_{im} = - \int_{0}^{+\infty} S(t) \left( \int_{0}^{+\infty} \int_{0}^{+\infty} \tilde{p}_m(\tau, v'_0, v_0) \tilde{\beta}_m(\tau, v'_0) i_m(t, \tau, v'_0) dv'_0 d\tau + \sum_{n=1}^{4} \int_{0}^{+\infty} \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} \tilde{p}_{nm}(\tau, v'_0, v_0) \tilde{\beta}_{nm}(\tau, v'_0, a_0) i_{nm}(t, \tau, v'_0, a_0) da_0 dv'_0 d\tau \right) dv_0.$$

We switch again the integrals and use the condition on probabilities, to obtain

$$T_{im} = -S(t) \left( \int_{0}^{+\infty} \int_{0}^{+\infty} \tilde{\beta}_{m}(\tau, v'_{0}) i_{m}(t, \tau, v'_{0}) dv'_{0} d\tau + \sum_{n=1}^{4} \int_{0}^{+\infty} \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} \tilde{\beta}_{nm}(\tau, v'_{0}, a_{0}) i_{nm}(t, \tau, v'_{0}, a_{0}) da_{0} dv'_{0} d\tau \right).$$

We recognize $\lambda_m(t)$ and $\lambda_{nm}(t)$, so

$$T_{im} = -S(t) \left( \lambda_m(t) + \sum_{n=1}^{4} \lambda_{nm}(t) \right).$$

We also obtain that

$$T_{imn} = - \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} \tilde{\alpha}_{mn}(\tau, a'_{0}) S_{m}(t, \tau, a'_{0}) \left( \lambda_{n}(t) + \sum_{l=1}^{4} \lambda_{ln}(t) \right) da'_{0} d\tau.$$
Adding up all these terms, with $m, n \in J$, we obtain

$$ N'(t) = \Lambda - \left( \delta S(t) + \delta R(t) + \sum_{m=1}^{4} \left[ \delta_m \int_{0}^{+\infty} \int_{0}^{+\infty} i_m(t, \tau, v_0) + \delta \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} S_m(t, \tau, a_0) + \sum_{n=1}^{4} \delta_{mn} \int_{0}^{+\infty} \int_{0}^{+\infty} \int_{\bar{a}}^{+\infty} i_{mn}(t, \tau, v_0, a_0) \right] \right), \\ \leq \Lambda - \underline{\delta} N(t). $$

We conclude that

$$ \limsup_{t \to +\infty} N(t) \leq \frac{\Lambda}{\underline{\delta}}. $$

□

### 4. Steady states of nested model with two serotypes

We consider only two fixed serotypes of $J := \{1, \dots, 4\}$. Without loss of generality, we call these two serotypes 1 and 2. We do not allow reinfection by the same serotype ($\tilde{\alpha}_{mm}(\tau, a_0) = 0$, for $m \in \{1, 2\}$ and all $\tau, a_0$) and we set $\tilde{\alpha}_{12}(\tau, a_0) = \tilde{\alpha}_{21}(\tau, a_0) = 1$ (see Figure 5). Furthermore, as $R$ is not involved in the other equations, it can be derived from $i_{12}$ and $i_{21}$ and we can therefore omit it from the analysis. As $i_1$ and $i_2$ (resp. $i_{12}$ and $i_{21}$) have similar expressions, we will only work with $i_m$ (resp. $i_{mn}$) where $m, n \in \{1, 2\}$ and $m \neq n$, when it is relevant. A steady state $\mathcal{E}^*$ of the epidemiological model is given by

$$ \mathcal{E}^* := (S^*, i_1^*(\tau, v_0), i_2^*(\tau, v_0), S_1^*(\tau, a_0), S_2^*(\tau, a_0), i_{12}^*(\tau, v_0, a_0), i_{21}^*(\tau, v_0, a_0)). $$

It satisfies the system, for $m, n \in \{1, 2\}$ and $m \neq n$,

$$ \left\{ \begin{array}{l} 0 = \Lambda - \delta S^* - S^* (\lambda_1^* + \lambda_2^* + \lambda_{12}^* + \lambda_{21}^*), \\ \frac{\partial i_m^*(\tau, v_0)}{\partial \tau} = -(\delta_m + \tilde{\mu}_m(\tau, v_0)) i_m^*(\tau, v_0), \\ \frac{\partial S_m^*(\tau, a_0)}{\partial \tau} = -\delta S_m^*(\tau, a_0) - S_m^*(\tau, a_0) (\lambda_n^* + \lambda_{mn}^*), \\ \frac{\partial i_{mn}^*(\tau, v_0, a_0)}{\partial \tau} = -(\delta_{mn} + \tilde{\mu}_{mn}(\tau, v_0, a_0)) i_{mn}^*(\tau, v_0, a_0), \end{array} \right. \quad (18) $$

where

$$ \left\{ \begin{array}{l} \lambda_m^* = \int_0^\infty \int_0^\infty \tilde{\beta}_m(\tau, v_0) i_m^*(\tau, v_0) d\tau dv_0, \\ \lambda_{mn}^* = \int_0^\infty \int_0^\infty \int_{\bar{a}}^\infty \tilde{\beta}_{mn}(\tau, v_0, a_0) i_{mn}^*(\tau, v_0, a_0) da_0 dv_0 d\tau. \end{array} \right. \quad (19) $$

We integrate the last three equations of (18) to obtain

$$ \left\{ \begin{array}{l} i_m^*(\tau, v_0) = i_m^*(0, v_0)\pi_m(\tau, v_0), \\ S_m^*(\tau, a_0) = S_m^*(0, a_0)e^{-(\delta+\lambda_n^*+\lambda_{mn}^*)\tau}, \\ i_{mn}^*(\tau, v_0, a_0) = i_{mn}^*(0, v_0, a_0)\pi_{mn}(\tau, v_0, a_0), \end{array} \right. \quad (20) $$
<figure><img src="image_6.png" /><figcaption>Figure 5: Diagram of the model with only two serotypes and no reinfection with the same serotype.</figcaption></figure>

with

$$ \pi_m(\tau, v_0) := e^{-(\delta_m \tau + \int_0^\tau \tilde{\mu}_m(\theta, v_0) d\theta)}, \quad \pi_{mn}(\tau, v_0, a_0) := e^{-(\delta_{mn} \tau + \int_0^\tau \tilde{\mu}_{mn}(\theta, v_0, a_0) d\theta)}. $$

$\pi_m(\tau, v_0)$ (respectively, $\pi_{mn}(\tau, v_0, a_0)$) represents the probability of still being infected, after a time $\tau$, with serotype $m$ (respectively, $n$) after a first infection (respectively, second infection following a first infection with serotype $m$).

For the sake of simplicity, we may sometimes make the following assumptions:

$$ (H1) \quad \tilde{p}_m(\tau, v'_0, v_0) = p_m(v_0), \text{ with } \int_0^{+\infty} p_m(v_0) dv_0 = 1. $$

$$ (H2) \quad \tilde{p}_{mn}(\tau, v'_0, v_0) = p_{mn}(v_0), \text{ with } \int_0^{+\infty} p_{mn}(v_0) dv_0 = 1. $$

$$ (H3) \quad \tilde{r}_m(\tau, a'_0, a_0) = r_m(a_0), \text{ with } \int_{\bar{a}}^{+\infty} r_m(a_0) da_0 = 1. $$

$$ (H4) \quad \tilde{p}_{mn}(\tau, v'_0, v_0) = \tilde{p}_m(\tau, v'_0, v_0). $$

Assumptions (H1) and (H2) imply that the viral load $v_0$ of a newly infected individual is independent of the source of infection. This is biologically reasonable, especially in the case of dengue, which is a vector-borne disease, meaning transmission does not occur directly from one individual to another but rather through the bite of a mosquito. Hypothesis (H3) states that the initial pre-existing antibody level $a_0$ of an individual during secondary infection is independent of their antibody level in the previously recovered compartment. This
hypothesis may hold when the time between two infections is sufficiently long when compared to the duration of infection, allowing antibody levels to approach equilibrium. Under assumption (H4), we assume that the level of viral load transmission does not depend on the nature of the individual who transmits the disease (i.e., whether it is primary or secondary). This also corresponds to the fact that dengue is transmitted indirectly. We also introduce the following notations, for $m, n \in \{1, 2\}$ and $z \in \mathbb{C}$:

$$ (\text{N1}) \quad C_m^{p_m} := \int_0^\infty \left( \int_0^\infty \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau \right) p_m(v_0) dv_0 \quad \text{and} \\ C_m^{p_{nm}} = \int_0^\infty \left( \int_0^\infty \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau \right) p_{nm}(v_0) dv_0. $$

$$ (\text{N2}) \quad C_{mn}^{r_m, p_n} = \int_0^\infty \int_0^\infty \int_{\bar{a}}^\infty \tilde{\beta}_{mn}(\tau, v_0, a_0) \pi_{mn}(\tau, v_0, a_0) r_m(a_0) p_n(v_0) dv_0 d\tau da_0 \quad \text{and} \\ C_{mn}^{r_m, p_{mn}} = \int_0^\infty \int_0^\infty \int_{\bar{a}}^\infty \tilde{\beta}_{mn}(\tau, v_0, a_0) \pi_{mn}(\tau, v_0, a_0) r_m(a_0) p_{mn}(v_0) dv_0 d\tau da_0. $$

$$ (\text{N3}) \quad \tilde{C}_m^{p_m} = \int_0^\infty \left( \int_0^\infty \tilde{\mu}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau \right) p_m(v_0) dv_0 \quad \text{and} \\ \tilde{C}_{mn}^{p_{mn}} = \int_0^\infty \left( \int_0^\infty \tilde{\mu}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau \right) p_{nm}(v_0) dv_0. $$

$$ (\text{N4}) \quad \mathcal{R}_0^{(m)} := \frac{\Lambda}{\delta} C_m^{p_m} \quad \text{and} \quad \mathcal{R}_0^{(mn)} := \frac{\Lambda}{\delta} C_{mn}^{r_m, p_n}. $$

$$ (\text{N5}) \quad A_{mn} := \frac{C_{mn}^{r_m, p_n} \tilde{C}_m^{p_m}}{C_m^{p_m}}. $$

$$ (\text{N6}) \quad M_m := \frac{1}{C_m^{p_m} A_{mn}} = \frac{1}{C_{mn}^{r_m, p_n} \tilde{C}_m^{p_m}}. $$

$$ (\text{N7}) \quad \begin{aligned} \mathcal{G}_m^{p_m}(z) &= \int_0^{+\infty} \left( \int_0^{+\infty} \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) e^{-z\tau} d\tau \right) p_m(v_0) dv_0 \quad \text{and} \\ \mathcal{G}_{mn}^{r_m p_n}(z) &= \int_0^{+\infty} \int_{\bar{a}}^{+\infty} \left( \int_0^{+\infty} \tilde{\beta}_{mn}(\tau, v_0, a_0) \pi_{mn}(\tau, v_0, a_0) e^{-z\tau} d\tau \right) r_m(a_0) p_n(v_0) da_0 dv_0. \end{aligned} $$

In notation (N1), the constant $C_m^{p_m}$ represents the average of an infective power $\int_0^\infty \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau$ when $v_0$ is randomly drawn according to the probability distribution $p_m$, that is, $C_m^{p_m} = \mathbb{E}[\int_0^\infty \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau]$. The same interpretation holds for the other constants $C_m^{p_{nm}}$, $C_{mn}^{r_m, p_n}$, $C_{mn}^{r_m, p_{mn}}$, and the functions $\mathcal{G}_m^{p_m}$ and $\mathcal{G}_{mn}^{r_m p_n}$. Similarly, $\tilde{C}_m^{p_m}$ and $\tilde{C}_m^{p_{mn}}$, in notation (N3), may be considered as the expected values of the proportions of individuals infected with serotype $m$ that recover. The expression for $\mathcal{R}_0^{(m)}$, as given by (N4), corresponds to the basic reproduction number of the model, considering only infections with serotype $m$. The constant $A_{mn}$ may be interpreted as a ratio between the infective power of an individual with serotype $m$ and the infective power of one infected with serotype $n$, and who had serotype $m$ previously. The notation (N6) is introduced to simplify the writing. Under the hypothesis (H4), the constants in (N1) are equal and independent of $n$. We then denote

$$ C_m := C_m^{p_m} = C_m^{p_{nm}} = \int_0^\infty \left( \int_0^\infty \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau \right) p_m(v_0) dv_0. \quad (21) $$
Similarly, we denote the constants in (N2) by

$$ 
\begin{aligned}
C_{mn} &:= C_{mn}^{r_m, p_n} = C_{mn}^{r_m, p_{mn}}, \\
&= \int_0^\infty \int_{\bar{a}}^\infty \left( \int_0^\infty \tilde{\beta}_{mn}(\tau, v_0, a_0) \pi_{mn}(\tau, v_0, a_0) d\tau \right) r_m(a_0) p_n(v_0) da_0 dv_0,
\end{aligned} 
\quad (22) $$

and the ones of (N3) by

$$ \tilde{C}_m := \tilde{C}_m^{p_m} = \tilde{C}_m^{p_{mn}} = \int_0^\infty \left( \int_0^\infty \tilde{\mu}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau \right) p_m(v_0) dv_0. \quad (23) $$

With these expressions, we obtain

$$ \mathcal{R}_0^{(m)} = \frac{\Lambda}{\delta} \int_0^\infty \left( \int_0^\infty \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau \right) p_m(v_0) dv_0, \quad (24) $$

and

$$ \mathcal{R}_0^{(mn)} = \frac{\Lambda}{\delta} \int_0^\infty \int_{\bar{a}}^\infty \left( \int_0^\infty \tilde{\beta}_{mn}(\tau, v_0, a_0) \pi_{mn}(\tau, v_0, a_0) d\tau \right) r_m(a_0) p_n(v_0) da_0 dv_0. \quad (25) $$

$\mathcal{R}_0^{(mn)}$ represents the basic reproduction number for individuals who were initially infected with serotype $m$ and later infected with serotype $n$.

### 4.1. Disease-free equilibrium and its stability

A disease-free equilibrium, denoted $\mathcal{E}^0$, describes a state where the disease is absent in the population. We introduce the following notation

$$ \mathcal{R}_0 := \max_{m \in \{1, 2\}} \left( \mathcal{R}_0^{(m)} \right), \quad (26) $$

where $\mathcal{R}_0^{(m)}$ is given by (24).

**Theorem 4.1.** Assume that (H1) is satisfied. Then, a disease-free equilibrium always exists and is given by $\mathcal{E}^0 = (\Lambda/\delta, 0_{\mathbb{R}^6})$. Moreover, $\mathcal{E}^0$ is locally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

*Proof.* Let $\mathcal{E}^0 := (S^0, 0, 0, S_1^0(\tau, a_0), S_2^0(\tau, a_0), 0, 0)$ be a disease-free equilibrium. We have from (19) that $\lambda_1^0 = \lambda_2^0 = \lambda_{12}^0 = \lambda_{21}^0 = 0$. Then, from the third equation of System (20) and the boundary conditions, we obtain $S_1^0(\tau, a_0) = S_2^0(\tau, a_0) = 0$, while the first equation gives $S^0 = \Lambda/\delta$. Thus, there exists a unique disease-free equilibrium $\mathcal{E}^0 := (\Lambda/\delta, 0_{\mathbb{R}^6})$. To prove the local asymptotic stability of $\mathcal{E}^0$, we linearize the system by putting $\hat{S}(t) = S(t) - \Lambda/\delta$.
Then, for $m \neq n \in \{1, 2\}$,

$$
\left\{
\begin{aligned}
\hat{S}'(t) &= -\delta\hat{S}(t) - \frac{\Lambda}{\delta}(\lambda_1(t) + \lambda_2(t) + \lambda_{12}(t) + \lambda_{21}(t)), \\
\frac{\partial i_m}{\partial t}(t, \tau, v_0) + \frac{\partial i_m}{\partial \tau}(t, \tau, v_0) &= -(\delta_m + \tilde{\mu}_m(\tau, v_0)) i_m(t, \tau, v_0), \\
\frac{\partial S_m}{\partial t}(t, \tau, a_0) + \frac{\partial S_m}{\partial \tau}(t, \tau, a_0) &= -\delta S_m(t, \tau, a_0), \\
\frac{\partial i_{mn}}{\partial t}(t, \tau, v_0, a_0) + \frac{\partial i_{mn}}{\partial \tau}(t, \tau, v_0, a_0) &= -(\delta_{mn} + \tilde{\mu}_{mn}(\tau, v_0, a_0)) i_{mn}(t, \tau, v_0, a_0),
\end{aligned}
\right.
$$

where $\lambda_m(t)$ and $\lambda_{mn}(t)$ are given by (11) and (12). Using (H1), the boundary conditions become

$$
\begin{aligned}
i_m(t, 0, v_0) &= \frac{\Lambda}{\delta} \times \left( p_m(v_0)\lambda_m(t) \right. \\
&\qquad \left. + \int_0^{+\infty} \int_0^{+\infty} \int_{\bar{a}}^{+\infty} \bar{p}_{nm}(\tau, v'_0, v_0) \tilde{\beta}_{nm}(\tau, v'_0, a_0) i_{nm}(t, \tau, v'_0, a_0) da_0 \, dv'_0 \, d\tau \right), \\
S_m(t, 0, a_0) &= \int_0^{+\infty} \int_0^{+\infty} \tilde{q}_m(\tau, v_0, a_0) \tilde{\mu}_m(\tau, v_0) i_m(t, \tau, v_0) dv_0 \, d\tau,
\end{aligned}
$$

and

$i_{mn}(t, 0, v_0, a_0) = 0.$

We consider the exponential solutions, with $z \in \mathbb{C}$,

$$
\begin{aligned}
\hat{S}(t) &= e^{zt}\hat{S}, \quad i_m(t, \tau, v_0) = e^{zt}i_m(\tau, v_0), \quad S_m(t, \tau, a_0) = e^{zt}S_m(\tau, a_0), \\
i_{mn}(t, \tau, v_0, a_0) &= e^{zt}i_{mn}(\tau, v_0, a_0),
\end{aligned}
$$

of the linearized system:

$$
\left\{
\begin{array}{l}
(z + \delta) \dot{S} = -\frac{\Lambda}{\delta} (\lambda_1 + \lambda_2 + \lambda_{12} + \lambda_{21}), \\
\\
\frac{\partial i_m}{\partial \tau}(\tau, v_0) + z i_m(\tau, v_0) = -(\delta_m + \tilde{\mu}_m(\tau, v_0)) i_m(\tau, v_0), \\
\\
\frac{\partial S_m}{\partial \tau}(\tau, a_0) + z S_m(\tau, a_0) = -\delta S_m(\tau, a_0), \\
\\
\frac{\partial i_{mn}}{\partial \tau}(\tau, v_0, a_0) + z i_{mn}(\tau, v_0, a_0) = -(\delta_{mn} + \tilde{\mu}_{mn}(\tau, v_0, a_0)) i_{mn}(\tau, v_0, a_0).
\end{array}
\right.
\qquad (27)
$$

We denote without ambiguity $\lambda_m(t) = e^{zt}\lambda_m$ and $\lambda_{mn}(t) = e^{zt}\lambda_{mn}$, with

$$
\begin{align}
\lambda_m &:= \int_0^{+\infty} \int_0^{+\infty} \tilde{\beta}_m(\tau, v_0) i_m(\tau, v_0) \, dv_0 \, d\tau, \nonumber \\
\lambda_{mn} &:= \int_0^{+\infty} \int_0^{+\infty} \int_{\bar{a}}^{+\infty} \tilde{\beta}_{mn}(\tau, v_0, a_0) i_{mn}(\tau, v_0, a_0) \, da_0 \, dv_0 \, d\tau. \tag{28}
\end{align}
$$
Noting that for the linearized system, we have the following null boundary conditions $i_{mn}(0, v_0, a_0) = 0$, we may integrate the last equations to obtain

$$
\begin{align*}
i_m(\tau, v_0) &= \pi_m(\tau, v_0)e^{-z\tau}i_m(0, v_0), \\
S_m(\tau, a_0) &= e^{-(z+\delta)\tau}S_m(0, a_0), \\
i_{mn}(\tau, v_0, a_0) &= \pi_{mn}(\tau, v_0, a_0)e^{-z\tau}i_{mn}(0, v_0, a_0) = 0.
\end{align*}
$$

Moreover, under the assumption $\tilde{p}_m(\tau, v'_0, v_0) = p_m(v_0)$, we can show that

$$
i_m(0, v_0) = \frac{\Lambda}{\delta} p_m(v_0) \lambda_m.
$$

We also have

$$
S_m(0, a_0) = \int_0^{+\infty} \int_0^{+\infty} \tilde{q}_m(\tau, v_0, a_0) \tilde{\mu}_m(\tau, v_0) i_m(\tau, v_0) dv_0 d\tau.
$$

Then, the first equation of (28) becomes

$$
\begin{align}
\lambda_m &= \int_0^{+\infty} \int_0^{+\infty} \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) e^{-z\tau} i_m(0, v_0) dv_0 d\tau, \nonumber \\
&= \frac{\Lambda}{\delta} \lambda_m \int_0^{+\infty} \left( \int_0^{+\infty} \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) e^{-z\tau} d\tau \right) p_m(v_0) dv_0, \tag{29} \\
&= \frac{\Lambda}{\delta} \mathcal{G}_m^{p_m}(z) \lambda_m, \nonumber
\end{align}
$$

where

$$
\mathcal{G}_m^{p_m}(z) = \int_0^{+\infty} \left( \int_0^{+\infty} \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) e^{-z\tau} d\tau \right) p_m(v_0) dv_0.
$$

The function $\mathcal{G}_m^{p_m}$ is strictly decreasing on $[0, +\infty)$, with $\mathcal{G}_m^{p_m}(0) = C_m$ (see (N1)) and $\lim_{z \to +\infty} \mathcal{G}_m^{p_m}(z) = 0$. If $\lambda_m \neq 0$ we have from Equation (29),

$$
\frac{\Lambda}{\delta} \mathcal{G}_m^{p_m}(z) = 1.
$$

Thus, a positive real solution exists if and only if $\Lambda/\delta \times \mathcal{G}_m^{p_m}(0) > 1$. This is equivalent to

$$
\mathcal{R}_0^{(m)} > 1.
$$

Hence, the disease-free steady state is unstable if $\max_{m \in \{1,2\}}(\mathcal{R}_0^{(m)}) > 1$, because there is a real positive eigenvalue. Now, suppose that $\mathcal{R}_0^{(m)} < 1$ and let $z = \mathcal{R}e(z) + \mathcal{I}m(z)i$ with $\mathcal{R}e(z) \ge 0$. Then,

$$
\begin{align*}
\left| \frac{\Lambda}{\delta} \mathcal{G}_m^{p_m}(z) \right| &\le \frac{\Lambda}{\delta} \int_0^{+\infty} \left( \int_0^{+\infty} \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) e^{-\mathcal{R}e(z)\tau} d\tau \right) p_m(v_0) dv_0, \\
&\le \frac{\Lambda}{\delta} \mathcal{G}_m^{p_m}(0) = \mathcal{R}_0^{(m)} < 1.
\end{align*}
$$
Hence $z$ with $\mathcal{R}e(z) \ge 0$ cannot be a solution of $\Lambda/\delta \times G_m^{p_m}(z) = 1$. Therefore, to satisfy Equation (29), we must have $\lambda_m = 0$. Then, with the first equation of System (27) we find $z = -\delta < 0$. This is in contradiction with $\mathcal{R}e(z) \ge 0$. We conclude that if $\max_{m \in \{1,2\}} (\mathcal{R}_0^{(m)}) < 1$ all eigenvalues have a negative real part and the disease-free steady state is locally asymptotically stable. $\square$

**Remark 4.2.** Thanks to Formula (26), we have $\mathcal{R}_0 < 1$ when both $\mathcal{R}_0^{(1)} < 1$ and $\mathcal{R}_0^{(2)} < 1$ are satisfied. In this case, no epidemic persists, and all infections disappear. We have $\mathcal{R}_0 > 1$ when at least one of $\mathcal{R}_0^{(1)}$ or $\mathcal{R}_0^{(2)}$ is greater than 1.

## 5. Endemic steady states

We now turn to the steady states with a non-null number of individuals infected with at least one of the two serotypes:

$$ \mathcal{E}^* := (S^*, i_1^*(\tau, v_0), i_2^*(\tau, v_0), S_1^*(\tau, a_0), S_2^*(\tau, a_0), i_{12}^*(\tau, v_0, a_0), i_{21}^*(\tau, v_0, a_0)). $$

In this section, we assume (H1), (H2), (H3) and (H4). We obtain the following expressions for $m, n \in \{1, 2\}, m \ne n$, at the endemic steady states:

$$ \left\{ \begin{aligned} \Lambda &= (\delta + \lambda_1^* + \lambda_2^* + \lambda_{12}^* + \lambda_{21}^*) S^*, \\ \lambda_m^* &= \int_0^\infty \int_0^\infty \tilde{\beta}_m(\tau, v_0) i_m^*(\tau, v_0) d\tau dv_0, \\ \lambda_{mn}^* &= \int_0^\infty \int_0^\infty \int_{\bar{a}}^\infty \tilde{\beta}_{mn}(\tau, v_0, a_0) i_{mn}^*(\tau, v_0, a_0) da_0 dv_0 d\tau, \end{aligned} \right. \quad (30) $$

and

$$ \left\{ \begin{aligned} i_m^*(\tau, v_0) &= i_m^*(0, v_0)\pi_m(\tau, v_0), \\ S_m^*(\tau, a_0) &= S_m^*(0, a_0)e^{-(\delta+\lambda_n^*+\lambda_{mn}^*)\tau}, \\ i_{mn}^*(\tau, v_0, a_0) &= i_{mn}^*(0, v_0, a_0)\pi_{mn}(\tau, v_0, a_0), \end{aligned} \right. $$

with the boundary conditions

$$ \left\{ \begin{aligned} i_m^*(0, v_0) &= S^* p_m(v_0) (\lambda_m^* + \lambda_{nm}^*), \\ i_{mn}^*(0, v_0, a_0) &= r_m(a_0) p_n(v_0) \bar{S}_m^*(\lambda_n^* + \lambda_{mn}^*), \\ S_m^*(0, a_0) &= \int_0^{+\infty} \int_0^{+\infty} \tilde{q}_m(\tau, v_0, a_0) \tilde{\mu}_m(\tau, v_0) i_m^*(\tau, v_0) dv_0 d\tau, \end{aligned} \right. $$

where $\bar{S}_m^*$ is the total population of recovered individuals from serotype $m$,

$$ \bar{S}_m^* = \int_0^\infty \int_{\bar{a}}^\infty S_m^*(\tau, a_0) d\tau da_0. $$

Using the boundary condition of $i_m^*$, and applying Fubini-Tonelli's theorem (which applies to all the integrals we consider, since they are positive), we obtain

$$ \lambda_m^* = \int_0^\infty \int_0^\infty \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) i_m^*(0, v_0) d\tau dv_0, $$
$$ 
\begin{aligned}
&= S^* \left[ \int_0^\infty \left( \int_0^\infty \tilde{\beta}_m(\tau, v_0) \pi_m(\tau, v_0) d\tau \right) p_m(v_0) dv_0 \right] (\lambda_m^* + \lambda_{nm}^*), \\
&= S^* C_m (\lambda_m^* + \lambda_{nm}^*),
\end{aligned}
$$

where $C_m$ is defined by (21). Doing the same for $\lambda_{mn}^*$, we obtain

$$ 
\begin{aligned}
\lambda_{mn}^* &= \int_0^\infty \int_0^\infty \int_{\bar{a}}^\infty \tilde{\beta}_{mn}(\tau, v_0, a_0) \pi_{mn}(\tau, v_0, a_0) i_{mn}^*(0, v_0, a_0) da_0 \, dv_0 \, d\tau, \\
&= C_{mn} (\lambda_n^* + \lambda_{mn}^*) \bar{S}_m^*.
\end{aligned}
$$

where $C_{mn}$ is given by (22). Now, using the second equation of (20), we can write

$$ 
\begin{aligned}
\bar{S}_m^* &= \left( \int_0^\infty e^{-(\delta+\lambda_n^*+\lambda_{nm}^*)\tau} d\tau \right) \left( \int_{\bar{a}}^\infty S_m^*(0, a_0) da_0 \right), \\
&= \frac{1}{\delta + \lambda_n^* + \lambda_{mn}^*} \int_{\bar{a}}^\infty S_m^*(0, a_0) da_0.
\end{aligned}
\quad (31) $$

The integral $\int_{\bar{a}}^\infty S_m^*(0, a_0) da_0$, corresponds to all the newly recovered individuals from serotype m. Using the boundary condition and the notations in (N3),

$$ 
\begin{aligned}
\int_{\bar{a}}^{\infty} S_{m}^{*}(0, a_{0}) da_{0} &= \int_{\bar{a}}^{\infty} \int_{0}^{\infty} \int_{0}^{\infty} q_{m}(\tau, v_{0}, a_{0}) \tilde{\mu}(\tau, v_{0}) i_{m}^{*}(\tau, v_{0}) d\tau \, dv_{0} \, da_{0}, \\
&= \int_{0}^{\infty} \int_{0}^{\infty} \tilde{\mu}(\tau, v_{0}) i_{m}^{*}(\tau, v_{0}) d\tau \, dv_{0}, \\
&= \tilde{C}_{m} S^{*}(\lambda_{m}^{*} + \lambda_{nm}^{*}).
\end{aligned}
$$

Combining all the previously obtained results yields the following system

$$ 
\left\{
\begin{array}{l}
[\delta + \lambda_m^* + \lambda_{nm}^* + \lambda_n^* + \lambda_{mn}^*] S^* = \Lambda, \\
\lambda_m^* = S^* C_m (\lambda_m^* + \lambda_{nm}^*), \\
\lambda_{mn}^* = S^* C_{mn}^n \frac{\lambda_n^* + \lambda_{mn}^*}{\delta + \lambda_n^* + \lambda_{mn}^*} \tilde{C}_m (\lambda_m^* + \lambda_{nm}^*), \\
\int_{\bar{a}}^{\infty} S_m^*(0, a_0) da_0 = S^* \tilde{C}_m (\lambda_m^* + \lambda_{nm}^*).
\end{array}
\right.
\quad (32) $$

Solving System (32) in $\mathbb{R}_+^6$ will yield a polynomial in $S^*$, whose roots will correspond one-to-one to the endemic steady states of the model, as we prove in lemma 5.5. There will be two types of endemic steady states: mono-endemic and bi-endemic. These two types are described and analyzed in the subsections 5.1 and 5.2. Since a complete study of bi-endemic steady states would be too laborious, we focus only on the number of bi-endemic equilibrium points, and list the possible bifurcation diagrams.

### 5.1. Mono-endemic steady state and its stability

**Definition 5.1.** A steady state is said to be mono-endemic for serotype m if it satisfies $\lambda_m^* \neq 0$ and $\lambda_n^* = 0$, for $n \neq m$, where $\lambda_m^*$ and $\lambda_n^*$ are given by (30), i.e., when only the serotype m is endemic in the population. Such an equilibrium is denoted $\mathcal{E}^m$.
In the following, and without loss of generality, we will consider a mono-endemic steady state for serotype $m = 1$.

**Theorem 5.2.** Assume (H1), (H2), (H3) and (H4). Then, the unique mono-endemic steady state for serotype 1 is

$$ \mathcal{E}^1 = (S^*, i_1^*(\tau, v_0), 0, S_1^*(\tau, a_0), 0, 0, 0), $$

with

$$ S^* = \frac{1}{C_1}, \quad i_1^*(\tau, v_0) = \frac{\delta}{C_1} (\mathcal{R}_0^{(1)} - 1) p_1(v_0) \pi_1(\tau, v_0) $$

and

$$ S_1^*(\tau, a_0) = \frac{\delta e^{-\delta \tau}}{C_1} (\mathcal{R}_0^{(1)} - 1) \int_0^{+\infty} \int_0^{+\infty} \tilde{q}_1(\tau, v_0, a_0) \tilde{\mu}_1(\tau, v_0) p_1(v_0) \pi_1(\tau, v_0) dv_0 d\tau, $$

where $C_1$ is given by (21). This mono-endemic steady state exists if and only if

$$ \mathcal{R}_0^{(1)} > 1. $$

Furthermore, suppose that

$$ \mathcal{R}_0^{(12)} < \frac{1}{\tilde{C}_1} \frac{\mathcal{R}_0^{(1)} - \mathcal{R}_0^{(2)}}{\mathcal{R}_0^{(1)} - 1}, \qquad (33) $$

where $\tilde{C}_1$, $\mathcal{R}_0^{(1)}$, $\mathcal{R}_0^{(2)}$ and $\mathcal{R}_0^{(12)}$ are given by (23), (24) and (25). Then, $\mathcal{E}^1$ is locally asymptotically stable. If the inequality in (33) is in the opposite direction (and strict), then the equilibrium $\mathcal{E}^1$ is unstable.

**Remark 5.3.** The two mono-endemic steady states

$$ \mathcal{E}^1 = (1/C_1, i_1^*(\tau, v_0), 0, S_1^*(\tau, a_0), 0, 0, 0), \quad \mathcal{E}^2 = (1/C_2, 0, i_2^*(\tau, v_0), 0, S_2^*(\tau, a_0), 0, 0), $$

exist simultaneously if and only if $\mathcal{R}_0^{(1)} > 1$ and $\mathcal{R}_0^{(2)} > 1$. Both steady states are unstable if

$$ \mathcal{R}_0^{(12)} > \frac{1}{\tilde{C}_1} \frac{\mathcal{R}_0^{(1)} - \mathcal{R}_0^{(2)}}{\mathcal{R}_0^{(1)} - 1} \quad \text{and} \quad \mathcal{R}_0^{(21)} > \frac{1}{\tilde{C}_2} \frac{\mathcal{R}_0^{(2)} - \mathcal{R}_0^{(1)}}{\mathcal{R}_0^{(2)} - 1}. $$

However, only the one with the largest basic reproduction number, say $\mathcal{R}_0^{(1)} > \mathcal{R}_0^{(2)}$, can be locally asymptotically stable (under the condition (33)). In this case, the other is necessarily unstable. In the symmetric case, where both serotypes have the same basic reproduction number, $\mathcal{R}_0^{(1)} = \mathcal{R}_0^{(2)} > 1$, we have $S^* = 1/C_1 = 1/C_2$. In this case, the condition (33) is never satisfied, and the two mono-endemic steady states $\mathcal{E}^1$ and $\mathcal{E}^2$ are both always unstable. If the two serotypes have the same basic reproduction number, either they both disappear or they persist together in the population in the form of a bi-endemic equilibrium (see Subsection 5.2). This situation is fundamentally different from
the classical framework of competitive exclusion. Furthermore, in this model, two different stable mono-endemic equilibrium points cannot coexist; instead, the one with the higher basic reproduction number will supplant the other or both will persist, forming a bi-endemic equilibrium (see Subsection 5.2).

*Proof.* Suppose that $\lambda_1^* \neq 0$ and $\lambda_2^* = 0$. This means that we are considering a mono-endemic steady state for serotype 1. In this case, from the second equation of System (32), we have $\lambda_{12}^* = 0$. Consequently, the last two equations lead to $S_2^* = \lambda_{21}^* = 0$. Equation (31) becomes

$$ \bar{S}_1^* = \frac{1}{\delta} \int_{\bar{a}}^{\infty} S_1^*(0, a_0) da_0. $$

The first two equations of System (32) give $(\delta + \lambda_1^*)S^* = \Lambda$ and $\lambda_1^* = S^*C_1\lambda_1^*$. Since $\lambda_1^* \neq 0$, we obtain

$$ S^* = \frac{1}{C_1}. $$

This gives

$$ \lambda_1^* = \delta (\mathcal{R}_0^{(1)} - 1) \quad \text{and} \quad \int_{\bar{a}}^{\infty} S_1^*(0, a_0) da_0 = \frac{\tilde{C}_1 \delta}{C_1} (\mathcal{R}_0^{(1)} - 1). $$

This mono-endemic steady state exists if and only if $\mathcal{R}_0^{(1)} > 1$. We linearize the system around this steady state by defining $\hat{S}(t) = S(t) - S^*$, $\hat{i}_1(t, \tau, v_0) = i_1(t, \tau, v_0) - i_1^*(\tau, v_0)$ and $\hat{S}_1(t, \tau, a_0) = S_1(t, \tau, a_0) - S_1^*(\tau, a_0)$. To simplify the notation, we retain $i_1$, $S_1$ and $\lambda_1$ instead of $\hat{i}_1$, $\hat{S}_1$ and $\hat{\lambda}_1$. Let $m, n \in \{1, 2\}$, with $m \neq n$. We have

$$ \left\{ \begin{aligned} \hat{S}'(t) &= -\delta \hat{S}(t) - S^*(\lambda_1(t) + \lambda_2(t) + \lambda_{12}(t) + \lambda_{12}(t)) - \hat{S}(t)\lambda_1^*, \\ \frac{\partial i_m}{\partial t}(t, \tau, v_0) + \frac{\partial i_m}{\partial \tau}(t, \tau, v_0) &= -(\delta_m + \tilde{\mu}_m(\tau, v_0)) i_m(t, \tau, v_0), \\ \frac{\partial S_1}{\partial t}(t, \tau, a_0) + \frac{\partial S_1}{\partial \tau}(t, \tau, a_0) &= -\delta S_1(t, \tau, a_0), \\ \frac{\partial S_2}{\partial t}(t, \tau, a_0) + \frac{\partial S_2}{\partial \tau}(t, \tau, a_0) &= -\delta S_2(t, \tau, a_0) - \lambda_1^* S_2(t, \tau, a_0), \\ \frac{\partial i_{mn}}{\partial t}(t, \tau, v_0, a_0) + \frac{\partial i_{mn}}{\partial \tau}(t, \tau, v_0, a_0) &= -(\mu + \tilde{\mu}_{mn}(\tau, v_0, a_0)) i_{mn}(t, \tau, v_0, a_0). \end{aligned} \right. $$

The linearized boundary conditions, under (H1), (H2), (H3) and (H4), for $m, n \in \{1, 2\}$ with $m \neq n$, are

$$ i_m(t, 0, v_0) = S(t)p_m(v_0)(\lambda_m^* + \lambda_{nm}^*) + S^*p_m(v_0)(\lambda_m(t) + \lambda_{nm}(t)), $$

$$ S_m(t, 0, a_0) = \int_0^{+\infty} \int_0^{+\infty} \tilde{q}_m(\tau, v_0, a_0) \tilde{\mu}_m(\tau, v_0) i_m(t, \tau, v_0) dv_0 d\tau, $$
and

$$i_{mn}(t, 0, v_0, a_0) = r_m(a_0)p_n(v_0)\bar{S}_m(t)[\lambda_n^* + \lambda_{mn}^*] \\ + r_m(a_0)p_n(v_0)\bar{S}_m^*[\lambda_n(t) + \lambda_{mn}(t)],$$

with

$$\bar{S}_m(t) = \int_0^\infty \int_{\bar{a}}^\infty S_m(t, \tau, a_0) da_0 d\tau \quad \text{and} \quad \bar{S}_m^* = \int_0^\infty \int_{\bar{a}}^\infty S_m^*(\tau, a_0) da_0 d\tau.$$

As in the study of the disease-free steady state, we plug eigenfunctions as an exponential solution of the form

$$\hat{S}(t) = e^{zt}\hat{S}, \quad i_m(t, \tau, v_0) = e^{zt}i_m(\tau, v_0), \quad S_m(t, \tau, a_0) = e^{zt}S_m(\tau, a_0), \\ i_{mn}(t, \tau, v_0, a_0) = e^{zt}i_{mn}(\tau, v_0, a_0), \quad \lambda_m(t) = e^{zt}\lambda_m, \quad \lambda_{mn}(t) = e^{zt}\lambda_{mn},$$

with $z \in \mathbb{C}$ an eigenvalue. Then, for $m, n \in \{1, 2\}$ with $m \neq n$, we have

$$\left\{
\begin{aligned}
z\hat{S} &= -\delta\hat{S} - S^*(\lambda_1 + \lambda_2 + \lambda_{21} + \lambda_{12}) - \hat{S}\lambda_1^*, \\
\frac{\partial i_m}{\partial \tau}(\tau, v_0) + z i_m(\tau, v_0) &= -(\delta_m + \tilde{\mu}_m(\tau, v_0)) i_m(\tau, v_0), \\
\frac{\partial S_1}{\partial \tau}(\tau, a_0) + z S_1(\tau, a_0) &= -\delta S_1(\tau, a_0), \\
\frac{\partial S_2}{\partial \tau}(\tau, a_0) + z S_2(\tau, a_0) &= -\delta S_2(\tau, a_0) - \lambda_1^* S_2(\tau, a_0), \\
\frac{\partial i_{mn}}{\partial \tau}(\tau, v_0, a_0) + z i_{mn}(\tau, v_0, a_0) &= -(\delta_{mn} + \tilde{\mu}_{mn}(\tau, v_0, a_0)) i_{mn}(\tau, v_0, a_0).
\end{aligned}
\right.
\quad (34)$$

Integrating the last equations, we obtain

$$\left\{
\begin{aligned}
i_m(\tau, v_0) &= \pi_m(\tau, v_0)e^{-z\tau}i_m(0, v_0), \\
S_1(\tau, a_0) &= e^{-(z+\delta)\tau}S_1(0, a_0), \\
S_2(\tau, a_0) &= e^{-(z+\delta+\lambda_1^*)\tau}S_2(0, a_0), \\
i_{mn}(\tau, v_0, a_0) &= \pi_{mn}(\tau, v_0, a_0)e^{-z\tau}i_{mn}(0, v_0, a_0).
\end{aligned}
\right.
$$

In particular, we have

$$i_m(0, v_0) = \hat{S}p_m(v_0)(\lambda_m^* + \lambda_{nm}^*) + S^*p_m(v_0)(\lambda_m + \lambda_{nm}), \\ i_{mn}(0, v_0, a_0) = r_m(a_0)p_n(v_0)[\bar{S}_m(\lambda_n^* + \lambda_{mn}^*) + \bar{S}_m^*(\lambda_n + \lambda_{mn})],$$

with $\lambda_m, \lambda_{mn}$ given by Equation (28). Then, using the same calculation methods as those applied to the disease-free steady state, we find, with notation (N7),

$$\begin{aligned}
\lambda_2 &= S^* \mathcal{G}_2^{p_2}(z) (\lambda_2 + \lambda_{12}), \\
\lambda_{12} &= \bar{S}_1^* \mathcal{G}_{12}^{r_1 p_2}(z) (\lambda_2 + \lambda_{12}), \\
\lambda_1 &= S^* \mathcal{G}_1^{p_1}(z) (\lambda_1 + \lambda_{21}) + \hat{S} \mathcal{G}_1^{p_1}(z) \lambda_1^*, \\
\lambda_{21} &= \bar{S}_2 \mathcal{G}_{21}^{r_2 p_1}(z) \lambda_1^*,
\end{aligned}
\quad (35)$$
where $G_m^{p_m}, G_{mn}^{r_m p_n}$ are decreasing functions of $z \in [0, +\infty)$ with $G_m^{p_m}(0) = C_m$ and $G_{mn}^{r_m p_n}(0) = C_{mn}$. Summing the two first equations of (35) yields the following relation

$$ \lambda_2 + \lambda_{12} = (S^* G_2^{p_2}(z) + \bar{S}_1^* G_{12}^{r_1 p_2}(z))(\lambda_{12} + \lambda_2) = (\lambda_2 + \lambda_{12})G(z). $$

The function

$$ G : z \mapsto S^* G_2^{p_2}(z) + \bar{S}_1^* G_{12}^{r_1 p_2}(z) $$

is a strictly decreasing function on $[0, +\infty)$, with $G(0) = S^*C_2 + \bar{S}_1^*C_{12}$ and $\lim_{z \to +\infty} G(z) = 0$. If $\lambda_2 + \lambda_{12} \neq 0$, we can simplify the previous equation, which gives the characteristic equation of the system:

$$ G(z) = 1, \quad z \in \mathbb{C}. \qquad (36) $$

Hence, a positive real solution of Equation (36) exists if and only if $G(0) > 1$. This condition may be rewritten, using $S^* = 1/C_1$, $\bar{S}_1^* = (\int_{\bar{a}}^{\infty} S_1^*(0, a_0) da_0) / \delta = \tilde{C}_1(R_0^{(1)} - 1) / C_1$ and notation (N4), as

$$ R_0^{(12)} > \frac{1}{\tilde{C}_1} \frac{R_0^{(1)} - R_0^{(2)}}{R_0^{(1)} - 1}. $$

Therefore, the mono-endemic equilibrium $\mathcal{E}^1$ has a real, positive eigenvalue if the inequality above is satisfied, making it unstable.

We now assume that the inequality (33) holds, i.e., $G(0) < 1$. Suppose that there exists an eigenvalue $z = \mathcal{R}(z) + i\mathcal{I}(z) \in \mathbb{C}$ with $\mathcal{R}(z) \ge 0$. We have proved in the study of the disease-free steady state that $|G_m^{p_m}(z)| \le G_m^{p_m}(\mathcal{R}(z))$. We have similar inequalities for $G_2^{p_2}$ and $G_{12}^{r_1 p_2}$. Thus,

$$ |G(z)| \le S^* G_2^{p_2}(\mathcal{R}(z)) + \bar{S}_1^* G_{12}^{r_1 p_2}(\mathcal{R}(z)) = G(\mathcal{R}(z)) \le G(0) < 1. $$

Then $G(z) \neq 1$. This leads to $\lambda_2 + \lambda_{12} = 0$. Then, using (35), we obtain that $\lambda_2 = \lambda_{12} = 0$, which implies that $i_2 = i_{12} = 0$. Moreover, we know that $S_2(\tau, a_0) = e^{-(z+\delta+\lambda_1^*)\tau} S_2(0, a_0)$, with

$$ S_2(0, a_0) = \int_0^\infty \int_0^\infty q_m(\tau, v_0, a_0) \mu_m(\tau, v_0) i_2(\tau, v_0) dv_0 d\tau. $$

Thus, we obtain $S_2 = 0$. Therefore $\bar{S}_2 = 0$ and from Equation (35) we obtain $\lambda_{21} = 0$. Then, the equation on $\lambda_1$ becomes

$$ \lambda_1 = S^* G_1^{p_1}(z) \lambda_1 + \hat{S} G_1^{p_1}(z) \lambda_1^*. $$

Using this equation and the first one of System (34), we find

$$ 1 = S^* G_1^{p_1}(z) - \frac{S^*}{z + \delta + \lambda_1^*} G_1^{p_1}(z) \lambda_1^* = S^* G_1^{p_1}(z) \frac{z + \delta}{z + \delta + \lambda_1^*}. $$
Using again $S^* = 1/C_1 = 1/G_1^{p_1}(0)$ and looking at the squared modulus, we find

$$1 = \left| \frac{G_1^{p_1}(z)}{G_1^{p_1}(0)} \times \frac{z + \delta}{z + \delta + \lambda_1^*} \right|^2 \le \left( \frac{G_1^{p_1}(\mathcal{R}e(z))}{G_1^{p_1}(0)} \right)^2 \frac{(\mathcal{R}e(z) + \delta)^2 + \mathcal{I}m(z)^2}{(\mathcal{R}e(z) + \delta + \lambda_1^*)^2 + \mathcal{I}m(z)^2} < 1.$$

This is because $\lambda_1^* > 0$ and $G_1^{p_1}(\mathcal{R}e(z))/G_1^{p_1}(0) \le 1$ since $\mathcal{R}e(z) \ge 0$ and $G_1^{p_1}$ is decreasing on $[0, +\infty)$. We can see that this is impossible. We conclude that the system has no eigenvalues in the right complex half-plane. $\square$

## 5.2. Bi-endemic steady states, backward bifurcation and invasion reproduction number

We are now searching for endemic steady states with coexisting serotypes, meaning steady states where both serotypes 1 and 2 are endemic in the population.

**Definition 5.4.** A steady state $\mathcal{E}^{mn}$ is said to be bi-endemic for serotypes $m, n \in \{1, 2\}$ with $m \neq n$, if it satisfies $\lambda_m^* \neq 0$ and $\lambda_n^* \neq 0$, i.e., when both serotypes are endemic in the population.

Note that, in the case of dengue, reaching a bi-endemic equilibrium is particularly dangerous, as it can lead to a significant increase in DHF cases due to the antibody-dependent enhancement phenomenon. Since the bi-endemic steady states cannot be expressed in a simple form for our model, we focus on determining the number of valid bi-endemic equilibrium points.

**Lemma 5.5.** The number of susceptible individuals in the bi-endemic steady states, denoted by $S^*$, corresponds to the roots within the interval

$$(r_+, \min\{1/C_1, 1/C_2\}),$$

of the cubic polynomial

$$P(x) = (x^2 - M_1 M_2 (1 - C_1 x)(1 - C_2 x)) \left(x - \frac{\Lambda}{\delta}\right) \\ + x \frac{1 - C_2 x}{C_1 A_{12}} (x + M_2(1 - C_1 x)) \\ + x \frac{1 - C_1 x}{C_2 A_{21}} (x + M_1(1 - C_2 x)),$$

where $A_{mn} = C_{mn} \tilde{C}_m / C_m$ is given by (N5) and $r_+$ is defined as the root, located in $(0, \min\{1/C_1, 1/C_2\})$, of the polynomial

$$Q(x) := x^2 - M_1 M_2 (1 - C_1 x)(1 - C_2 x).$$

Furthermore, a bi-endemic equilibrium is uniquely characterized by its susceptible component $S^*$, and all its other components are non-zero. In particular, $\lambda_{nm}^* \neq 0$, for $m, n \in \{1, 2\}$ with $m \neq n$.
Proof. We have from System (32), for $m, n \in \{1, 2\}$ with $m \neq n$,

$$ 
\left\{ 
\begin{aligned}
\lambda_m^* &= S^* C_m (\lambda_m^* + \lambda_{nm}^*), \\
\Lambda &= \delta S^* + \frac{\lambda_1^*}{C_1} + \frac{\lambda_2^*}{C_2}, \\
\lambda_{mn}^* &= A_{mn} \frac{\lambda_n^* \lambda_m^*}{C_n \delta S^* + \lambda_n^*}.
\end{aligned}
\right. 
$$

The first equation may be rewritten as

$$ \lambda_{nm}^* = \frac{1 - S^* C_m}{S^* C_m} \lambda_m^*. \qquad (37) $$

To obtain a non-negative steady state, we must ensure that

$$ S^* \leq \min\{1/C_1, 1/C_2\}. $$

Plugging (37) into the third equation of the system, and assuming $\lambda_2^* \neq 0$, yields

$$ \lambda_1^* = \frac{\delta(1 - S^*C_2)}{A_{12}} + \frac{\lambda_2^*}{A_{12}} \frac{1 - S^*C_2}{S^*C_2}. $$

By using a similar expression for $\lambda_2^*$, we obtain

$$ \lambda_1^* = \frac{\delta(1 - S^*C_2)}{A_{12}} + \frac{1 - S^*C_2}{A_{12}S^*C_2} \left( \frac{\delta(1 - S^*C_1)}{A_{21}} + \frac{\lambda_1^*}{A_{21}} \frac{1 - S^*C_1}{S^*C_1} \right). $$

This equation is equivalent to

$$ \left(1 - M_1 M_2 \frac{(1 - S^* C_1)(1 - S^* C_2)}{S^{*2}}\right) \lambda_1^* = \frac{\delta(1 - S^* C_2)}{A_{12}} \left(1 + M_2 \frac{(1 - S^* C_1)}{S^*}\right). $$

We obtain a similar equation for $\lambda_2^*$. Therefore, $\lambda_1^*$ and $\lambda_2^*$ can be uniquely determined from $S^*$. In addition, using the expression (37), we find $\lambda_{12}^*$ and $\lambda_{21}^*$ uniquely in terms of $S^*$. Since the right-hand sides of the last expressions are non-negative, and they are equal to 0 if and only if $S^*C_m = 1$, for $m \in \{1, 2\}$, we require

$$ 1 - M_1 M_2 \frac{(1 - S^* C_1)(1 - S^* C_2)}{S^{*2}} > 0 $$

to obtain a bi-endemic steady state. We may also see that if $S^* = 1/C_m$, then $\lambda_m^* = \lambda_{nm}^* = 0$, and the equilibrium is mono-endemic. Namely, we see that $\lambda_{12}^* \neq 0$.

Using the equation

$$ \Lambda = \delta S^* + \frac{\lambda_1^*}{C_1} + \frac{\lambda_2^*}{C_2}, $$

we get

$$ \Lambda = \delta S^* + \frac{\delta(1 - S^* C_2)S^* (S^* + M_2(1 - S^* C_1))}{C_1 A_{12}(S^{*2} - M_1 M_2(1 - S^* C_1)(1 - S^* C_2))} + \frac{\delta(1 - S^* C_1)S^* (S^* + M_1 M_2 A_{21}(1 - S^* C_2))}{C_2 A_{21}(S^{*2} - M_1 M_2(1 - S^* C_2)(1 - S^* C_1))}. $$
We may rearrange this equation to show that if $S^*$ is a disease-free or mono-endemic steady state, then it must be a (positive) root of the polynomial

$$P(x) = (x^2 - M_1 M_2 (1 - C_1 x)(1 - C_2 x)) \left(x - \frac{\Lambda}{\delta}\right) + x \frac{1 - C_2 x}{C_1 A_{12}} (x + M_2(1 - C_1 x)) + x \frac{1 - C_1 x}{C_2 A_{21}} (x + M_1(1 - C_2 x)).$$

The root will yield a bi-endemic equilibrium if and only if $0 < S^* < \min\{1/C_1, 1/C_2\}$ and $Q(S^*) > 0$, where

$$Q(x) = (1 - M_1 M_2 C_1 C_2)x^2 + M_1 M_2 (C_1 + C_2)x - M_1 M_2. \quad (38)$$

We have $Q(0) = -M_1 M_2 < 0$ and $Q(1/C_m) = 1/C_m^2 > 0$, $m \in \{1, 2\}$. Since $Q(x)$ is a quadratic polynomial, it must have exactly one root in the interval $(0, \min\{1/C_1, 1/C_2\})$, which we denote by $r_+$. For $x \in (0, r_+)$, $Q(x)$ is negative, and for $x \in (r_+, \min\{1/C_1, 1/C_2\})$, $Q(x)$ is positive. Thus, we search for the roots of $P(x)$ located in the interval $(r_+, \min\{1/C_1, 1/C_2\})$. $\square$

### 5.2.1. Symmetric case

Here, we assume that the infective power of the serotypes are similar, which means that

$$C_1 = C_2 =: C, C_{12} = C_{21} =: \hat{C}, \tilde{C}_1 = \tilde{C}_2 =: \tilde{C}, A_{12} = A_{21} =: A, M_1 = M_2 =: M. \quad (39)$$

Notice that we still may have $\beta_1 \neq \beta_2$, $\mu_1 \neq \mu_2$, $p_1 \neq p_2$, and so on. Moreover, even if $\lambda_1^* = \lambda_2^*$, we still may have $i_1^* \neq i_2^*$, namely through the distribution in $v_0$.

In this case, we have

$$MC = \frac{1}{A}, \quad \mathcal{R}_0 = \mathcal{R}_0^{(1)} = \mathcal{R}_0^{(2)} = \frac{\Lambda}{\delta}C \quad \text{and} \quad \tilde{C}\hat{C} = \frac{1}{M}.$$

A primary result can be obtained in the following theorem (see Figure 6). First, we define the function $\rho: [0, +\infty) \to [0, 1]$ by

$$\rho(x) = \begin{cases} \frac{x}{(1+x)^2} (3 - x + 2\sqrt{2\sqrt{1-x}}), & x < \frac{1}{2}, \\ 1, & x \ge \frac{1}{2}. \end{cases} \quad (40)$$

It is straightforward to verify that $\rho$ is an increasing and continuously differentiable function (see Figure 6).

**Theorem 5.6.** The model exhibits the following equilibrium behavior depending on $\mathcal{R}_0$ and MC:

(i) If $\mathcal{R}_0 < \rho(MC)$, no mono-endemic or bi-endemic equilibrium exists.
(ii) If $\rho(MC) < \mathcal{R}_0 < 1$, two bi-endemic equilibrium points exist, while no mono-endemic equilibrium is present. This scenario occurs exclusively when $MC < 1/2$, indicating the presence of a backward bifurcation.

(iii) If $\mathcal{R}_0 > 1$, there is a unique bi-endemic equilibrium and two mono-endemic equilibrium points.

**Remark 5.7.** • As we are in the case $\mathcal{R}_0 = \mathcal{R}_0^{(1)} = \mathcal{R}_0^{(2)}$, the two mono-endemic steady states, when they exist, are always unstable.

• If, in addition to the symmetric case, i.e., under the assumptions (39), we also assume that the recovery and transmission rates are constant, we find that

$$C = \frac{\beta_m}{\delta_m + \mu_m}, \quad \hat{C} = \frac{\beta_{mn}}{\delta_{mn} + \mu_{mn}} \quad \text{and} \quad \tilde{C} = \frac{\mu_m}{\delta_m + \mu_m}, \quad \text{for all } m, n \in \{1, 2\}.$$

In this case, we may write

$$2MC = \frac{2C}{\hat{C}\tilde{C}} = \frac{1}{\tilde{C}} \times \frac{\mathcal{R}_0^{(1)} + \mathcal{R}_0^{(2)}}{\mathcal{R}_0^{(12)}}.$$

The number $2MC$ may be interpreted as the ratio between the total expected number of new cases generated by two primary infections (each from a different serotype) to the expected number of new cases generated by a single secondary infection. Theorem 5.6 indicates that, if the secondarily infected individuals are more contagious than both primary serotypes combined, i.e., $2MC < 1$, then a phenomenon of backward bifurcation takes place (see Figures 6 and 7). We may also show that, in the case of $N$ serotypes, even with the same serotype reinfection, the condition for backward bifurcation of a steady state with $N$ coexisting endemic serotypes is

$$MC < 1 - \frac{1}{N}.$$

*Proof.* We will prove Theorem 5.6 by performing a bifurcation analysis on the parameter related to the demography of the susceptible individuals:

$$\alpha := \frac{\Lambda}{\delta}.$$

The bifurcation parameter $\alpha$ is independent of the infections, and in particular, independent of $C$ and $M$. In fact, we have

$$\mathcal{R}_0 = \alpha C.$$

The polynomial $P(x)$ may be rewritten as

$$P_{\alpha}(x) = (x - \alpha)Q(x) + 2MxR(x),$$
<figure><img src="image_7.png" /><figcaption>Figure 6: The different possible cases for the endemic steady states as a function of $x = MC$ in the symmetric case. The black curve corresponds to $\rho$ as a function of $x = MC$.</figcaption></figure>

where $Q(x)$ (given by Equation (38)) and $R(x)$ are defined as follows

$$Q(x) = x^2 - M^2(1 - Cx)^2 = (1 - M^2C^2)x^2 + 2M^2Cx - M^2,$$ $$R(x) = (1 - Cx)((1 - MC)x + M).$$

We recall that we are looking for the roots of $P_\alpha(x)$ located in $(r_+, 1/C)$. The polynomial $Q(x)$ has the two roots

$$r_+ := \frac{-M^2C + M}{1 - M^2C^2} = \frac{M}{1 + MC} \quad \text{and} \quad r_- := \frac{-M^2C - M}{1 - M^2C^2} = \frac{M}{1 - MC}.$$

Both roots are positive if $MC > 1$. In the case $MC < 1$, only one of these roots is positive. In fact, if $MC < 1$, we have

$$Cr_+ = \frac{MC}{1 + MC} < 1.$$

Thus, $r_- < 0 < r_+ < 1/C$. If $MC > 1$, we have $0 < r_+ < 1/C < r_-$. The roots of the polynomial $R(x)$ are $1/C$ and $r_-$. Thus, $R(x)$ and $Q(x)$ have a common root, $r_- \notin (0, 1/C)$. We may factor $(x - r_-)$ from $P_\alpha(x)$, yielding

$$P_{\alpha}(x) = (x - r_{-})[(1 - M^{2}C^{2})(x - \alpha)(x - r_{+}) + 2(1 - MC)Mx(1 - Cx)],$$

or equivalently,

$$P_{\alpha}(x) = (x - r_{-})(1 - MC)\tilde{P}_{\alpha}(x),$$

with

$$\tilde{P}_{\alpha}(x) := (1 + MC)(x - \alpha)(x - r_{+}) + 2Mx(1 - Cx).$$

We need only to study the roots $l_1(\alpha)$ and $l_2(\alpha)$ of the quadratic polynomial $\tilde{P}_\alpha(x)$. The functions $\alpha \mapsto l_i(\alpha)$, $i = 1, 2$, are holomorphic in $\mathbb{C}$ and may only leave the real line when they merge together. When one root is non-real, the
other must be its conjugate. Moreover, for every $\alpha > 0$, we have $\tilde{P}_{\alpha}(r_+) > 0$, so no root may cross $r_+$ through the real line. Finally, we have

$$ \tilde{P}_{\alpha}(1/C) = (1 + MC)(1/C - \alpha)(1/C - r_{+}), $$

so a root must cross $1/C$ exactly once, when $\alpha = 1/C$. We may also check that $\tilde{P}_0(x)$ has roots $l_1(0) = 0$ and $l_2(0) = r_-$.

Suppose $MC > 1$. Then, $l_2(0) = r_- > 1/C > r_+ > l_1(0) = 0$. So the roots cannot merge. The root $l_1(\alpha)$ must stay in $(-\infty, r_+)$ and $l_2(\alpha)$ must then cross $1/C$ towards the left when $\alpha = 1/C$. If $\alpha < 1/C$ (which is equivalent to $\mathcal{R}_0 < 1$), none of the three roots of $P_{\alpha}(x)$ is located in $(r_+, 1/C)$, so there is no bi-endemic equilibrium. If $\alpha > 1/C$, then $l_2(\alpha) \in (r_+, 1/C)$, so there is exactly one bi-endemic equilibrium.

Suppose $MC < 1$. Then, the roots start at $l_2(0) < l_1(0) < r_+$. Since a root must cross $1/C$ when $\alpha = 1/C$, and the roots cannot cross $r_+$ through the real line, they must merge before becoming complex conjugates and then merge again to reappear in the interval $(r_+, \infty)$. To know whether they appear in $(r_+, 1/C)$ or in $(1/C, +\infty)$, we calculate the sign of $l'_i(\alpha)$ (i.e., the direction through which it will cross) for $i = 1, 2$. By differentiating the function $\alpha \mapsto \tilde{P}_{\alpha}(l_i(\alpha)) = 0$, and taking $\alpha = 1/C$, we have $l_i(1/C) = 1/C$, which yields:

$$ l'_{i}(1/C) = (1 + MC) \frac{1/C - r_{+}}{\tilde{P}'_{1/C}(1/C)}. $$

The sign of $l'_i(1/C)$ is the same as the sign of $\tilde{P}'_{1/C}(1/C)$. We can check that $\tilde{P}'_{1/C}(1/C) = \frac{1}{C} - 2M$ is negative if and only if $MC > 1/2$. In this case, the roots must merge and reappear in the interval $(1/C, +\infty)$. Hence, the system has no bi-endemic steady states if $\alpha < 1/C$, because $l_i(\alpha) > l_i(1/C) = 1/C$ for $i = 1, 2$, and exactly one bi-endemic equilibrium if $\alpha > 1/C$, because one of the roots has crossed $1/C$ and lies in $(r_+, 1/C)$.

If $MC < 1/2$, then $l'_i(1/C) > 0$, so the roots merge and reappear in the interval $(r_+, 1/C)$. The roots reappear when $\alpha$ exceeds a threshold noted $\alpha_0$. Therefore, there is initially no bi-endemic equilibrium, then two appear when $\alpha > \alpha_0$, and exactly one appears when $\alpha > 1/C$. We can calculate $\alpha_0$ by analyzing the discriminant of $\tilde{P}_{\alpha}(x)$:

$$ \Delta(\alpha) = (1 + MC)^2\alpha^2 + 2M(MC - 3)\alpha + M^2. $$

The double root of $\tilde{P}_{\alpha}(x)$ appears when $\Delta(\alpha) = 0$ for the second time. Indeed, at first, both roots are real ($\Delta(0) = M^2 > 0$), then they merge ($\Delta(\alpha) = 0$), and they rejoin the real line when $\Delta(\alpha) = 0$ again. Finally, $\Delta(\alpha)$ is, in turn, a quadratic polynomial in $\alpha$, and its larger root is

$$ \alpha_0 = M \times \frac{3 - MC + 2\sqrt{2\sqrt{1 - MC}}}{(1 + MC)^2}. \qquad (41) $$

Moreover, $\alpha = \Lambda/\delta > \alpha_0$ is equivalent to $\mathcal{R}_0 > C\alpha_0 = \rho(MC)$, where $\rho$ is given by Equation (40). The result is thus proven. $\square$
### 5.2.2. Non symmetric case

The general case, without symmetry, yields the same types of bifurcation diagrams, but with less neat and interpretable conditions than previously. For the remainder of this section, we assume that serotype 1 is more infectious than serotype 2, i.e.,

$$ \mathcal{R}_0^{(1)} > \mathcal{R}_0^{(2)}. \qquad (42) $$

This is equivalent to say that $C_1 > C_2$. We could assume the opposite without loss of generality. We define the polynomial

$$ L(x) = (1 + M_1 M_2 C_1 C_2 - M_1 C_2 - M_2 C_1)x^2 \qquad (43) $$

$$ + M_1 M_2 (C_1 - C_2)(M_1 C_2 - 1)x - M_1 M_2 C_1 \alpha^*, \qquad (44) $$

where

$$ \alpha^* := \frac{1}{C_1} + M_1 C_2 \left( \frac{1}{C_2} - \frac{1}{C_1} \right) > 0 \quad \text{and} \quad M_m := (C_m A_{mn})^{-1}. \qquad (45) $$

**Theorem 5.8.** Assume that the condition (42) is satisfied. This is also equivalent to $C_1 > C_2$. Then, the mono-endemic equilibrium for serotype 1 is locally asymptotically stable if and only if

$$ \frac{1}{C_1} < \frac{\Lambda}{\delta} < \alpha^*. $$

Moreover, exactly one of the following two bifurcation profiles can occur:

* If $L(1/C_1) < 0$, where L is defined by (43), then:

- No bi-endemic steady state exists if $\Lambda/\delta < \alpha^*$.

- A unique bi-endemic steady state exists if $\Lambda/\delta > \alpha^*$.

* If $L(1/C_1) > 0$, then a backward bifurcation occurs. There exists $0 \le \rho < \alpha^*$, where $\alpha^*$ is given by (45), such that:

- No bi-endemic equilibrium points exist when $\Lambda/\delta < \rho$.

- Two bi-endemic equilibrium points exist when $\rho < \Lambda/\delta < \alpha^*$.

- A unique bi-endemic equilibrium point exists when $\Lambda/\delta > \alpha^*$.

*Proof.* We consider the polynomial previously used, $P_\alpha(x)$, which can be written as

$$ P_{\alpha}(x) = (x - \alpha)Q(x) + xM_1(1 - C_2x)(x + M_2(1 - C_1x)) \\ + xM_2(1 - C_1x)(x + M_1(1 - C_2x)). $$

We recall that $Q(x) = (1 - M_1 M_2 C_1 C_2)x^2 + M_1 M_2(C_1 + C_2)x - M_1 M_2$ is given by (38), and has roots $r_+$ and $r_-$ as seen previously. We have $P_\alpha(0) = M_1 M_2 \alpha > 0$ and

$$ P_{\alpha} \left( \frac{1}{C_1} \right) = \frac{1}{C_1^2} \left( \frac{1}{C_1} - \alpha + M_1 C_2 \left( \frac{1}{C_2} - \frac{1}{C_1} \right) \right). $$
We know that the complex roots $r_0(\alpha)$, $r_1(\alpha)$ and $r_2(\alpha)$ of $P_\alpha$ are analytic functions in $\mathbb{C}$, verifying the property that a root can only leave the real line $\mathbb{R}$ if it merges with another root, and if one of them is non-real, then the other must be its conjugate. Moreover, we see that there exists

$$ \alpha^* := \frac{1}{C_1} + M_1 C_2 \left( \frac{1}{C_2} - \frac{1}{C_1} \right) > \frac{1}{C_1} $$

such that $P_\alpha(1/C_1)$ is positive for $\alpha < \alpha^*$ and negative for $\alpha > \alpha^*$. Thus, a root must cross $1/C_1$ through the real line at $\alpha^*$. Finally, a root cannot cross $r_+$, since $P_\alpha(r_+) > 0$. Let us now consider $P_0(x) = x(q_2x^2 + q_1x + q_0)$, with

$$ q_2 = 1 + M_1 M_2 C_1 C_2 - M_1 C_2 - M_2 C_1, \quad q_1 = M_1 + M_2 - M_1 M_2 (C_1 + C_2), \quad q_0 = M_1 M_2. $$

We may check that $P_0(r_+) > 0$ and $P_0(1/C_1) > 0$. Moreover we have $P_0(0) = 0$ and $q_0 > 0$. Consider first the case $q_2 < 0$. Then, $P_0$ has one negative root, one at 0, and one in $(1/C_1, +\infty)$. We distinguish between the cases $\alpha < \alpha^*$ and $\alpha > \alpha^*$. We first prove that a double root of $P_\alpha$ cannot appear in the interval $(r_+, 1/C_1)$ if $\alpha < \alpha^*$. Suppose, by contradiction, that such a double root exists. Since the mapping $\alpha \mapsto P_\alpha(x)$ is decreasing for every $x \in (r_+, 1/C_1)$, it follows that $P_{\alpha^*}$ would have at least one root in this interval. Therefore, it is sufficient to prove that $P_{\alpha^*}$ has no roots in $(r_+, 1/C_1)$. To do so, we factorize $P_{\alpha^*}$ as

$$ P_{\alpha^*}(x) = (x - 1/C_1)(l_2x^2 + l_1x + l_0) := (x - 1/C_1)L(x), $$

where the coefficients of $L$ are given by

$$ l_2 = q_2 < 0, \quad l_1 = M_1 M_2 (C_1 - C_2)(M_1 C_2 - 1), \quad l_0 = -C_1 \alpha^* M_1 M_2 < 0. $$

The polynomial $L$ attains its maximum at $-l_1/(2l_2)$. We denote by $g$ the function, defined for $C_1 \ge C_2$, that associates the argmax with $C_1$:

$$ g(C_1) := -\frac{l_1}{2l_2} = \frac{M_1 M_2 (C_1 - C_2)(1 - M_1 C_2)}{2(1 + M_1 M_2 C_1 C_2 - M_1 C_2 - M_2 C_1)}. $$

As a ratio of two affine functions, $g$ is monotonic and we have $g(C_2) = 0$. Furthermore,

$$ \lim_{C_1 \to +\infty} g(C_1) = \frac{M_1 M_2 (1 - M_1 C_2)}{2 M_2 (M_1 C_2 - 1)} = -\frac{M_1}{2} < 0. $$

Thus, $g$ is non-positive and the maximum of $L$ is attained in $\mathbb{R}_-$. Since $L(0) = l_0 < 0$, this means that $L$ and thus $P_{\alpha^*}$, has either two complex roots or two negative real roots. Therefore, $P_{\alpha^*}$ has no root in $(r_+, 1/C_1)$. Hence, if $\alpha < \alpha^*$, there is no bi-endemic equilibrium. For $\alpha > \alpha^*$, we will show that a double root cannot appear in $(r_+, 1/C_1)$. To do so, we use Descartes' rule of signs. The coefficients of $P_\alpha$ are:

$$ a_3 := q_2 < 0, \quad a_2(\alpha) := q_1 - \alpha(1 - M_1 M_2 C_1 C_2), $$
$$a_1(\alpha) := M_1 M_2(1 - \alpha(C_1 + C_2)), \text{ and } a_0(\alpha) := \alpha M_1 M_2 > 0.$$

Since $\alpha > \alpha^* > 1/C_1 > 1/(C_1 + C_2)$, we have $a_1(\alpha) < 0$. Thus, $P_\alpha$ can have three positive roots only if $a_2(\alpha) > 0$. However, $a_2(\alpha)$ is an affine function of $\alpha$, so it can only change signs once. Therefore, if two roots appear on the positive real line, they cannot leave it. Now, the equation $P_\alpha(x) = 0$ can be rewritten as $\alpha = P_0(x)/Q(x)$. If $Q(x)$ and $P_0(x)$ have the same sign, then there exists a unique $\alpha$ such that $P_\alpha(x) = 0$. Since $P_0(x)$ is cubic and $Q(x)$ is quadratic, there exists an infinite interval, either $(-\infty, \kappa)$ or $(\kappa, +\infty)$, where they both share the same sign. Thus, there must be a root that explores the entire interval, which implies that one of the roots is unbounded. Moreover, using the continuity of $\alpha \mapsto r(\alpha)$, we may conclude that a root satisfies $r(\alpha) \to \pm\infty$ as $\alpha \to +\infty$. Suppose a double root appears in $(r_+, 1/C_1)$ at $\alpha > \alpha^*$: we know that there was already one root (that appeared through $1/C_1$) in the interval, and thus, all three roots must be in this interval. They cannot merge (by Descartes' rule of signs), and they cannot cross $r_+$ or $1/C_1$; they must stay bounded. This leads to a contradiction with the fact that $|r(\alpha)| \to +\infty$ as $\alpha \to +\infty$. Therefore, there can only be one bi-endemic equilibrium for $\alpha > \alpha^*$.

Suppose $q_2 > 0$. Then, $P_\alpha$ always has a negative root (it has a positive leading coefficient and is positive at 0). Moreover, since $\alpha \mapsto P_\alpha(x)$ is decreasing for every $x \in (r_+, 1/C_1)$, the number of roots in $(r_+, 1/C_1)$ can only go from 0 to 2, or from 0 to 1, with the last switch occurring when $\alpha = \alpha^*$. The first case happens when the root $r_i$ passing through $1/C_1$ at $\alpha = \alpha^*$ satisfies $r'_i(\alpha^*) > 0$, or equivalently, $L(1/C_1) > 0$. Thus, in this case, there exists $0 \le \rho < \alpha^*$ such that no bi-endemic equilibrium exists for $\alpha < \rho$. The two roots merge in $(r_+, 1/C_1)$ when $\alpha = \rho$, resulting in two bi-endemic equilibrium points for $\rho < \alpha < \alpha^*$. Finally, one of the roots exits the interval $(r_+, 1/C_1)$ when $\alpha > \alpha^*$, leaving only one bi-endemic equilibrium. In the other case, where $L(1/C_1) < 0$, there is no bi-endemic equilibrium for $\alpha > \alpha^*$ and only one for $\alpha < \alpha^*$. Thus, the only possible bifurcation diagrams are the ones listed previously. In conclusion, we have proven that only one of the two types of bifurcation diagrams presented can occur. $\square$

Figure 8 provides a clearer view of the possible bifurcation profiles listed in Theorem 5.8.

**Remark 5.9.** • The condition $\alpha > \alpha^*$ in Theorem 5.8, may be rewritten as

$$\mathcal{R}_0^{(1)} - 1 > \frac{C_1 - C_2}{\tilde{C}_1 C_{12}} = \frac{1}{\tilde{C}_1} \frac{\mathcal{R}_0^{(1)} - \mathcal{R}_0^{(2)}}{\mathcal{R}_0^{(12)}}.$$

It is the same condition for the stability of the mono-endemic equilibrium. In Theorem 5.8, choosing $\mathcal{R}_0^{(1)}$ as the bifurcation parameter, as might initially be envisaged, would not be appropriate, since $\mathcal{R}_0^{(1)}$ and $C_1$ are not independent. Varying $\mathcal{R}_0^{(1)}$ would modify other parameters that also depend on $C_1$, making it difficult to determine the occurrence of steady
<figure><img src="image_8.png" /><figcaption>Figure 7: Bifurcation diagrams of $S^*$ for two different values of $MC$ in the symmetric case. The bifurcation parameter is $\alpha = \Lambda/\delta$. Green lines represent the disease-free steady state, orange lines correspond to the combined mono-endemic equilibrium points, and red lines denote the bi-endemic equilibrium points. Note that the two mono-endemic steady states have the same value of $S^*$ in the symmetric case, but have distinct components. Dashed lines indicate unstable steady states, while solid lines represent locally asymptotically stable steady states. Left: $MC = 0.15 < 1/2$. Right: $MC = 0.75 > 1/2$. The vertical dash-dotted line marks $\alpha = \alpha_0$, where $\alpha_0$ is given by (41), and the vertical dotted lines indicate $\alpha = 1/C$, i.e., $\mathcal{R}_0 = 1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: The two possible bifurcation diagrams for $S^*$, when we vary $\alpha = \Lambda/\delta$ in the non symmetric case, $\mathcal{R}_0^{(1)} > \mathcal{R}_0^{(2)}$. Note that, in the non symmetric case, a given value of $S^*$ leads to a unique equilibrium. The vertical dash-dotted line corresponds to $\alpha = \rho$. The first vertical dotted line is $1/C_1$ (i.e., $\mathcal{R}_0^{(1)} = 1$); the second represents $\alpha^*$; the last one corresponds to $1/C_2$. Notice that we could have $1/C_2 < \alpha^*$ and so the existence of five steady states. The parameters are: $M_1 = 0.25$, $C_1 = 1.2$, $C_2 = 0.6$ and $M_2 = 0.7$ (left) or $M_2 = 0.1$ (right).</figcaption></figure>

states. The parameter $\Lambda/\delta$, which is independent of $C_1$, provides a simpler and more appropriate choice.

* We can see that a bi-endemic equilibrium exists if $\Lambda/\delta > \alpha^*$. This is not always the case when there are multiple serotypes in competition, as shown by [25], i.e., in a model without age structure. Reinfection allows serotypes to coexist, since the same resources (represented here by
the susceptible class) can be used by two different serotypes. The inva-
sion reproduction number of serotype $n$ when it invades a population at
mono-endemic equilibrium with serotype $m$, denoted $\mathcal{R}_n^{\text{inv}}$, is the effective
reproduction number of serotype $n$ when it is introduced in a population
where serotype $m$ is already established at its endemic equilibrium. This
number is defined by

$$
\mathcal{R}_n^{\text{inv}} := \frac{\Lambda}{\delta \alpha^*} = \frac{\mathcal{R}_0^{(m)} - 1}{\mathcal{R}_0^{(m)} - \mathcal{R}_0^{(n)}} \tilde{C}_m \mathcal{R}_0^{(mn)}. \quad (46)
$$

If $\mathcal{R}_n^{\text{inv}} > 1$, then serotype $n$ "invades" the mono-endemic equilibrium of serotype $m$. This number may also be rewritten as follows

$$
\mathcal{R}_n^{\text{inv}} = \bar{S}_m^* C_{mn} \left( 1 + \frac{C_n}{C_m - C_n} \right),
$$

where $\bar{S}_m^*$ is the number of individuals who recovered from serotype $m$, associated to the mono-endemic equilibrium for serotype $m$. This is similar to a classic basic reproduction number, but where the susceptible population considered is $\bar{S}_m^*$ instead of $S^*$.

* The second bifurcation diagram (on the right of Figure 8) shows a backward bifurcation. Previous studies, such as that presented in [12], have also shown that this phenomenon occurs when the epidemiological model takes reinfection into account. Although our model does not explicitly include reinfection by the same serotype, but rather cross-infection between different serotypes, we still observe this backward bifurcation, particularly for bi-endemic equilibrium points. If we were to consider reinfection by the same serotype, we could imagine a backward bifurcation for mono-endemic equilibrium points. In this case, the constant $\rho$, defined in Theorem 5.8, could be interpreted, as in the article [15], like a threshold of reinfection.

The backward bifurcation condition, $L(1/C_1) > 0$, in Theorem 5.8, is par-
ticularly interesting to analyze. After a simple calculation of $L(1/C_1)$, we see
that a backward bifurcation occurs if and only if

$$
\sqrt{M_1 M_2 C_1} < \frac{1}{2} \quad \text{and} \quad C_2 \in (0, C_1) \cap (c_-, c_+),
$$

with

$$
c_{-} = \frac{2M_{2}C_{1} + 2M_{1}M_{2}C_{1}^{2} - 1 - \sqrt{1 - 4M_{1}M_{2}C_{1}^{2}}}{2M_{1}M_{2}C_{1}},
$$

$$
\text{and } c_+ = \frac{2M_2 C_1 + 2M_1 M_2 C_1^2 - 1 + \sqrt{1 - 4M_1 M_2 C_1^2}}{2M_1 M_2 C_1}.
$$

We then have the following proposition, which presents some simple cases of the
presence or absence of a backward bifurcation.
**Proposition 5.10.** Assume that the condition (42) is satisfied (this is equivalent to $C_1 > C_2$), and let $x_1 = M_1C_1$ and $x_2 = M_2C_1$.

1. If

$$ \max\{x_1 + x_2, x_2(1 + x_1)^2\} < 1, \quad (47) $$

then there is a backward bifurcation for any $C_2 \in (0, C_1)$.

2. If

$$ 4x_1x_2 > 1, \quad (48) $$

then there is no backward bifurcation for any $C_2 \in (0, C_1)$.

*Proof.* We define, for $C_2 \in [0, C_1]$,

$$ \mathcal{L}(C_2) := 1 + M_1 M_2 C_1 C_2 - M_1 C_2 - M_2 C_1 + C_1 M_1 M_2 (C_1 - C_2)(M_1 C_2 - 1) - M_1 M_2 C_1^3 \alpha^*. $$

This corresponds to $\mathcal{L}(C_2) = C_1^2 L(1/C_1)$. The function $C_2 \mapsto \mathcal{L}(C_2)$ is a quadratic polynomial, with a negative quadratic coefficient. Moreover,

$$ \mathcal{L}(0) = 1 - x_2(1 + x_1)^2, \quad \mathcal{L}(C_1) = 1 - x_1 - x_2, $$

and the discriminant of $\mathcal{L}(C_2)$ is

$$ \Delta = M_1^2(1 - 4x_1x_2). $$

Then, the sign of the polynomial $\mathcal{L}(C_2)$, for $C_2 \in (0, C_1)$, is directly determined, and the proposition is therefore proved. $\square$

**Remark 5.11.** Note that the conditions (47) and (48) are already included in Theorem 5.6. Indeed, the condition $x_1 + x_2 < 1$ can be rewritten as

$$ \frac{M_1 + M_2}{2} C_1 < \frac{1}{2}, $$

while the condition (48) becomes

$$ \sqrt{M_1 M_2} C_1 > \frac{1}{2}. $$

These expressions are similar to the condition on $MC$ in Theorem 5.6, now involving the arithmetic and geometric means of $M_1$ and $M_2$. This confirms once again that the existence of a backward bifurcation is directly related to a high average infectious potential of secondary infections, when compared to primary infections.
<figure><img src="image_10.png" /><figcaption>Figure 9: Zones where backward bifurcation occurs: (i) for every $C_2 \in (0, C_1)$, (ii) for every $C_2 \in (c_-, C_1)$, (iii) for every $C_2 \in (0, c_+)$, and (iv) for no $C_2 \in (0, C_1)$, in the $(x_1 = M_1C_1, x_2 = M_2C_1)$-plane. In the grey zones, either there is backward bifurcation for $C_2 \in (c_-, c_+) \subset (0, C_1)$, or there is no backward bifurcation for any value of $C_2 \in (0, C_1)$.</figcaption></figure>

# 6. Relationship between immunological and epidemiological dynamics

Examining the relationship between immunological and epidemiological parameters can help in developing treatment or vaccine strategies to control an epidemic. The age structure of the model plays a crucial role, as it links the full immunological dynamics and parameters to key epidemiological thresholds, such as the basic reproduction number and the invasion reproduction number. The previous analysis of steady states indicates that the system is mainly characterized by three key parameters: $R_0 = \max\{\mathcal{R}_0^{(1)}, \mathcal{R}_0^{(2)}\}$, $\mathcal{R}_2^{\text{inv}}$ given by (46) and $\rho$ defined by (40) (as well as that defined in Theorem 5.8). The basic reproduction number $R_0$ determines the stability of the trivial equilibrium and the existence of mono-endemic equilibrium points. It depends only on the immunological parameters of the primary infection caused by a single serotype. In what follows, we assume that $C_1 \ge C_2$, which is equivalent to $\mathcal{R}_0^{(1)} \ge \mathcal{R}_0^{(2)}$. So, $R_0 = \mathcal{R}_0^{(1)}$.

## 6.1. Relationship between $R_0$ and immunological parameters

We consider that the non-trivial immunological equilibrium $(v_1^*, a_1^*)$ for serotype 1, System (7), exists and that a newly infected individual begins with an initial condition equal to this equilibrium. More precisely, this implies that $p_1(v_0) = \delta_{v_1^*}$, where $p_1$ is the probability defined by (H1), for serotype 1, and $\delta$ represents the Dirac-delta measure. In simpler terms, we assume that during the primary infection with serotype 1, the immunological dynamics remain constant and are equal to the non-trivial equilibrium: $v_1(\tau) = v_1^*$ and $a_1(\tau) = a_1^*$,
for all $\tau$. Recall that System (7) is defined as

$$ egin{cases} \frac{dv}{d\tau} = r(a)v\left(1 - \frac{v}{K}\right) - \gamma_v v, \\ \frac{da}{d\tau} = \eta(v) - \gamma_a a. \end{cases} $$

For the simplicity, we use the notations $\gamma_a$, $\gamma_v$, $K$, $\eta$, $r$, for serotype 1. Under the previous assumption, $\beta_1$ and $\mu_1$ depend only on $(v_1^*, a_1^*)$ and the basic reproduction number becomes

$$ \mathcal{R}_0 = \frac{\Lambda}{\delta} C_1 = \frac{\Lambda}{\delta} \int_0^{\infty} \beta_1(v_1^*, a_1^*) \pi_1(\tau, v_1^*) d\tau = \frac{\Lambda}{\delta} \times \frac{\beta_1(v_1^*, a_1^*)}{\delta_1 + \mu_1(v_1^*, a_1^*)} $$

We know from Theorem 2.1 that $a_1^* = \eta(v_1^*)/\gamma_a$. Then, we may write

$$ \bar{\beta}_1(v_1^*, \gamma_a) := \beta_1 \left( v_1^*, \frac{\eta(v_1^*)}{\gamma_a} \right) \quad \text{and} \quad \bar{\mu}_1(v_1^*, \gamma_a) := \mu_1 \left( v_1^*, \frac{\eta(v_1^*)}{\gamma_a} \right). $$

Equivalently, we can express $\beta_1$ and $\mu_1$ as functions of $(a_1^*, \gamma_a)$. So,

$$ \mathcal{R}_0(v^*, \gamma_a) = \frac{\Lambda}{\delta} \times \frac{\bar{\beta}_1(v_1^*, \gamma_a)}{\delta_1 + \bar{\mu}_1(v_1^*, \gamma_a)}. $$

The following result shows the variation of $\mathcal{R}_0$ as a function of the immunological parameters $K$, $\gamma_a$ and $\gamma_v$ of serotype 1.

**Proposition 6.1.** 1. Assume that the force of infection $\bar{\beta}_1$ and the recovery rate $\bar{\mu}_1$ are increasing and decreasing with respect to the viral load, respectively. Then, $\mathcal{R}_0$ is increasing with respect to $K$ and decreasing with respect to $\gamma_v$.

2. Furthermore, if $\bar{\beta}_1(v_1^*, \gamma_a) = \bar{\beta}_1(v_1^*)$ and $\bar{\mu}_1(v_1^*, \gamma_a) = \bar{\mu}_1(v_1^*)$, then $\mathcal{R}_0$ is increasing with respect to $\gamma_a$.

*Proof.* Let $\psi_1 \in \{K, \gamma_a, \gamma_b\}$. If $\psi_1 \neq \gamma_a$, then we have

$$ \frac{\partial \mathcal{R}_0}{\partial \psi_1} = \frac{\Lambda}{\delta} \left( \frac{\frac{\partial \bar{\beta}_1}{\partial v}(v_1^*, \gamma_a)(\delta_1 + \bar{\mu}_1(v_1^*, \gamma_a)) - \bar{\beta}_1(v_1^*, \gamma_a) \frac{\partial \bar{\mu}_1}{\partial v}(v_1^*, \gamma_a)}{(\delta_1 + \bar{\mu}_1(v_1^*, \gamma_a))^2} \right) \frac{\partial v_1^*}{\partial \psi_1}. $$

Using the assumptions on $\bar{\beta}_1$ and $\bar{\mu}_1$, we get that $\partial \mathcal{R}_0/\partial \psi_1$ and $\partial v_1^*/\partial \psi_1$ have the same sign. To compute $\partial v_1^*/\partial \psi_1$, with $\psi_1 \in \{K, \gamma_a, \gamma_b\}$, we recall that

$$ L_1(v_1^*, \psi_1) = \frac{1}{\gamma_v} r \left( \frac{1}{\gamma_a} \eta(v_1^*) \right) \left( 1 - \frac{v_1^*}{K} \right) = 1 \quad \text{with} \quad a_1^* = \frac{1}{\gamma_a} \eta(v_1^*). $$

Then, by differentiating, we obtain

$$ \frac{\partial v_1^*}{\partial \psi_1} = - \frac{\partial L_1}{\partial \psi_1} / \frac{\partial L_1}{\partial v_1^*}. $$
It is clear that

$$ \frac{\partial L_1}{\partial v_1^*} = \frac{1}{\gamma_v \gamma_a} \eta'(v_1^*) r' \left( \frac{1}{\gamma_a} p(v_1^*) \right) \left( 1 - \frac{v_1^*}{K} \right) - \frac{1}{K \gamma_v} r \left( \frac{1}{\gamma_a} \eta(v_1^*) \right) < 0. $$

Moreover,

$$ \frac{\partial L_1}{\partial \gamma_v}(v_1^*, \gamma_v) < 0, \quad \frac{\partial L_1}{\partial \gamma_a}(v_1^*, \gamma_a) > 0 \quad \text{and} \quad \frac{\partial L_1}{\partial K}(v_1^*, K) > 0. $$

We conclude that

$$ \frac{\partial v_1^*}{\partial \gamma_v} < 0, \quad \frac{\partial v_1^*}{\partial \gamma_a} > 0 \quad \text{and} \quad \frac{\partial v_1^*}{\partial K} > 0. $$

Then, $\mathcal{R}_0$ is decreasing with respect to the virus mortality rate $\gamma_v$ and increasing with the virus carrying capacity $K$. With regard to the antibody clearance $\gamma_a$, there is no immediate conclusion. However, if the functions $\tilde{\beta}_1$ and $\tilde{\mu}_1$ depend only on $v_1^*$, we can use the same calculations to show that $\mathcal{R}_0$ is increasing with respect to $\gamma_a$. $\square$

## 6.2. Relationship between the invasion reproduction number and the immunological parameters

Recall that we defined the invasion reproduction number $\mathcal{R}_2^{\text{inv}}$ in (46): serotype 2 invades serotype 1, i.e., the mono-endemic equilibrium, associated with serotype 1, destabilizes and the bi-endemic equilibrium becomes the only stable steady state if and only if $\mathcal{R}_2^{\text{inv}} > 1$. We can see from the expression

$$ \mathcal{R}_2^{\text{inv}} = \frac{\mathcal{R}_0^{(1)} - 1}{\mathcal{R}_0^{(1)} - \mathcal{R}_0^{(2)}} \tilde{C}_1 \mathcal{R}_0^{(12)}, $$

that $\mathcal{R}_2^{\text{inv}}$ depends on the immunological parameters of the first infection with serotype 1 through $\mathcal{R}_0^{(1)}$ and $\tilde{C}_1$, the immunological parameters of the first infection with serotype 2 through $\mathcal{R}_0^{(2)}$, and the immunological parameters of the second infection with serotype 2 after a first infection with serotype 1 through $\mathcal{R}_0^{(12)}$. If we assume that all transmission and recovery rates depend only on the viral load at non-trivial immunological equilibrium, i.e., $\beta_m(v_m^*, a_m^*) = \beta_m(v_m^*)$ and $\mu_m(v_m^*, a_m^*) = \mu_m(v_m^*)$, then $\mathcal{R}_2^{\text{inv}}$ may be rewritten as a function of $v_1^*, v_2^*, v_{12}^*$ using, for $m \in \{1, 2\}$, the expressions

$$ \mathcal{R}_0^{(m)} = \frac{\Lambda}{\delta} \times \frac{\beta_m(v_m^*)}{\delta_m + \mu_m(v_m^*)}, \quad \mathcal{R}_0^{(12)} = \frac{\Lambda}{\delta} \times \frac{\beta_{12}(v_{12}^*)}{\delta_{12} + \mu_{12}(v_{12}^*)} \quad \text{and} \quad \tilde{C}_1 = \frac{\mu_1(v_1^*)}{\delta_1 + \mu_1(v_1^*)}. $$

We refer here to $\gamma_v^m, \gamma_a^m, K_m, m \in \{1, 2\}$, as the immunological parameters associated with a primary infection by the serotype $m$, and $\gamma_v^{mn}, \gamma_a^{mn}, \gamma_b^{mn}, K_{mn}, m \neq n \in \{1, 2\}$, as those corresponding to a secondary infection by the serotype $n$ following a primary infection by the serotype $m$.
**Proposition 6.2.** We assume that all the functions $\beta$ and $\mu$ depend only on the viral load. Suppose also that the immunological trajectories $(v_{12}(\tau), a_{12}(\tau), b_{12}(\tau))$ start at a non-trivial equilibrium $(v_{12}^*, a_{12}^*, b_{12}^*)$, and thus are constant. Moreover, assume that $v_m \mapsto \beta_m(v_m)$ and $v_{mn} \mapsto \beta_{mn}(v_{mn})$ are increasing, and $v_m \mapsto \mu_m(v_m)$ and $v_{nm} \mapsto \mu_{nm}(v_{nm})$ are decreasing. Then, there exists a threshold $\kappa_2 < 1$ such that:

(i) $\mathcal{R}_2^{\text{inv}}$ is increasing with respect to $K_2$ and $\gamma_a^2$, and decreasing with respect to $\gamma_v^2$.

(ii) If $\mathcal{R}_0^{(2)} < \kappa_2$, then $\mathcal{R}_2^{\text{inv}}$ is increasing with respect to $K_1$ and $\gamma_a^1$, and decreasing with respect to $\gamma_v^1$.

(iii) If $\mathcal{R}_0^{(2)} > \kappa_2$, then $\mathcal{R}_2^{\text{inv}}$ is decreasing with respect to $K_1$ and $\gamma_a^1$, and increasing with respect to $\gamma_v^1$.

(iv) If $\partial L_2(v_{12}^*, \psi_{12}) / \partial v_{12}^* < 0$, where $L_2$ is given by (4), then $\mathcal{R}_2^{\text{inv}}$ is increasing with respect to $K_{12}$ and decreasing with respect to $\gamma_v^{12}$.

*Proof.* (i). Clearly, the variations of $\mathcal{R}_2^{\text{inv}}$ as a function of the parameters $K_2$, $\gamma_a^2$ and $\gamma_v^2$ are the same as those of $\mathcal{R}_0^{(2)}$ (the same as for $\mathcal{R}_0$ in Proposition 6.1). This proves the result stated.

For (ii) and (iii), we have, for a parameter $\psi_1$ associated with a primary infection of serotype 1,

$$ \frac{\partial \mathcal{R}_2^{\text{inv}}}{\partial \psi_1}(\psi_1, \psi_2, \psi_{12}) = \mathcal{R}_0^{(12)} \frac{(\mathcal{R}_0^{(1)} - 1) \tilde{C}'_1 (\mathcal{R}_0^{(1)} - \mathcal{R}_0^{(2)}) + \tilde{C}_1 \mathcal{R}_0^{(1)'} (1 - \mathcal{R}_0^{(2)})}{(\mathcal{R}_0^{(1)} - \mathcal{R}_0^{(2)})^2}, $$

with $\mathcal{R}_0^{(1)'} = \partial \mathcal{R}_0^{(1)} / \partial \psi_1$ and $\tilde{C}'_1 = \partial \tilde{C}_1 / \partial \psi_1$. Using the same method as in the proof of Proposition 6.1 for $\mathcal{R}_0$, we show that $\partial \tilde{C}_1 / \partial \psi_1$ has a sign opposite to the one of $\partial v_1^* / \partial \psi_1$:

$$ \frac{\partial \tilde{C}_1}{\partial \psi_1} = \frac{\delta_1}{(\delta_1 + \mu_1(v_1^*))^2} \frac{d\mu_1}{dv} \frac{\partial v_1^*}{\partial \psi_1}. $$

If $\mathcal{R}_0^{(2)}$ is larger than a certain threshold less than 1, noted $\kappa_2$, we can see that $\partial \mathcal{R}_2^{\text{inv}} / \partial \psi_1$ and $\partial v_1^* / \partial \psi_1$ are of opposite signs. This proves (ii) and (iii).

To prove (iv), we can again see that $\partial \mathcal{R}_0^{(12)} / \partial \psi_{12}$ has the same sign as $\partial v_{12}^* / \partial \psi_{12}$, and, using the relation $L_2(v_{12}^*, \psi_{12}) = 0$, with

$$ L_2(v_{12}^*, \psi_{12}) = r \left( \frac{\eta_{12}(v_{12}^*)}{\gamma_a^{12}}, \frac{\bar{\eta}_{12}(v_{12}^*)}{\gamma_b^{12}} \right) \left( 1 - \frac{v_{12}^*}{K_{12}} \right) - \gamma_v^{12}, $$

we find

$$ \frac{\partial v_{12}^*}{\partial \psi_{12}} = - \frac{\partial L_2}{\partial \psi_{12}} / \frac{\partial L_2}{\partial v_{12}^*}. $$
If $\partial L_2(v_{12}^*, \psi_{12})/\partial v_{12}^* < 0$, then $\partial v_{12}^*/\partial \psi_{12}$ has the same sign as $\partial L_2/\partial \psi_{12}$. This
yields the same results for parameters $K_{12}$ and $\gamma_v^{12}$ as obtained in Proposition
6.1. For the parameters $\gamma_a^{12}$ and $\gamma_b^{12}$, $R_0^{(12)}$ has more complex variations, since
the function $r(a_{12}, b_{12})$ is not monotonic with respect to $a_{12}$ or $b_{12}$. $\square$

**Remark 6.3.** If the immunological basic reproduction number of the second infection, $R_{02}$, defined by (3), is greater than 1 then, by Remark 2.6, there is at least one non-trivial equilibrium $(v_{12}^*, a_{12}^*, b_{12}^*)$ such that $L_2'(v_{12}^*) < 0$. Moreover, if $\gamma_a^{12} = \gamma_b^{12}$, the condition $L_2'(v_{12}^*) < 0$ from (iv) in Proposition 6.2 implies that the equilibrium $(v_{12}^*, a_{12}^*, b_{12}^*)$ is locally asymptotically stable.

**7. Discussion and Numerical Simulations**

7.1. *Trajectories and Bistability*

The complete description of the stability of endemic steady states of the
nested immuno-epidemiological system (7)-(10) and (13)-(17) is theoretically
challenging, particularly when multiple equilibrium points coexist. To complete
this study, we perform numerical simulations. Assuming that all functions $\beta$
and $\mu$ are constant simplifies the model into an ODE system, which is easier
to simulate while still capturing the key phenomena. Although developing a
numerical scheme for our model is theoretically feasible, the ODE simplifica-
tion reduces computational time and avoids numerical instabilities, while still
capturing the essential qualitative behavior of the full PDE model. The simu-
lations were carried out in the symmetric case described in Theorem 5.6, where
$\rho(MC) < \mathcal{R}_0 < 1$, with $\rho(MC)$ given by (40). In this case, two bi-endemic
steady states coexist.

In Figure 10, we observe that even if $\mathcal{R}_0 < 1$ and the disease-free steady state is locally asymptotically stable, an attractive bi-endemic equilibrium exists with a large basin of attraction. This suggests that one of the bi-endemic equilibrium points is locally asymptotically stable. This phenomenon of bistability implies that, depending on the initial conditions of the epidemic, either both serotypes become endemic or they are eradicated. Bistability arises from a backward bifurcation, a phenomenon observed in many models with partial immunity, due to reinfection, such as in [3, 4, 11, 12]. In our model, a backward bifurcation occurs when the high infectivity of secondary infections can drive the system to a bi-endemic equilibrium, even if primary infections are not highly contagious. However, this effect becomes significant only if a sufficient number of individuals are infected by two different serotypes, which explains the observed bistability.

One of the consequences of this bistability is that a large and sudden influx
of individuals, i.e., when $\Lambda/\delta$ becomes large, can trigger a sudden transition
from a stable disease-free steady state to a stable bi-endemic equilibrium. A
bi-endemic equilibrium in the case of secondary dengue infection, for example,
is particularly worrying because it significantly increases the number of cases of
DHF. Another consequence of the hysteresis-type bifurcation diagram in Figure
8 is that, in the case of a backward bifurcation, eradicating just one of the
<figure><img src="image_11.png" /><figcaption>Figure 10: Trajectories of $S, I_1, I_2, S_1, S_2, I_{12}, I_{21}$, and $R$ in the symmetric case (Theorem 5.6) are shown. The initial values were set to $S(0) = 20, I_1(0) = I_2(0) = 0.001, S_1(0) = S_2(0) = 0, I_{12}(0) = I_{21}(0) = 0.1$ (left) or $0.09$ (right). Parameters were chosen such that the system is in the symmetric case of Theorem 5.6, with $\mathcal{R}_0 < 1$ and $MC < 1/2: \Lambda = 30, \delta = 1, \delta_1 = \delta_2 = 1.2, \delta_{12} = \delta_{21} = 1.5, \beta_1 = \beta_2 = 0.107, \beta_{12} = \beta_{21} = 1, \mu_1 = \mu_2 = 4$, and $\mu_{12} = \mu_{21} = 0.5$. Depending on the initial conditions, the system reaches either the bi-endemic equilibrium (left) or the disease-free equilibrium (right).</figcaption></figure>

serotypes becomes much more difficult when the population is already in a bi-endemic equilibrium favorable to DHF. The presence of two dengue serotypes in a population will therefore considerably increase the difficulty of controlling the epidemic.

When the bi-endemic equilibrium is reached and $\mathcal{R}_0^{(1)} > \mathcal{R}_0^{(2)}$ (serotype 1 is more transmissible than serotype 2), Figure 11 shows that $I_1^* > I_2^*$. However, we cannot determine a priori the order between $I_{12}^*$ and $I_{21}^*$. In the left part of Figure 11, we assume that $\beta_1 > \beta_2$ and $\beta_{12} < \beta_{21}$. Under these conditions, the equilibrium satisfies $I_{12}^* > I_{21}^*$. This result follows from the assumption that individuals are more likely to be firstly infected with serotype 1. Since the recovery rates are identical, the number of semi-recovered individuals $S_1$ is higher, leading to an increase in $I_{12}$. Conversely, if the recovery rate for primary serotype 1 infection is lower ($\mu_1 < \mu_2$), we obtain $I_{12}^* < I_{21}^*$ (see the purple lines in the right part of Figure 11). In this case, individuals firstly infected with serotype 2 recover more easily, leading to a higher number of $S_2$, which in turn facilitates secondary infections with serotype 1.

## 7.2. Relation between immunological and epidemiological parameters

Although our theoretical study provides an overview of immuno-epidemiological dynamics, it has certain limitations. In particular, we assume that the immuno-logical basic reproduction numbers are greater than 1 and that the viral load is constant (at non-trivial equilibrium). Even though this is a reasonable assumption for chronic diseases, dengue does not fit this framework. Furthermore, we assume that the force of infection $\beta$ and the recovery rate $\mu$ only depend on
<figure><img src="image_12.png" /><figcaption>Figure 11: Trajectories of $S, I_1, I_2, S_1, S_2, I_{12}, I_{21}$, and $R$ in non-symmetric cases are shown. The initial values were set to $S(0) = 20, I_1(0) = I_2(0) = 0.001, S_1(0) = S_2(0) = 0, I_{12}(0) = I_{21}(0) = 0.1$ (left) or $0.09$ (right). Parameters were chosen such that $\mathcal{R}_0 < 1$. For both figures, we have $\Lambda = 30, \delta = 1, \delta_1 = \delta_2 = 1.2, \delta_{12} = \delta_{21} = 1.5$, and $\mu_{12} = \mu_{21} = 0.5$. In the left figure, serotype 1 is more likely to be transmitted, with $\beta_1 = 0.12 > \beta_2 = 0.107$ and $\beta_{12} = 1 < \beta_{21} = 1.2$, while the recovery rates are identical: $\mu_1 = \mu_2 = 4$. In the right figure, we consider the same forces of infection $\beta_1 = \beta_2 = 0.107$ and $\beta_{12} = \beta_{21} = 1.2$, but different recovery rates: $\mu_1 = 3 < \mu_2 = 5$. In both cases, the bi-endemic equilibrium is reached. However, at equilibrium, we observe that $I_{12}^* > I_{21}^*$ on the left and $I_{12}^* < I_{21}^*$ on the right.</figcaption></figure>
the viral load $v$, which may be an oversimplification. To make our results more general, we supplement our study with numerical simulations.

For a given set of within-host parameters, we numerically compute the immunological trajectories $(v(\tau), a(\tau))$. Then, by calculating the integral in Expression (N1), we obtain the numerical values of $\mathcal{R}_0^{(m)}$ and $\mathcal{R}_0^{(nm)}$. By repeating this calculation for different sets of within-host parameters, we investigate the elasticity of the epidemiological invasion reproduction number $\mathcal{R}_2^{\text{inv}}$. The results of this method are shown in Figures 12 and 13.

In Figure 12, we observe that our simulations confirm the intuition derived from our theoretical study, despite the restrictive assumptions in Proposition 6.2. When $\mathcal{R}_0^{(2)}$ is sufficiently small, $\mathcal{R}_2^{\text{inv}}$ increases with respect to $\gamma_a^1$ and decreases with respect to $\gamma_v^1$. However, when $\mathcal{R}_0^{(2)}$ is large, the variations are reversed. This can be understood as follows: a reduction in the steady viral load of infection 1 leads to a higher proportion of individuals in $i_1$ recovering and transitioning to compartment $S_1$. Then, if serotype 2 is infectious enough, i.e., $\mathcal{R}_0^{(2)}$ is large, these individuals will be reinfected, entering compartment $i_{12}$, and subsequently infecting other susceptible individuals with serotype 2, which increases the chances of an invasion by serotype 2. On the other hand, if $\mathcal{R}_0^{(2)}$ is not high enough, there will still be a large number of individuals susceptible to serotype 2, but its infectivity will be too low to increase the chances of invasion. Meanwhile, serotype 1 will become more contagious, thus increasing its dominance. These results suggest that prioritizing treatment for one serotype while neglecting the other may inadvertently promote a transition to a bi-endemic equilibrium, leading to the persistence of both serotypes in the population and an increase in the number of DHF cases.

In Figure 13, we plot the invasion reproduction number $\mathcal{R}_2^{\text{inv}}$, as a function of the immunological parameters of the secondary infection. As expected, it is increasing as a function of $\gamma_b^{12}$, since the antibody b is neutralizing. However, we notice that the behavior of $\mathcal{R}_2^{\text{inv}}$ with respect to $\gamma_a^{12}$ changes depending on the value of $\xi_{12}$, i.e., on the strength of the ADE phenomenon. Indeed, $\xi_{12}$ was identified in the paper [10] as one of the main parameters for simulating the increase in pre-existing antibodies, thus promoting ADE. Concerning the effect of these immunological parameters on $\mathcal{R}_2^{\text{inv}}$, when antibody a is neutralizing, we find that $\mathcal{R}_2^{\text{inv}}$ simply increases with $\gamma_a^{12}$, just as for the antibody b. If $\xi_{12}$ becomes large, the antibody a becomes enhancing, and we see that the invasion reproduction number $\mathcal{R}_2^{\text{inv}}$ increases up to a certain maximum value, before decreasing towards a positive equilibrium. The same behavior has already been observed for the immunological reproduction number in [10]. Namely, there is a range of $\gamma_a^{12}$ where invasion becomes much more likely, increasing the number of secondary infections, and thus, of the number of DHF cases. This therefore indicates that the ADE phenomenon increases not only the probability of DHF in the host, but also the number of DHF cases in the population. In the paper [10], immunological data from patients with DHF during a secondary infection were analyzed to establish their relationship with the ADE phenomenon. This analysis could be supplemented and extended to epidemiological data, using the
<figure><img src="image_13.png" /><figcaption>Figure 12: Plot of $\mathcal{R}_2^{\text{inv}}$ as a function of $\gamma_a$ (left) or $\gamma_v$ (right), for different values of $\mathcal{R}_0^{(2)}$. The other functions and parameters were chosen as follows: $K = 1$, $\xi = 0.8$, $\mathcal{R}_0^{(12)} = 1$, $\Lambda = 1$, $\delta = 0.1$ and $\gamma_v = 1$ (left) or $\gamma_a = 1$ (right). The probability density $p_1$ of $v_0$ was chosen as uniform on $[0, K]$, $r(a) = (1 + 2a + \xi a^2)/(1 + 2a + a^2)$, $\eta(v) = 1 + v$, $\beta_1(v, a) = v/(1 + a)$, $\mu_1(v, a) = a/(1 + v)$. All parameter sets present in the plots are such that $\mathcal{R}_0^{(1)} > 1$ and $\mathcal{R}_0^{(2)} > \mathcal{R}_0^{(1)}$, so that the invasion reproduction number has meaning.</figcaption></figure>

immuno-epidemiological approach of this study. Other expected results include, for example, how to avoid promoting a bi-endemic equilibrium during treatment or vaccination, or how to prevent treatment that results in the clearance of antibodies from a primary infection and drives the invasion reproduction number close to its maximum.

## Acknowledgments

The authors thank the Inria international associate teams DIFaire and MoCoVec.

## References

[1] B. Adams and M. Boots, Modelling the relationship between antibody-dependent enhancement and immunological distance with application to dengue, *Journal of Theoretical Biology*, **242** (2006), 337–346, URL https://doi.org/10.1016/j.jtbi.2006.03.002.

[2] M. Adimy, P. F. A. Mancera, D. S. Rodrigues, F. L. P. Santos and C. P. Ferreira, Maternal passive immunity and dengue hemorrhagic fever in infants, *Bulletin of Mathematical Biology*, **82** (2020), 24, URL https://doi.org/10.1007/s11538-020-00699-x.
<figure><img src="image_14.png" /><figcaption>Figure 13: Plot of $\mathcal{R}_2^{\text{inv}}$ as a function of $\gamma_a^{12}$, for different values of $\xi_{12}$ (left) or $\gamma_b^{12}$ (right). We have chosen $\gamma_b^{12} = 1$ (on the left) and $\xi_{12} = 5$ (on the right). The other functions and parameters were chosen as follows: $K_{12} = 100$, $\Lambda = 50$, $\delta = 1$ and $\gamma_v^{12} = 2$. The probability densities of $v_0$ and $a_0$ were chosen as Dirac measures at $v_0 = 0.1$ and $a_0 = \eta_a(0)/\gamma_a^{12}$. The immunological functions are $r(a, b) = (1+2\xi_{12}a+0.95a^2)/(1+2a+a^2) \times (1+2b+0.85b^2)/(1+2b+b^2)$, $\eta_a(v) = \eta_b(v) = 1+v$, $\beta_{12}(v, a, b) = v/(1+a+b)$, $\mu_{12}(v, a, b) = (a+b)/(1+v)$. We have also chosen $\tilde{C}_1(\mathcal{R}_0^{(1)} - 1)/(\mathcal{R}_0^{(1)} - \mathcal{R}_0^{(2)}) = 1$.</figcaption></figure>

[3] R. Anguelov, S. M. Garba and S. Usaini, Backward bifurcation analysis of epidemiological model with partial immunity, *Computers and Mathematics with Applications*, **68** (2014), 931–940, URL https://doi.org/10.1016/j.camwa.2014.06.010.

[4] J. Arino, C. C. McCluskey and P. van den Driessche, Global Results for an Epidemic Model with Vaccination that Exhibits Backward Bifurcation, *SIAM Journal on Applied Mathematics*, **64** (2003), 260–276, URL https://doi.org/10.1137/S0036139902413829.

[5] L. Billings, I. B. Schwartz, L. B. Shaw, M. McCrary, D. S. Burke and D. A. Cummings, Instabilities in multiserotype disease models with antibody-dependent enhancement, *Journal of Theoretical Biology*, **246** (2007), 18–27, URL https://doi.org/10.1016/j.jtbi.2006.12.023.

[6] F. d. A. Camargo, M. Adimy, L. Esteva, C. Métayer and C. P. Ferreira, Modeling the relationship between antibody-dependent enhancement and disease severity in secondary dengue infection, *Bulletin of Mathematical Biology*, **83** (2021), 85, URL https://doi.org/10.1007/s11538-021-00919-y.

[7] H. E. Clapham, T. H. Quyen, D. T. H. Kien, I. Dorigatti, C. P. Simmons and N. M. Ferguson, Modelling virus and antibody dynamics during dengue virus infection suggests a role for antibody in virus clearance, *PLoS Com-
*putational Biology, 12 (2016), 1–15, URL https://doi.org/10.1371/journal.pcbi.1004951.*

[8] H. E. Clapham, V. Tricou, N. Van Vinh Chau, C. P. Simmons and N. M. Ferguson, Within-host viral dynamics of dengue serotype 1 infec- tion, *Journal of The Royal Society Interface*, **11** (2014), 20140094, URL https://doi.org/10.1098/rsif.2014.0094.

[9] O. Diekmann, J. A. P. Heesterbeek and J. A. J. Metz, On the definition and the computation of the basic reproduction ratio $r_0$ in models for infectious diseases in heterogeneous populations, *Journal of Mathematical Biology*, **28** (1990), 365–382, URL https://doi.org/10.1007/BF00178324.

[10] C. Dugourd-Camus, C. P. Ferreira and M. Adimy, Modelling the mech- anisms of antibody mixtures in viral infections: the cases of sequential homologous and heterologous dengue infections, *Journal of The Royal So- ciety Interface*, **21** (2024), 20240182, URL https://doi.org/10.1098/ rsif.2024.0182.

[11] M. Fang, Mathematical modelling, observation and identification of epidemiological models with reinfection, *Thesis at Laboratoire Jacques-Louis Lions of Sorbonne University*, (2024), URL https://theses.hal. science/tel-04848875v1

[12] Z. Feng, C. Castillo-Chavez and A. F. Capurro, A model for tuberculosis with exogenous reinfection, *Theoretical Population Biology*, **57** (2000), 235–247, URL https://doi.org/10.1006/tpbi.2000.1451.

[13] A. Gandolfi, A. Pugliese and C. Sinisgalli, Epidemic dynamics and host immune response: a nested approach, *Journal of Mathematical Biology*, **70**, 399–435, URL https://doi.org/10.1007/s00285-014-0769-8.

[14] M. A. Gilchrist and A. Sasaki, Modeling host–parasite coevolution: A nested approach based on mechanistic models, *Journal of Theoretical Biol- ogy*, **218** (2002), 289–308, URL https://doi.org/10.1006/jtbi.2002. 3076.

[15] M. G. M. Gomes, L. J. White and G. F. Medley, The reinfection threshold, *Journal of Theoretical Biology*, **236** (2005), 111–113, URL https://doi.org/10.1016/j.jtbi.2005.03.001.

[16] A. P. Goncalvez, R. E. Engle, M. St. Claire, R. H. Purcell and C.-J. Lai, Monoclonal antibody-mediated enhancement of dengue virus infection in vitro and in vivo and strategies for prevention, *Proceedings of the National Academy of Sciences*, **104** (2007), 9422–9427, URL https://doi.org/10.1073/pnas.0703498104.

[17] H. Gulbudak and C. J. Browne, Infection severity across scales in multi-strain immuno-epidemiological dengue model structured by host antibody
level, *Journal of Mathematical Biology*, **80** (2020), 1803–1843, URL https://doi.org/10.1007/s00285-020-01480-3.

[18] C. Gupta, N. Tuncer and M. Martcheva, A network immuno-epidemiological HIV model, *Bulletin of Mathematical Biology*, **83** (2021), 18, URL https://doi.org/10.1007/s11538-020-00855-3.

[19] S. B. Halstead, Pathogenesis of dengue: challenges to molecular biology, *Science*, **239** (1988), 476–481, URL https://doi.org/10.1126/science.3277268.

[20] S. B. Halstead, Dengue antibody-dependent enhancement: knowns and unknowns, *Microbiology Spectrum*, **2** (2014), URL https://doi.org/10.1128/microbiolspec.aid-0022-2014.

[21] M. Iannelli and F. Milner, *The Basic Approach to Age-Structured Population Dynamics: Models, Methods and Numerics*, Lecture Notes on Mathematical Modelling in the Life Sciences, Springer, 2017, URL https://doi.org/10.1007/978-94-024-1146-1.

[22] H. Inaba, *Age-Structured Population Dynamics in Demography and Epidemiology*, Springer Singapore, 2017, URL https://doi.org/10.1007/978-981-10-0188-8.

[23] M. Martcheva, *An Introduction to Mathematical Epidemiology*, vol. 61 of Texts in Applied Mathematics, Springer US, 2015, URL https://doi.org/10.1007/978-1-4899-7612-3.

[24] M. Martcheva, S. Lenhart, S. Eda, D. Klinkenberg, E. Momotani and J. Stabel, An immuno-epidemiological model for johne's disease in cattle, *Veterinary Research*, **46** (2015), 69, URL https://doi.org/10.1063/1.3659918.

[25] M. Martcheva, S. S. Pilyugin and R. D. Holt, Subthreshold and superthreshold coexistence of pathogen variants: The impact of host age-structure, *Mathematical Biosciences*, **207** (2007), 58–77, URL https://doi.org/10.1016/j.mbs.2006.09.010.

[26] N. Reich, S. Shrestha, A. King, P. Rohani, J. Lessler, S. Kalayanarooj, I. Yoon, R. Gibbons, D. Burke and D. Cummings, Interactions between serotypes of dengue highlight epidemiological impact of cross-immunity, *Journal of The Royal Society Interface*, **10** (2013), 20130414, URL https://doi.org/10.1098/rsif.2013.0414.

[27] V. Tricou, N. N. Minh, J. Farrar, H. T. Tran and C. P. Simmons, Kinetics of viremia and NS1 antigenemia are shaped by immune status and virus serotype in adults with dengue, *PLoS Neglected Tropical Diseases*, **5** (2011), 1–10, URL https://doi.org/10.1371/journal.pntd.0001309.