# Modelling the mechanisms of antibody mixtures in viral infections: the cases of sequential homologous and heterologous dengue infections

Charlotte Dugourd-Camus, Claudia Ferreira, Mostafa Adimy

► To cite this version:

Charlotte Dugourd-Camus, Claudia Ferreira, Mostafa Adimy. Modelling the mechanisms of antibody mixtures in viral infections: the cases of sequential homologous and heterologous dengue infections. Journal of the Royal Society Interface, 2024, 21 (219), pp.20240182. 10.1098/rsif.2024.0182. hal-04777129

HAL Id: hal-04777129

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

Submitted on 12 Nov 2024

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
# Modeling the mechanisms of antibody mixtures in viral infections: the cases of sequential homologous and heterologous dengue infections

Charlotte Dugourd-Camus<fnref n="1" />, Claudia P. Ferreira<fnref n="2" />, Mostafa Adimy<fnref n="1" />

## Abstract

Antibodies play an essential role in the immune response to viral infections, vaccination, or antibody therapy. Nevertheless, they can be either protective or harmful during the immune response. Moreover, competition or cooperation between mixed antibodies can enhance or reduce this protective or harmful effect. Using the laws of chemical reactions, we propose a new approach to modeling the antigen-antibody complex activity. The resulting expression covers not only purely competitive or purely independent binding but also synergistic binding which, depending on the antibodies, can promote either neutralization or enhancement of viral activity. We then integrate this expression of viral activity in a within-host model and investigate the existence of steady-states and their asymptotic stability. We complete our study with numerical simulations to illustrate different scenarios: firstly, where both antibodies are neutralizing, and secondly, where one antibody is neutralizing and the other enhancing. The results indicate that efficient viral neutralization is associated with purely independent antibody binding, whereas strong viral activity enhancement is expected in the case of purely competitive antibody binding. Finally, data collected during a secondary dengue infection were used to validate the model. The data set includes sequential measurements of virus and antibody titers during viremia in patients. Data fitting shows that the two antibodies are in strong competition, as the synergistic binding is low. This contributes to the high levels of virus titers and may explain the Antibody-Dependent Enhancement phenomenon. Besides, the mortality of infected cells is almost twice as high as that of susceptible cells, and the heterogeneity of viral kinetics in patients is associated with variability in antibody responses between individuals. Other applications of the model may be considered, such as the efficacy of vaccines and antibody-based therapies.

**Subject:** Biomathematics

**Keywords:** Humoral immunity response, Infectious diseases, Chemical reactions, Ordinary differential equations, Basic reproduction number, Local and global asymptotic stability.
34C60, 34D05, 37N25, 92C40, 92C45.

**Author for correspondence:** Charlotte Dugourd-Camus: charlotte.dugourd@inria.fr
Claudia P. Ferreira: claudia.pio@unesp.br
Mostafa Adimy: mostafa.adimy@inria.fr

## 1 Introduction

Active immunity is triggered by the entrance of a pathogen or by vaccination and involves the production of antibodies by the immune system. On the other hand, passive immunity involves the vertical transmission of antibodies from mother to offspring or the horizontal transmission of antibodies from humans or animals to susceptible individuals through antibody donation (antibody therapy). The interaction between the antibody and the antigen occurs through a specific chemical reaction, to form an antigen-antibody complex.

This immune complex can either neutralize or enhance the activity of the pathogen. The strength of the antigen-antibody complex depends firstly on the affinity of the antibody for the antigen, secondly on the number of antigen-antibody binding sites, and thirdly on the structural arrangement of the interacting parts, [1, 2]. In the case of secondary infection by a homologous or heterologous virus, pre-existing antibodies and newly produced antibodies may coexist, [3]. In most cases, an active infection leads to the production of neutraliz-

<footnote n="1">¹Inria, ICJ UMR5208, CNRS, Ecole Centrale de Lyon, INSA Lyon, Université Claude Bernard Lyon 1, Université Jean Monnet, 69603 Villeurbanne, France</footnote>

<footnote n="2">²São Paulo State University (UNESP), Institute of Biosciences, 18618-689 Botucatu, SP, Brazil</footnote>
ing antibodies. However, cross-reaction of pre-existing
antibodies may or may not help neutralize this ac-
tive infection; in fact, it has been observed that cross-
reacting antibodies can enhance the infection, [4–8]. In
the case of influenza, dengue, or COVID-19, for exam-
ple, it has been observed that two successive homolo-
gous infections can increase the neutralizing effect of
antibodies, [1, 9]. On the other hand, secondary in-
fection with a heterologous dengue serotype may in-
crease the risk of developing a severe form of the dis-
ease, due to Antibody-Dependent Enhancement (ADE)
phenomenon, [4–8, 10–12]. Furthermore, in the pres-
ence of a mixture of two antibodies and a pathogen,
two other factors come into play: (i) the competition
between antibodies to bind to the antigen, and (ii) the
synergistic interactions between antibodies, [1,13]. This
leads to three distinct types of antigen-antibody bind-
ing: (a) purely independent binding, where the binding
of one antibody does not affect the binding of other an-
tibodies to the same receptor; (b) purely competitive
binding, where the two antibodies cannot bind simulta-
neously to the same receptor; and (c) an intermediate
situation, where the two antibodies can bind simulta-
neously to many sites of the same receptor, with syner-
gistic interactions, i.e., the binding of one has an effect
on the binding of the other, either by making it weaker
or stronger, [1, 13].

Several mathematical within-host infectious disease
models have been developed to study the dynamics
of viral infections and their interaction with the im-
mune system, [14–26]. They involve, according to each
paper, healthy and infected target cells, intracellular
pathogen replication, T-cells, B-cells, cytokine produc-
tion, and antibodies. Most of them resemble epidemio-
logical models and use the mass action law or a satura-
tion function to model the force of infection. In general,
only one of the immune responses - humoral or cellu-
lar - is modeled, and very few models address interac-
tion among components of the immune response. As an
example, [15, 18, 20] took into account only one compo-
nent of the immune system: T-cells, the only ones capa-
ble of killing infected cells. In particular, [15] used the
Beddington-DeAngelis incidence rate to model interac-
tion between susceptible cells and free viruses. They
argued that they could thus reduce the time needed for
the immune response to clean the dengue virus, com-
pared with using the mass action law of the paper [20].
In both articles, depending on parameter values, more
than one endemic steady-state was found, and thresh-
olds for the stability of each of them have been attained.
In [18], the authors used data coming from hospitalized
primary and secondary dengue patients - virus viremia
- to estimate model parameters that confirmed the role
of the immune response in shaping variation between
individuals consistent with the hypothesis of ADE. In
[14], the authors developed deterministic and stochas-
tic within-host models to explore different dengue infec-
tion scenarios, taking into account individual immuno-

logical variability. Their models are calibrated using
empirical data on viral load and antibody concentra-
tions (IgM and IgG), incorporating confidence intervals
derived from stochastic realizations. In [19], the au-
thors considered heterologous antibodies that interact
through a Heaviside step function that triggers neutral-
ization or enhancement of the infection depending on
the amount of the antibody from the first infection. The
model also considered a delay in the immune response
which switches the stability of the system through a
Hopf bifurcation. The authors argued that they could
qualitatively replicate the humoral immune responses
observed in primary and secondary infections. The au-
thors of [16] have shown that the risk of developing a
severe form of dengue may be related to increased cy-
tokine production due to the interaction of the immune
system (T-cells) with the dengue virus. The model was
parameterized to reproduce qualitatively the data set
described in [18]. The virological indicators used were
the level of peak viremia, the time to peak viremia, and
the viral clearance rate. Although the model did not
consider the humoral immune response explicitly, the
authors argue that the reparametrization of the infec-
tivity rate can be biologically interpreted in the context
of the ADE phenomenon.

Other original approaches have been developed in
[17,27,28]. All, explicitly considered interactions among
susceptible and infected target cells, dengue virus, and
dengue antibodies through bilinear and trilinear terms.
In [17], two thresholds were obtained: the mean num-
ber of virions produced by one invading virus in the
very early stage of secondary infection, and the rate at
which one infected macrophage dies. The existence and
stability of the endemic equilibrium depend on both
parameters, and one of them acts as a weakening fac-
tor for ADE. In [23], the authors proposed a within-
host mathematical model to evaluate the role of mem-
ory B and T cells in heterologous secondary dengue
infection. They showed that memory T cells play an
essential role in eliminating the possibility of ADE oc-
currence. In [24], the authors developed a within-host
model to study the impact of ADE phenomena on dis-
ease severity of Zika virus and dengue virus sequential
or co-infection. In [26], the authors developed a within-
host model for primary and secondary dengue infec-
tions, taking into account the IgM and IgG antibody
response. The model suggests that a faster rate of clear-
ance of antibody-virus complexes may lead to a higher
peak viral load and may explain the ADE phenomena
in heterogeneous dengue infections. In [27], the authors
explored the occurrence of Dengue Hemorrhagic Fever
(DHF) in infants born from dengue-immune mothers,
during their first dengue infection. The neutralizing
and enhancing activities of maternal antibodies against
the virus are represented by a function derived from
experimental data. They were able to fit the model
to data on the amount of maternal antibodies and the
age of the infant at which DHF was reported. The
authors were thus able to reproduce the time delay observed between the end of the protective level of maternal antibodies and the onset of hemorrhagic fever, in agreement with data from the literature. The article [28] focused on the ADE hypothesis and developed a mathematical model of secondary dengue infection by a different viral serotype. It considered indirect competition between neutralizing and enhancing antibodies. Here, the functions of enhancement and neutralization - which depend on the amount and type of antibody - are derived from basic concepts of chemical reactions and used to model virus-antibody complexes binding formed by distinct populations of antibodies, classified as cross-reactive or type-specific ones. The authors of this article concluded that virus-antibody immune complexes may promote viral clearance or enhancement of infection depending on the amount of cross-reacting antibodies and the rapid activation of the neutralizing antibodies.

The authors of [22, 25] considered the modeling of within-host dynamics of HIV infection and therapy. Both are review articles and can give a broad understanding of how the immune system interacts with viruses and the state of the art of HIV modeling which includes virus evolution. The authors claimed that HIV viruses have several different epitopes that can be recognized by immune response. Furthermore, different populations of HIV viruses and antibodies when mixed, may have cross-reactive responses. In [21], the proposed model includes both innate and adaptive immune responses, during an influenza virus infection. Different from the other models, a class of uninfected cells that are refractory to infections was included. Modeling predictions were compared with both interferon and viral kinetic data. The first post-peak viral decline is explained by the lysis of infected cells during the innate immune response, the subsequent viral plateau/second peak is generated by the loss of the IFN-induced antiviral effect and the increased availability of target cells.

None of the cited works either considered the direct competition of the antibodies for the virus epitope or the synergistic interaction between antibodies. Therefore, the work presented here aims to introduce a new formalism, inspired by [1, 13], to efficiently describe the formation of antigen-antibody complexes, and to use it in a mathematical model to describe the interaction between a virus, target cells and two antibodies competing to bind to virus receptors, taking into account their synergistic interaction, and the effectiveness of the complex formed in neutralizing or enhancing the virus. To achieve this, we first consider a monoclonal antibody and define the viral activity function of the antigen-

antibody complex. We then consider a mixture of two antibodies and generalize the viral activity function obtained previously, which now depends on the amount of the two antibodies and the interaction between them. Finally, we use this viral activity function in a within-host infectious disease model. We derive some basic properties of the within-host model obtained, in particular the existence of steady-states (disease-free and endemic), and investigate the local and global asymptotic stability of these steady-states. We complete our study with numerical simulations to highlight various scenarios. Although we have parameterized the model to study homologous and heterologous secondary infection with the dengue virus, it can be easily adapted to other viruses. It can also be extended to vaccination or antibody therapy, [29, 30]. To validate the proposed formalism, we fitted the within-host infection model to the data from [4]. The data include virus and antibody titer measurements recorded sequentially during a secondary infection from dengue patients. To understand the mechanisms underlying the heterogeneity of viral kinetics in patients, we focus on estimating parameters linked to the force of infection. The results show that the IgG and IgM antibodies generated during the first and second infections, respectively, have low synergistic binding. This means that they are in strong competition to bind to the virus which contributes to the enhancement of the infection. In addition, the strength of the interaction between antibodies varies from one individual to another, which explains the differences in the temporal evolution of virus and antibody populations and, therefore, the spectrum of dengue infection that includes the phenomenon of antibody-dependent enhancement (ADE).

# 2 Modeling the antigen-antibody complex

## 2.1 Monoclonal antibody binding to a receptor

Consider a monoclonal antibody A that binds to an antigen R (viral receptor) through a paratope-epitope bond and, depending on its affinity, inhibits or enhances viral activity. An antibody can recognize a number n of epitopes of the same antigen, [12]. Antigen-antibody binding occurs through a chemical reaction, [1], in which n paratopes N of A bind to n epitopes of R. According to the standard law of mass action, the process leading to this bond is given by the following sequential chemical reaction:

$$R + nN \underset{k^-}{\overset{nk^+}{\rightleftharpoons}} R_{1N} + (n-1)N \underset{2k^-}{\overset{(n-1)k^+}{\rightleftharpoons}} \dots \underset{ik^-}{\overset{(n-i+1)k^+}{\rightleftharpoons}} R_{iN} + (n-i)N \underset{(i+1)k^-}{\overset{(n-i)k^+}{\rightleftharpoons}} \dots \underset{nk^-}{\overset{k^+}{\rightleftharpoons}} R_{nN}, \quad (1)$$
where $R_{iN}$ stands for the complex formed by $i$ epitopes of $R$ bound by $i$ paratopes $N$ of $A$, $R = R_{0N}$ means free antigen, $k^+$ and $k^-$ denotes the forward and backward reaction rates, [31]. When the chemical reaction reaches equilibrium, we obtain the following equation

$$(n - i)k^+[R_{iN}][N]^{n-i} = (i + 1)k^-[R_{(i+1)N}][N]^{n-i-1},$$

valid for $i \in \{0, 1, \dots, n-1\}$, where $[\cdot]$ is the concentration of the chemical species. After some algebraic manipulation, we get

$$[R_{iN}] = \frac{n - i + 1}{i} \frac{[N]}{K} [R_{(i-1)N}], \quad \text{for } i \in \{1, \dots, n\},$$

where $K = k^-/k^+$ is the dissociation constant of the reaction.

Solving this recurrence equation yields to

$$[R_{iN}] = C_n^i ([N]/K)^i [R],$$

with $C_n^i = \frac{n!}{i!(n-i)!}$, the binomial coefficient. Given the total antigen concentration

$$[R_{\text{total}}] = \sum_{i=0}^{n} [R_{iN}] = \left(1 + \frac{[N]}{K}\right)^n [R],$$

and knowing that the antibody concentration $A$ is proportional to the concentration of paratope $N$, $A = \theta[N]$, the fraction of occupied epitopes of the antigen is given by

$$\frac{[R_{iN}]}{[R_{\text{total}}]} = \frac{C_n^i (A/K_D)^i}{(1 + A/K_D)^n}, \quad (2)$$

with $K_D = \theta K$, where $\theta$ is a constant. The antigen-antibody complex modifies viral activity. Let's define the relative activity of a virus belonging to this complex as $\xi_i \ge 0$, where $i$ represents the number of occupied epitopes of the antigen. The value $0 < \xi_i < 1$ means that the complex antigen-antibody partially inhibits virus activity, while $\xi_i > 1$ means that it increases virus activity. The case $\xi_i = 1$ means that the relative virus activity remains unchanged. In particular, $\xi_0 = 1$, because no antigen-antibody complex is formed. As an example, in [32], it was estimated that below 130 IgG or 30 IgM bound per virions, infectivity is totally preserved. Furthermore, binding to defective entry-mediating proteins on virions would not be directly relevant for neutralization, [33]. Therefore, the activity of the virus is obtained by adding the products of the fractions (2) by their associated relative activities $\xi_i$,

$$\begin{aligned} \text{Activity} := G(A) &= \sum_{i=0}^{n} \xi_i \times \frac{C_n^i (A/K_D)^i}{(1 + A/K_D)^n}, \\ &= \frac{1 + \sum_{i=1}^{n} \xi_i C_n^i (A/K_D)^i}{1 + \sum_{i=1}^{n} C_n^i (A/K_D)^i}. \end{aligned} \quad (3)$$

<figure><img src="image_2.png" /><figcaption>Figure 1: Viral activity as a function of monoclonal antibody $A$, $G(A)$ (Equation (3)). Each curve is associated with a specific set of relative viral activities $(\xi_1, \xi_2, \xi_3)$, $K_D = 0.7$ and $n = 3$. Red curves highlight the antigen-antibody complex leading to enhancement, while blue curves are those leading to neutralization. The black curve means that virus activity is neither neutralized nor enhanced. Depending on the amount of antibody, virus activity can change from enhancement to neutralization (red-blue dotted line).</figcaption></figure>

The binding affinity between the antigen's epitope and the antibody's paratope at a single binding site can be interpreted in different ways. For example, the higher the affinity, the lower the values of $K_D$ and $\xi_i$. This explains the property of neutralizing antibodies, which have a high affinity for the virus and, after binding to the receptor, neutralize it, $G(A) < 1$ (blue lines in Figure 1). Non-neutralizing cross-reactive antibodies, on the other hand, only partially recognize the receptor (low affinity), [12]. The constant $K_D$ is therefore high, and $\xi_i$ are close to 1. In particular, if all relative activities $\xi_i$ are equal to 1, the virus activity remains unaffected, $G(A) = 1$ (black line in Figure 1). For some viral infections, such as secondary dengue infection, non-neutralizing cross-reactive antibodies may enhance the infection. This means an increase in viral activity, $G(A) > 1$ (red lines in Figure 1). As the interaction antigen-antibody is a multi-hit phenomenon, [34], enhancement can occur for an intermediate amount of antibody, while neutralization occurs in the presence of a sufficiently large amount of antibody, [35, 36]. This suggests that if the number of occupied epitopes is high enough, epitope-paratope binding will neutralize antigen activity; and if the number of occupied epitopes is intermediate, this binding will tend to enhance antigen activity (red-blue dotted line in Figure 1).
## 2.2 Binding of a mixture of two antibodies to a receptor

In [1], the authors have developed a statistical mechanical model that predicts the collective efficacy of a mixture of antibodies whose constituents are assumed to bind to a single site on a receptor. We generalize their method to the case where the mixture can bind to multiple sites on a receptor, [12]. Let's consider now a

mixture of two different antibodies $A_1, A_2$ that interact with a virus receptor $R$. This mixture can result in an antigen-antibody complex with purely independent bindings, purely competitive bindings, or synergistic bindings. This binding classification is based on the interactions (competition/cooperation) between the two antibodies. In this case, the activity of the virus is defined by the function $G(A_1, A_2)$, given by,

$$G(A_1, A_2) = \frac{1 + \sum_{i=1}^{n_1} \xi_{1i} C_{n_1}^i \left( \frac{A_1}{K_{D_1}} \right)^i + \sum_{j=1}^{n_2} \xi_{2j} C_{n_2}^j \left( \frac{A_2}{K_{D_2}} \right)^j + \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} f_{ij} \tilde{\xi}_{1i} \tilde{\xi}_{2j} C_{n_1}^i C_{n_2}^j \left( \frac{A_1}{\tilde{K}_{D_1}} \right)^i \left( \frac{A_2}{\tilde{K}_{D_2}} \right)^j}{1 + \sum_{i=1}^{n_1} C_{n_1}^i \left( \frac{A_1}{K_{D_1}} \right)^i + \sum_{j=1}^{n_2} C_{n_2}^j \left( \frac{A_2}{K_{D_2}} \right)^j + \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} f_{ij} C_{n_1}^i C_{n_2}^j \left( \frac{A_1}{\tilde{K}_{D_1}} \right)^i \left( \frac{A_2}{\tilde{K}_{D_2}} \right)^j} \quad (4)$$

As in (3), $n_1$ (resp. $n_2$) is the number of epitopes of the antigen that can be recognized by the antibody $A_1$ (resp. $A_2$), $\xi_{1i}$ (resp. $\xi_{2j}$) is the relative activity of the virus when bound to $i$ paratopes of $A_1$ (resp. $j$ paratopes of $A_2$), and $K_{D_1}$ (resp. $K_{D_2}$) is the dissociation constant associated to the antibody $A_1$ (resp. $A_2$). The new parameters $\tilde{\xi}_{1i}, \tilde{\xi}_{2j}$ are the relative virus activities modified by a synergistic binding (the relative virus activity may decrease in the presence of the other antibody because of competition, $\tilde{\xi}_{kl} < \xi_{kl}$, or it can increase because of cooperation, $\tilde{\xi}_{kl} > \xi_{kl}$). Similarly, a synergistic interaction between two antibodies can modify their binding to the receptor and thus their dissociation constants $K_{D_k}$, $k=1,2$. We then introduce

modified dissociation constants $\tilde{K}_{D_k}$, with $\tilde{K}_{D_k} < K_{D_k}$ when antigen-antibody binding becomes stronger and $\tilde{K}_{D_k} > K_{D_k}$ when antigen-antibody binding becomes weaker. The coefficient $0 \le f_{ij} \le 1$ corresponds to the fraction of simultaneous binding of both antibodies. Thus, three scenarios are possible: (i) purely competitive binding, with $f_{ij} = 0$; (ii) synergistic binding, where $0 < f_{ij} \le 1$; and (iii) purely independent binding, where $f_{ij} = 1$, $\tilde{\xi}_{ki} = \xi_{ki}$ and $\tilde{K}_{D_k} = K_{D_k}$, for $k=1,2$, $i=1, \dots, n_1$ and $j=1, \dots, n_2$. In the case of [1], $n_1 = n_2 = 1$. The purely competitive binding and the purely independent binding lead respectively to the following expressions of the virus activity

$$G(A_1, A_2) = \frac{1 + \sum_{i=1}^{n_1} \xi_{1i} C_{n_1}^i \left( \frac{A_1}{K_{D_1}} \right)^i + \sum_{j=1}^{n_2} \xi_{2j} C_{n_2}^j \left( \frac{A_2}{K_{D_2}} \right)^j}{1 + \sum_{i=1}^{n_1} C_{n_1}^i \left( \frac{A_1}{K_{D_1}} \right)^i + \sum_{j=1}^{n_2} C_{n_2}^j \left( \frac{A_2}{K_{D_2}} \right)^j}, \quad (5)$$

and

$$G(A_1, A_2) = \left( \frac{1 + \sum_{i=1}^{n_1} \xi_{1i} C_{n_1}^i \left( \frac{A_1}{K_{D_1}} \right)^i}{\left( 1 + \frac{A_1}{K_{D_1}} \right)^{n_1}} \right) \times \left( \frac{1 + \sum_{j=1}^{n_2} \xi_{2j} C_{n_2}^j \left( \frac{A_2}{K_{D_2}} \right)^j}{\left( 1 + \frac{A_2}{K_{D_2}} \right)^{n_2}} \right). \quad (6)$$

To study the synergistic binding effect on viral activity, we introduce the following notations $\tilde{\nu}_{ij} = \tilde{\xi}_{1i}\tilde{\xi}_{2j}/(\xi_{1i}\xi_{2j})$ and $\tilde{\kappa}_{ij} = K_{D_1}^i K_{D_2}^j/(\tilde{K}_{D_1}^i \tilde{K}_{D_2}^j)$, and explore the values of these relative quantities. Therefore, the expression (4) becomes

$$G(A_1, A_2) = \frac{1 + \sum_{i=1}^{n_1} \xi_{1i} C_{n_1}^i \left( \frac{A_1}{K_{D_1}} \right)^i + \sum_{j=1}^{n_2} \xi_{2j} C_{n_2}^j \left( \frac{A_2}{K_{D_2}} \right)^j + \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} \tilde{\nu}_{ij} \tilde{\kappa}_{ij} f_{ij} \xi_{1i} \xi_{2j} C_{n_1}^i C_{n_2}^j \left( \frac{A_1}{K_{D_1}} \right)^i \left( \frac{A_2}{K_{D_2}} \right)^j}{1 + \sum_{i=1}^{n_1} C_{n_1}^i \left( \frac{A_1}{K_{D_1}} \right)^i + \sum_{j=1}^{n_2} C_{n_2}^j \left( \frac{A_2}{K_{D_2}} \right)^j + \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} \tilde{\kappa}_{ij} f_{ij} C_{n_1}^i C_{n_2}^j \left( \frac{A_1}{K_{D_1}} \right)^i \left( \frac{A_2}{K_{D_2}} \right)^j} \quad (7)$$
The new parameters $\tilde{\nu}_{ij}$ and $\tilde{\kappa}_{ij}$ can be interpreted as follows. If $\tilde{\nu}_{ij} < 1$, the relative viral activity due to synergistic binding decreases. On the other hand, if $\tilde{\nu}_{ij} > 1$, relative viral activity increases. Furthermore, if $\tilde{\kappa}_{ij} < 1$, the binding between the antibodies and the virus is enhanced by synergistic binding, and if $\tilde{\kappa}_{ij} > 1$ this binding becomes weaker. In theory, many combinations are possible, but we will explore the most relevant ones: when both antibodies neutralize the virus activity, and when one antibody neutralizes and the other enhances the virus activity.

Figure 2 shows the behavior of viral activity as a function of the two antibodies $A_1$ and $A_2$. The case of purely independent binding is explored. In Figure 2(a), the two antibodies act to neutralize the virus. Whatever the direction, horizontal or vertical from the left or bottom border, we always move towards lower viral activity, which means more neutralization. This is coherent with the assertion of [1,9] that two antibodies are more effective than one in neutralizing a virus. In the right panel of Figure 2, i.e. Figure 2(b), it is shown the variation in viral activity in the case of enhancing antibody $A_1$ and neutralizing antibody $A_2$. For each fixed amount of neutralizing antibody, $A_2$, viral activity first increases very rapidly as the amount of enhancing antibody $A_1$ increases, reaching a peak in the red zone, then slowly decreases (moving from red to blue) until reaches a threshold in the blue (same behavior of the red/blue line of Figure 1). On the other hand, if we fix the amount of enhancing antibody, $A_1$, and increase the amount of the neutralizing one, $A_2$, we obtain different scenarios: (i) If $A_1$ is fixed at very low or high levels, viral activity will decrease rapidly at first, then reach a certain value (remain in the blue zone once entered, meaning that viral activity is neutralized); (ii) If $A_1$ is fixed at an intermediate level, viral activity will start at a high level and change very quickly from red to blue or from dark red to light red, then return to an intermediate value and remain at this level. The main information to be deduced from Figure 2(b) is that infection is more severe when the amount of enhancing antibodies is in the intermediate range. This is compatible with the assertion of [28,35,36] that the risk of severe dengue in a secondary heterologous infection is higher when there is an intermediate amount of pre-existing antibodies from the first infection. The common values of the parameters in Figure 2 are $K_{D_1} = K_{D_2} = 0.7$ mol ml$^{-1}$ and $\beta = 8 \times 10^{-9}$ ml RNA copies$^{-1}$. In Figure 2(a), we take $(\xi_{11}, \xi_{12}, \xi_{13}) = (\xi_{21}, \xi_{22}, \xi_{23}) = (1, 0.95, 0.65)$, and in Figure 2(b), we take $(\xi_{11}, \xi_{12}, \xi_{13}) = (1, 4, 0.85)$ and $(\xi_{21}, \xi_{22}, \xi_{23}) = (1, 0.4, 0.85)$.

The immune response to a viral infection can be extremely complex, orchestrated by a complex interaction between elements of the innate and adaptive immune response that depends on a number of factors such as host genetics, heterologous immunity, virus structure and receptor recognition, and avoidance of the immune response by the virus. Mechanisms preventing

virion attachment to target cell receptors include binding at or near the viral receptor binding site and prevention of attachment by steric obstruction, disassembly or conformational modification of viral surface entry proteins, and virion aggregation. In addition, the involvement of antibodies in virus dynamics also includes post-attachment neutralization, [37]. Nevertheless, in this article we try to describe the interactions between antibodies in a simple way. For example, we assume that viral surface molecules are static, but it is increasingly accepted that these molecules undergo conformational relaxation over time that exposes new epitopes or increases the exposure of existing epitopes. Cellular factors limit or enhance neutralization and indicate whether combinations produce additive, synergistic or antagonistic net effects. The list of antibody-virus interactions includes incomplete neutralization, synergistic neutralization, additive neutralization and enhancement; it does not exclude the existence of other interactions, [33].

## 3 Within-host (re)infection dynamics in the presence of a mixture of two antibodies

As an example, let us consider a viral (re)infection at the cellular and immune level with a mixture of two antibodies. We note by $A_1$ the concentration of pre-existing antibody due to the first infection and by $A_2$ the concentration of the new antibody generated by the secondary infection with a homologous or heterologous virus. The quantities $X$ and $Y$ are the susceptible and infected target cells respectively, and $V$ is the free virus. Target cells may be epithelial cells in the case of influenza, [38] or macrophages and dendritic cells in the case of dengue, [39]. The model is described by the ordinary differential system given by

$$ 
\begin{cases}
\frac{dA_1}{dt} = p_1(Y) - \gamma_{A_1} A_1, \\
\frac{dA_2}{dt} = p_2(Y) - \gamma_{A_2} A_2, \\
\frac{dX}{dt} = \Omega - \gamma_X X - \beta G(A_1, A_2) V X, \\
\frac{dY}{dt} = \beta G(A_1, A_2) V X - \gamma_Y Y, \\
\frac{dV}{dt} = \Phi Y - (\delta + \gamma_V) V.
\end{cases}
\quad (8) $$

The antibodies $A_i$ can be produced by memory B-cells or by plasmablasts upon stimulation by the virus or by infected cells, [40]. In the case of a secondary viral infection, the first antibody $A_1$ is generated mainly by memory B-cells, while the second antibody $A_2$ is generated mainly by plasmablasts. It is reasonable to assume that the production of both antibodies depends on the concentration of infected cells. The production of antibody $A_i$ is then considered with a rate equal to $p_i(Y)$,
<figure><img src="image_3.png" /><figcaption>Figure 2: Viral activity $G(A_1, A_2)$, in the case of Neutralizing/Neutralizing (left panel) and Enhancing/Neutralizing (right panel) antibodies. The case of purely independent binding is explored (Equation (6)). The palette of colors, from blue (neutralization) to red (enhancement), reflects the variation in viral activity on antigen-antibody complexes.</figcaption></figure>

$p_i$ being a nonnegative, continuously differentiable, and nondecreasing function on $[0, +\infty)$. For all numerical simulations, we use the function $p_i(Y) = \Lambda_i + \alpha_i Y$. The natural mortality rates are denoted by $\gamma_{A_i}$, $\gamma_X$, $\gamma_Y$ and $\gamma_V$. The parameter $\Omega$ is the rate of production of susceptible target cells, which are generally produced in the bone marrow, [41]. $\Phi$ is the rate of virus production by infected cells, [42, 43]. The parameter $\delta$ represents the rate of free virus loss by other means than natural mortality, such as neutralization or entry into target cells. When the two antibodies $A_1$ and $A_2$ interact with the virus, its activity may decrease or increase depending on the nature of these antibodies. This activity has been defined in Subsection 2.2 by the function $G(A_1, A_2)$ given by (7). Then, the force of infection is $\beta G(A_1, A_2) V$, with $\beta > 0$. This system can model a homologous or heterologous secondary viral infection, such as dengue fever or influenza. It can also be adapted to antiretroviral therapy in the case of HIV, provided we change the functions $p_i$ and adapt the parameter values.

Obviously, a simple mathematical model cannot fully reflect the real behavior of all the examples cited above. However, we are interested in an important part of the immune response, common to all these diseases, which concerns the mixture of competing or cooperating antibodies during infection.

Throughout this paper, we will need to make certain assumptions.

$$ (H1) \quad \gamma := \gamma_X = \gamma_Y. $$

$$ (H2) \quad \gamma_A := \gamma_{A_1} = \gamma_{A_2}. $$

$$ (H3) \quad G(A_1, A_2) \le G\left(\frac{p_1(0)}{\gamma_{A_1}}, \frac{p_2(0)}{\gamma_{A_2}}\right), \quad \text{for all } (A_1, A_2) \in \mathbb{R}_+^2. $$

D, with

$$ \mathbb{D} := \left\{ (A_1, A_2) \in \mathbb{R}_+^2 : A_1 \ge \frac{p_1(0)}{\gamma_{A_1}}, A_2 \ge \frac{p_2(0)}{\gamma_{A_2}} \right\}. $$

$$ (H4) \quad p_i(Y) - p_i(0) \le p'_i(0)Y, \quad \text{for } i = 1, 2, \text{ and all } Y \ge 0. $$

$$ (H5) \quad p'_1(0) > 0 \text{ or } p'_2(0) > 0. $$

The assumptions (H1) and (H2), mean that the disease does not affect cell mortality and that both antibodies have the same mortality rate. They are used, for the sake of simplicity, in the proof of the global asymptotic stability of $P_0$ and/or in the proof of the local asymptotic stability of $P^*$, the disease-free equilibrium and the endemic one, respectively. Under (H3), quantities of antibodies $A_1$ and $A_2$ higher than those of the disease-free equilibrium reduce viral activity. This gives an advantage to neutralization. Hypothesis (H4) means that the dynamic of antibody production follows, at best, a linear growth. If (H5) is satisfied, then $p_1$ or $p_2$ is strictly increasing in a neighborhood of 0, which implies a production strictly positive of $A_1$ or $A_2$ around 0. The hypotheses (H3), (H4) and (H5) are assumed satisfied in Theorem 4.3 and Theorem 4.5.

# 4 Mathematical analysis of the within-host (re)infection model

In this section, we analyze the ordinary differential system (8). In particular, we establish the existence, uniqueness, and positivity of solutions, determine the existence of steady-states, calculate the basic reproduction number, and investigate the local and global
asymptotic stability of the disease-free steady-state and the local asymptotic stability of the endemic equilibrium. For convenience, all the proofs of the theoretical results (Theorems, Propositions and Lemma) are given in the appendix, Section A.

## 4.1 Basic properties, steady-states, and basic reproduction number

**Proposition 1.** The solution of the initial value problem (8), associated with a nonnegative initial condition, is unique, nonnegative, and bounded on $[0, +\infty)$.

*Proof.* The proof is given in the appendix (Subsection AA.2). □

Let $P = (A_1^*, A_2^*, X^*, Y^*, V^*)$ be an equilibrium:

$$ 
\left\{ 
\begin{array}{lll}
p_1(Y^*) - \gamma_{A_1} A_1^* & = 0, & \\
p_2(Y^*) - \gamma_{A_2} A_2^* & = 0, & \\
\Omega - \gamma_X X^* - \beta G(A_1^*, A_2^*) V^* X^* & = 0, & (9) \\
\beta G(A_1^*, A_2^*) V^* X^* - \gamma_Y Y^* & = 0, & \\
\Phi Y^* - (\gamma_V + \delta) V^* & = 0. & 
\end{array}
\right. 
$$

By solving (9), we obtain: the disease-free equilibrium $P_0$ and endemic equilibrium points $P^*$. The disease-free equilibrium is given by

$$ P_0 = \left( \frac{p_1(0)}{\gamma_{A_1}}, \frac{p_2(0)}{\gamma_{A_2}}, \frac{\Omega}{\gamma_X}, 0, 0 \right). \quad (10) $$

The next-generation matrix, [44], is used to define the bifurcation parameter $\mathcal{R}_0$. For this, we consider the two dimensions infected subsystem - $(Y, V)^T$ in (8) - that describe the production of new infections and changes in the state of the infected individuals

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{dY}{dt} = \beta G(A_1, A_2) VX - \gamma_Y Y, \\
\displaystyle \frac{dV}{dt} = \Phi Y - (\delta + \gamma_V)V.
\end{array}
\right. 
$$

The rates of appearance of the newly infected individuals are given by

$$ W = \begin{pmatrix} 0 & \beta G \left( \frac{p_1(0)}{\gamma_{A_1}}, \frac{p_2(0)}{\gamma_{A_2}}, \frac{\Omega}{\gamma_X} \right) \\ 0 & 0 \end{pmatrix}, $$

and the rates of transfer of individuals within the infected compartment, by any means other than the appearance of newly infected individuals from the uninfected compartment, [44], are given by

$$ T = \begin{pmatrix} -\gamma_Y & 0 \\ \Phi & -(\gamma_V + \delta) \end{pmatrix}. $$

These two matrices are obtained by decomposing the Jacobian matrix into two parts, the transmission $W$ and the transition $T$, which are evaluated at the disease-free

equilibrium $P_0 := (p_1(0)/\gamma_{A_1}, p_2(0)/\gamma_{A_2}, \Omega/\gamma_X, 0, 0)$. Finally, the bifurcation parameter $\mathcal{R}_0$ is the spectral radius of the matrix product $-WT^{-1}$, [44],

$$ \mathcal{R}_0 = \rho(-WT^{-1}) = \frac{\Phi\beta\Omega}{\gamma_Y\gamma_X(\gamma_V + \delta)} G \left( \frac{p_1(0)}{\gamma_{A_1}}, \frac{p_2(0)}{\gamma_{A_2}} \right). $$

$\mathcal{R}_0$ measures the expected number of infected cells (or viral particles) generated by an infected cell (or viral particle) in a cell population where all cells are assumed to be susceptible. As in epidemiology, it can be proven that if $\mathcal{R}_0 < 1$ the disease-free equilibrium $P_0$ is locally asymptotically stable, this means that the infection is acute, otherwise if $\mathcal{R}_0 > 1$, $P_0$ becomes unstable and the infection is chronic.

## 4.2 Local and global asymptotic stability of the disease-free equilibrium

For local asymptotic stability and instability of the disease-free equilibrium $P_0$, only the conditions $\mathcal{R}_0 < 1$ and $\mathcal{R}_0 > 1$ are required.

**Theorem 4.1.** The disease-free equilibrium point $P_0$ is locally asymptotically stable if $\mathcal{R}_0 < 1$ and it is unstable if $\mathcal{R}_0 > 1$.

*Proof.* The proof is given in the appendix (Subsection AA.1). □

To analyze the global asymptotic stability of the disease-free steady-state $P_0$, we introduce the following subset $\mathcal{D}$ of $\mathbb{R}_+^2$,

$$ \mathcal{D} := \left\{ (A_1, A_2) \in \mathbb{R}_+^2 : A_1 \ge \frac{p_1(0)}{\gamma_{A_1}}, A_2 \ge \frac{p_2(0)}{\gamma_{A_2}} \right\}. $$

To prove the global asymptotic stability of $P_0$, we have to check two cases, whether the system starts with a nonnegative initial condition inside or outside the set $\mathcal{D} \times \mathbb{R}_+^3$. Before that, we need the following invariance result.

**Lemma 4.2.** Assume (H1). Then, the subset $\mathcal{D} \times \mathbb{R}_+^3$ is invariant under System (8).

*Proof.* The proof is given in the appendix (Subsection AA.3). □

We use a comparison result, Lemma 1 of [28] and [45], to prove the global asymptotic stability of $P_0$.

**Theorem 4.3.** Assume (H1), (H3), (H4), (H5) and $\mathcal{R}_0 < 1$. Then, the disease-free equilibrium $P_0$ is globally asymptotically stable.

*Proof.* The proof is given in the appendix (Subsection AA.1). □

**Remark 1.** We give here an example of a function $G$ that satisfies Condition (H3). We consider the case of a purely independent binding (6) and we assume that one of the following points is satisfied.

* $\xi_{1,i} \le \xi_{1,i-1} \le 1$ and $\xi_{2,j} \le \xi_{2,j-1} \le 1$, for all $2 \le i \le n_1$, $2 \le j \le n_2$,
• $\xi_{1,n_1} \le \xi_{1,n_1-1}$, $\xi_{2,n_2} \le \xi_{2,n_2-1}$, and all others equal to 1,

$$ \frac{p_1(0)}{\gamma_{A_1}} \ge (n_1 - 1) \frac{\xi_{1,n_1-1} - 1}{\xi_{1,n_1-1} - \xi_{1,n_1}} $$

and

$$ \frac{p_2(0)}{\gamma_{A_2}} \ge (n_2 - 1) \frac{\xi_{2,n_2-1} - 1}{\xi_{2,n_2-1} - \xi_{2,n_2}}. $$

Then, Condition (H3) is satisfied.

Indeed, in the case of a purely independent binding, we have $G(A_1, A_2) = G(A_1, 0) \times G(0, A_2)$ and if one of the two points of Remark 1 is satisfied, then the two functions $A_1 \mapsto G(A_1, 0)$ and $A_2 \mapsto G(0, A_2)$ are decreasing on $\mathbb{D}$. This means that Condition (H3) is satisfied on $\mathbb{D}$. Let's give a biological interpretation of the two conditions in Remark 1. Suppose that $n_1 = n_2 = 3$. We only need to explain for $G(A_1, 0)$ as the other is similar. The first condition becomes $\xi_{1,3} \le \xi_{1,2} \le \xi_{1,1} \le 1$. This means that the increase in the number of epitopes recognized and bound by antibodies promotes the decrease of the virus activity. In any case, the formed immune complexes result in the neutralization of virus activity. In this case, the function $A_1 \mapsto G(A_1, 0)$ is always decreasing. The second condition becomes $\xi_{1,3} \le \xi_{1,2}$, $\xi_{1,1} = 1$ and $\frac{p_1(0)}{\gamma_{A_1}} \ge 2 \frac{\xi_{1,2} - 1}{\xi_{1,2} - \xi_{1,3}}$. With this condition, even if the relative activities are not less than 1, the associated viral activity $A_1 \mapsto G(A_1, 0)$ is decreasing on the set $\mathbb{D}$ because the amount of antibodies before the virus invasion is big enough to stop the infection.

## 4.3 Existence of endemic equilibrium points

**Theorem 4.4.** Assume that $\mathcal{R}_0 > 1$. Then, there exists (at least) one endemic steady-state and all endemic steady-states belong to the subset $\mathbb{D} \times \mathbb{R}_+^3$.

*Proof.* The proof is given in the appendix (Subsection AA.1). □

As mentioned before, Condition (H3) promotes neutralization and it is not surprising in this case to have no endemic equilibrium if $\mathcal{R}_0 < 1$. We prove this in the following result.

**Proposition 2.** Assume that $\mathcal{R}_0 < 1$ and the function $G$ satisfies Condition (H3) on the subset $\mathbb{D}$. Then, there is no endemic steady-state.

*Proof.* The proof is given in the appendix (Subsection AA.2). □

**Proposition 3.** Suppose that the function $Y^* \mapsto G(p_1(Y^*)/\gamma_{A_1}, p_2(Y^*)/\gamma_{A_2})$ is decreasing on the interval $[0, \Omega/\gamma_Y]$. Then, if $\mathcal{R}_0 < 1$ there is no endemic steady-state, otherwise if $\mathcal{R}_0 > 1$ there is a unique endemic steady-state.

*Proof.* The proof is given in the appendix (Subsection AA.2). □

**Remark 2.** Suppose that the function $Y^* \mapsto G(p_1(Y^*)/\gamma_{A_1}, p_2(Y^*)/\gamma_{A_2})$ is not decreasing on the interval $[0, \Omega/\gamma_Y]$. Then, we may have the existence of:

• an endemic equilibrium, even in the case where $\mathcal{R}_0 < 1$,

• more than one endemic equilibrium in the case where $\mathcal{R}_0 > 1$.

## 4.4 Local asymptotic stability analysis of an endemic equilibrium point

**Theorem 4.5.** Assume that (H1), (H2), (H3), (H4) are satisfied, $\mathcal{R}_0 > 1$, and the following hypothesis

$$ p_1'(Y^*)G_{A_1}(A_1^*, A_2^*) + p_2'(Y^*)G_{A_2}(A_1^*, A_2^*) \le 0, \quad (11) $$

where $P^* = (A_1^*, A_2^*, X^*, Y^*, V^*)$ is an endemic steady-state, $G_{A_1}$ and $G_{A_2}$ are the partial derivatives of the function $G$ with respect to the first and second variables. Then, $P^* = (A_1^*, A_2^*, X^*, Y^*, V^*)$ is locally asymptotically stable.

*Proof.* The proof is given in the appendix (Subsection AA.1). □

**Remark 3.** Condition (11) can be interpreted geometrically by writing it as a scalar product $\nabla p \cdot (\nabla G)^T \le 0$ at the endemic equilibrium, where $p$ is the vector $(p_1, p_2)^T$, $T$ is for transpose and $\nabla$ is the gradient. As $p_i'(Y^*)$ is nonnegative, the vector $\nabla p$ is located in the first quadrant, $(+, +)$. Then, $(\nabla G)^T$ has to be at least in the second or fourth quadrant, $(-, +)$ or $(+, -)$. If $(\nabla G)^T$ is in the third quadrant, $(-, -)$, Condition (11) is always satisfied. It should also be pointed out that Condition (11) is only sufficient but not necessary for the result of Theorem 4.5 to be valid.

## 5 Numerical simulations

We will now numerically study the system (8) in two different scenarios: (i) a mixture of two neutralizing antibodies (neutralizing/neutralizing); (ii) a mixture of enhancing and neutralizing antibodies (enhancing/neutralizing). Our main goal is to investigate how competition/cooperation between the two antibodies $A_1$ and $A_2$ can affect disease progression during secondary infection. Competition or cooperation between the two antibodies is primarily determined by their ability to bind to the same receptor. This is indicated by the parameters $f_{ij}$. Once binding to the same receptor, the antibodies can interact with synergistic binding through the parameters $\tilde{\nu}_{ij}$ and $\tilde{\kappa}_{ij}$.

For the numerical simulations, we use the functions $p_i(Y) = \Lambda_i + \alpha_i Y$. The baseline parameters are taken from the literature (see [27, 28], for further details) and are summarized in Table 1. They are associated with
dengue infection. The number of parameters $f_{ij}$ is 9. In our analysis, we will consider the case where $f_{ij}$ are independent of $i, j$. The initial conditions are $(A_1(0), A_2(0), X(0), Y(0), V(0)) = (0, 0, 10^7, 0, 1)$. The units for $A_1, A_2$ are mol ml$^{-1}$, for $X, Y$ cells ml$^{-1}$, and for $V$ RNA copies ml$^{-1}$.

In our approach, we fix the two parameters $\tilde{\nu}_{ij}$ and $\tilde{\kappa}_{ij}$ at a value corresponding to each of the two scenarios, neutralizing/neutralizing and enhancing/neutralizing, and vary $0 \le f_{ij} \le 1$. Similar results can be obtained if we fix $f_{ij}$ and either $\tilde{\kappa}_{ij}$ or $\tilde{\nu}_{ij}$, and vary the other.

## 5.1 Homologous sequential DENV infection (Neutralizing/Neutralizing)

We consider the case where both antibodies $A_1$ and $A_2$ have the same characteristics to neutralize the virus, $K_{D_1} = K_{D_2}$ and $(\xi_{11}, \xi_{12}, \xi_{13}) = (\xi_{21}, \xi_{22}, \xi_{23})$. We choose $K_{D_1} = 0.6$ and $(\xi_{11}, \xi_{12}, \xi_{13}) = (1, 0.95, 0.8)$. We assume that their neutralizing effects are reinforced by the same synergistic binding, $\tilde{K}_{D_1} = \tilde{K}_{D_2}$ and $(\tilde{\xi}_{11}, \tilde{\xi}_{12}, \tilde{\xi}_{13}) = (\tilde{\xi}_{21}, \tilde{\xi}_{22}, \tilde{\xi}_{23})$. We take the values $\tilde{K}_{D_1} = 0.8$ and $(\tilde{\xi}_{11}, \tilde{\xi}_{12}, \tilde{\xi}_{13}) = (\sqrt{0.65}, 0.95 \times \sqrt{0.65}, 0.8\sqrt{0.65})$. This corresponds to $\tilde{\nu}_{ij} = 0.65 < 1$ and $\tilde{\kappa}_{ij} = 1$, for all $i, j = 1, 2, 3$.

We first consider the case where there is only one neutralizing antibody $A_2$, taking in Figure 3(a), $\Lambda_1 = 0$ and $\alpha_1 = 0$. We then examine what happens when a second neutralizing antibody is introduced, taking in Figure 3, (b)-(d), $\Lambda_1 = 7.5 \times 10^3$ and $\alpha_1 = 0.4$. Finally, we vary the parameter $f_{ij}$, which simulates the fraction of simultaneous binding of the two antibodies to the same receptor, in each sub-figure (b), (c), and (d) of Figure 3, by taking: (b) $f_{ij} = 0$, (c) $f_{ij} = 10^{-6}$ and (d) $f_{ij} = 1$. Figures 3, (a)-(d), show that the presence of two neutralizing antibodies is at least more efficient than a single one. Figure 3(b) shows that purely competitive binding, $f_{ij} = 0$, of two neutralizing antibodies, while causing a very rapid decrease in viral activity, does not significantly increase neutralization, compared with a single neutralizing antibody, Figures 3(a). As the synergistic binding is increased, neutralization is improved. Indeed, as the parameter $f_{ij}$ increases, the proportion of infected cells is reduced and the amount of virus at peak is lower, while infection also occurs later.

## 5.2 Heterologous sequential DENV infection (Enhancing/Neutralizing)

We consider one enhancing antibody, $A_1$, and one neutralizing antibody, $A_2$. We choose $K_{D_1} = \tilde{K}_{D_1} = K_{D_2} = \tilde{K}_{D_2} = 0.8$ and $(\xi_{11}, \xi_{12}, \xi_{13}) = (\tilde{\xi}_{11}, \tilde{\xi}_{12}, \tilde{\xi}_{13}) = (3, 3, 3)$. We assume that the relative viral activities associated with the second antibodies are such that $\xi_{1i}\xi_{2j} = \tilde{\xi}_{1i}\tilde{\xi}_{2j} = 1$. With this choice, we have $\tilde{\kappa}_{ij} = 1$ for all $i, j$. We also choose the parameter $\tilde{\nu}_{ij}$ with the

same value for all $i, j$. Two situations can be considered, one with $\tilde{\nu}_{ij} < 1$, meaning that the sum of the effects of the two antibodies leads to neutralization, and the other with $\tilde{\nu}_{ij} > 1$, meaning that the sum of the effects of the two antibodies produces enhancement. We focus here on the case with $\tilde{\nu}_{ij} = 0.4 < 1$.

Let's first consider the case of a single neutralizing antibody $A_2$, by taking in Figure 4(a), $\Lambda_1 = 0$ and $\alpha_1 = 0$. The parameters were chosen to obtain a high degree of neutralization. We then present the enhancing antibody $A_1$, by considering in Figure 4, (b)-(d), $\Lambda_1 = 7.5 \times 10^3$ and $\alpha_1 = 0.4$. As in the previous case, we study the effect of synergistic binding by increasing the parameter $f_{ij}$. More precisely, for $f_{ij} = 0$, a purely competitive binding, Figure 4(b) shows that the antibody $A_1$ strongly increases infection, with a peak of the virus reached very quickly and viral activity remaining in the red zone (refereed to Figure 2). As synergistic binding increases, Figure 4: (b) $f_{ij} = 0$, (c) $f_{ij} = 10^{-3}$, (d) $f_{ij} = 1$, neutralization becomes more effective than enhancement. Indeed, the proportion of infected cells is reduced and the amount of virus at peak is lower, while infection also occurs later. In the same way, viral activity falls into the blue zone (referred to Figure 2). Based on these comparisons, we can affirm that purely competitive binding, $f_{ij} = 0$, gives an advantage to the enhancing antibody while purely independent binding, $f_{ij} = 1$, confers a greater competitive advantage to the neutralizing antibody.

## 5.3 Fitting dengue data

Here we compare the numerical simulations with data from measurements of virus and antibody titers collected sequentially during infection of dengue patients, Figure 5. These data were previously published in [4]. We selected part of the data corresponding to secondary infection by DENV-2 and, among these, we considered only eight individuals for whom the size of the temporal series was long enough to permit estimating the parameters related to the force of infection. To compare the data with the numerical simulations, the viral load data are offset by an elapsed time $t_V$ which corresponds to the incubation period of the dengue virus of 2-6 days, [4, 46]. IgG and IgM antibody levels are multiplied by scaling factors, $SF_{IgG}$ and $SF_{IgM}$, because the ELISA assay cannot measure them directly, [4]. In addition, the IgG/IgM measurements are shifted by $t_{IgG}$ and $t_{IgM}$ to take into account the time delay between exposure to the antigen and the development of an effective immune response leading to clearance of the virus, [47]. As this is a secondary, the IgG titers correspond to the antibodies generated by the first infection and the IgM titers to those generated by the second infection. The scaling factors $SF_{IgG}$, $SF_{IgM}$, and the antibody offsets $t_{IgG}$ and $t_{IgM}$ are defined by hand. Table 1 shows the values of the parameters used in the simulations. Four parameters of the model $f_{ij}$, $\tilde{\nu}_{ij}$, $\delta$,
<figure><img src="image_4.png" /><figcaption>(b) Two homologous antibodies with purely competitive binding ($f_{ij} = 0$)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Two homologous antibodies with low synergistic binding ($f_{ij} = 10^{-6}$)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) Two homologous antibodies with high synergistic binding ($f_{ij} = 1$)</figcaption></figure>

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

Figure 3: Neutralizing/Neutralizing. Here we compare the temporal dynamics of antibodies, cells, and viral populations during a primary viral infection and a secondary infection by a homologous virus, by varying the competition between the two antibodies to bind to the virus, (b) $f_{ij} = 0$, (c) $f_{ij} = 10^{-6}$ and (d) $f_{ij} = 1$. From left to right, the first column of each panel shows the temporal evolution of the neutralizing antibody $A_1$ (dotted blue line), the neutralizing antibody $A_2$ (solid blue line), and the virus $V$ (solid red line). The second column shows the evolution of the proportion of susceptible and infected cells $X$ and $Y$ (solid blue and red lines). The third column shows the temporal activity of the virus (solid blue line). The gray dotted line indicates the threshold of 1, i.e., when the complex antigen-antibody does not affect virus activity. The last column shows the force of infection, FOI (red line). The parameter values are given in Table 1.
<figure><img src="image_8.png" /><figcaption>(b) Two heterologous antibodies with purely competitive binding ($f_{ij} = 0$)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Two heterologous antibodies with low synergistic binding ($f_{ij} = 10^{-3}$)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Two heterologous antibodies with high synergistic binding ($f_{ij} = 1$)</figcaption></figure>

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

Figure 4: Enhancing/Neutralizing. Here we compare the temporal dynamics of antibodies, cells, and viral populations during a primary viral infection and a secondary infection by a heterologous virus, by varying the competition between the two antibodies to bind to the virus, (b) $f_{ij} = 0$, (c) $f_{ij} = 10^{-3}$ and (d) $f_{ij} = 1$. From left to right, the first column of each panel shows the temporal evolution of the enhancing antibody $A_1$ (dotted red line), the neutralizing antibody $A_2$ (solid blue line), and the virus $V$ (solid red line). The second column shows the evolution of the proportion of susceptible and infected cells $X$ and $Y$ (solid blue and red lines). The gray dotted line indicates the threshold of 1, i.e., when the complex antigen-antibody does not affect virus activity. The third column shows the temporal virus activity (solid blue-red line). The last column shows the force of infection, FOI (red line). The parameter values are given in Table 1.
<figure><img src="image_12.png" /><figcaption>Figure 5: Temporal evolution of enhancing antibodies (red curves), neutralizing antibodies (blue curves), and viral populations (black curves) during secondary infection with a heterologous virus. Each plot corresponds to an individual infected with dengue virus serotype II. Patient data are represented by symbols: × IgG titers, □ IgM titers and ★ virus titers. The set of parameters used in each simulation is given in Tables 2 and 3, and the patient identification number on the top right of each plot.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 6: Boxplots of estimated model parameters *f*<sub>ij</sub>, *ν̃*<sub>ij</sub>, *δ*, *γ*<sub>Y</sub> and virus offset *t*<sub>V</sub> for each patient. The results comprise a set of 200 parameters obtained after running the Genetic Algorithm (GA) for 50 generations. The set of parameters corresponding to the best solution of the GA is represented by the symbol ×, and the corresponding values are given in Table 2.</figcaption></figure>
$\gamma_Y$ plus the virus offset $t_V$ are estimated by a Genetic Algorithm (GA). For each individual patient, the best solution corresponds to the minimum of the function $\sum_{i=1}^{n_m} |(\log_{10}(V(t_i) + 10) - \log_{10}(V_i + 10))| / \log_{10}(\sigma_V^2 + 10)$, where $\sigma_V^2$ is the variance of the viral load measurements, $V_i$ the $i^{th}$ measurement, $n_m$ the number of measurements for the patient, $V(t_i)$ the virus load given by the simulation at time $t_i$. The results are plotted in Figure 5. The parameter ranges used in the GA are as follows $f_{ij} \in [0, 1]$, $\tilde{v}_{ij} \in [0, 1]$, $\delta \in [0, 5]$, $\gamma_Y \in [2.5, 3]$ and $t_V \in [0, 10]$. In addition, the other parameters used to run the GA are 200, 0.1, respectively, the size of the GA population and the crossover rate. The mutation function randomly generates directions that are adaptive with respect to the last successful or unsuccessful generation. The code was run in Matlab. Tables 2, 3 and Figure 6 summarize the scaling factors, antibody shifts, and estimated model parameters as well as virus offsets. The main objective is either to validate the model and to discuss the heterogeneity of patients' viral kinetics. We note that our curves have the same shape as the data, Figure 5. In addition, parameter estimation provides very low values of $f_{ij}$, Figure 6 and Tables 2, 3, showing that the two antibodies are in very strong competition (as the fraction of simultaneous binding is very low), which contributes to high virus levels. Whereas for patient 102, the values of $\tilde{v}_{ij}$ range from 0 to 0.5, for patient 233 they range from 0.6 to 1. In both cases, the relative viral activity due to synergistic binding is low, but more so for the first patient. In particular, the set of parameters corresponding to the best solution of the GA is represented by the symbol $\times$ in Figure 6. The model assumes a single immunity variable that starts immediately after virus exposure. But, natural infection by dengue virus triggers a fast antiviral innate immune response and a later adaptive immune response that resolves infection and leads to long-term immunological memory.

# 6 Discussion

We developed a formula that gives the virus activity when bound to two antibodies, using the principle of chemical reactions. It generalizes the work of [1] by assuming that the antibodies can bind to multiple sites on a receptor. We have thus been able to show that this antigen-antibody complex can neutralize or enhance virus activity, given the interaction - competition or cooperation - between antibodies. The viral activity was then integrated into a model that emulates the interaction between healthy target cells, infected cells, viruses, and antibodies. This allows us to study the dynamic process of immune response when challenged by sequential homologous or heterologous viruses. We highlighted two scenarios: a mixture of two neutralizing antibodies (Figure 3) and a mixture of enhancing and neutralizing antibodies (Figure 4). Other situa-

tions not considered in this study, such as the mixture of two types of enhancing antibodies, could be envisaged. This can happen, for example, when using therapeutic antibodies. In this case, $p_i$, $i = 1, 2$, can also be a function of time. We could also generalize the model to the case where there are more than two antibody types, as in polyclonal antibody therapy, but the number of parameters would increase rapidly, making parametrization, analysis, and interpretation of the model a challenge. By keeping the model as simple as possible, we were able to explore it analytically and numerically and assert that virus neutralization is associated with purely independent binding (Figures 3(d) and 4(d)), while enhanced virus activity is expected when purely competitive binding occurs (Figure 4(b)). This opens up promising modeling prospects for the study of second viral infections, vaccine efficacy, or antibody treatment. Other potential applications, such as transcriptional gene regulation, [48], or any field where there is an interaction between activation and inhibition processes, could be favored by this type of new modeling. If the condition (H3) and/or the assumption in Proposition 3 are relaxed, then we could have several endemic equilibrium points even in the parameter space where $R_0 < 1$. It would be interesting to study their biological significance and asymptotic behavior, given that the latter depends on initial conditions such as the size of the virus inoculum and the amount of each antibody present in the mixture. Another challenge for our model is to confront it with biological data for dengue, influenza, Covid-19, or HIV. Preliminary tests carried out here on data from dengue patients are promising. The data concern secondary heterologous infection, and cover 5 or 6 days after the onset of symptoms; they therefore start around the time when virus titers reach their peak. For each patient, the model captures the range of virus titers and their rapid decrease after the peak. Rapid increases in IgG and IgM titers are also observed (Figure 5). Furthermore, for the set of data explored here, the antigen-antibody complex formed includes antibodies with a very low synergistic neutralizing effect (Figure 6). This corresponds to high viral activity. Finally, the estimated parameters show that the mortality of infected cells $\gamma_Y$ is almost twice as high as that of susceptible cells and that the heterogeneity of patients' viral kinetics is associated with variability of antibody responses between individuals, measured by the parameters $f_{ij}$, $\tilde{v}_{ij}$ and $\delta$ (Table 3). This is in agreement with previous works showing that antibody kinetics shape virus dynamics acting either on free virus or infected cells [4]. Other mechanisms not considered here, such as cell-mediated immune responses, may also be associated with the lysis of virus-infected cells. In addition, although there are only four antigenically distinct serotypes of the dengue virus, each serotype is divided into several genotypes. This factor, combined with genetics and the immune status of the host, may explain the kinetics of viral load and antibodies, [46].
<table><caption>Table 1: Parameter of the mathematical model with their descriptions, values and units. Some parameters are varied among simulations and are distinguished: Figure 3 □, Figure 4 ⊗, and Figure 5 □. The baseline parameters are taken from the literature (see [27, 28], for further details).</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value and Units</th></tr></thead><tbody><tr><td>n<sub>1</sub>, n<sub>2</sub></td><td>Number of epitopes that can be recognized by A<sub>1</sub>, A<sub>2</sub></td><td>3</td></tr><tr><td>Ω</td><td>Production rate of susceptible target cells</td><td>6.7 × 10<sup>3</sup> cells ml<sup>-1</sup> days<sup>-1</sup></td></tr><tr><td>Φ</td><td>Production rate of virus by infected cells</td><td>1 × 10<sup>4</sup> RNA copies cells<sup>-1</sup> days<sup>-1</sup></td></tr><tr><td>β</td><td>Virus transmission rate</td><td>8.5 × 10<sup>-10</sup> ml RNA copies<sup>-1</sup> days<sup>-1</sup></td></tr><tr><td>γ<sub>X</sub></td><td>Natural mortality rate of target cells</td><td>log(2)/5 days<sup>-1</sup></td></tr><tr><td>γ<sub>Y</sub></td><td>Natural mortality rate of infected cells</td><td>log(2) days<sup>-1</sup> ; □</td></tr><tr><td>γ<sub>V</sub></td><td>Natural mortality rate of virus</td><td>24 log(2)/2.7 days<sup>-1</sup></td></tr><tr><td>γ<sub>A<sub>1</sub></sub> γ<sub>A<sub>2</sub></sub></td><td>Natural mortality rate of antibodies</td><td>log(2)/40 days<sup>-1</sup></td></tr><tr><td>α<sub>1</sub>, α<sub>2</sub></td><td>Production rate of antibodies A<sub>1</sub>, A<sub>2</sub> by infected cells</td><td>0.4, 0.44 mol cells<sup>-1</sup> days<sup>-1</sup></td></tr><tr><td>Λ<sub>1</sub>, Λ<sub>2</sub></td><td>Natural production rate of antibodies A<sub>1</sub>, A<sub>2</sub></td><td>7.5 × 10<sup>3</sup>, 0 mol ml<sup>-1</sup> days<sup>-1</sup></td></tr><tr><td>(ξ<sub>11</sub>, ξ<sub>12</sub>, ξ<sub>13</sub>)</td><td>Relative activity of the virus-antibody A<sub>1</sub> complex</td><td>(1, 0.95, 0.8) □; (3, 3, 3) ⊗</td></tr><tr><td>(ξ<sub>21</sub>, ξ<sub>22</sub>, ξ<sub>23</sub>)</td><td>Relative activity of the virus-antibody A<sub>2</sub> complex</td><td>(1, 0.95, 0.8) □; (1/3, 1/3, 1/3) ⊗</td></tr><tr><td>ν̃<sub>ij</sub></td><td>Ratio between modified and normal relative activities of the virus-antibody complex</td><td>0.65 □; 0.4 ⊗ ; □</td></tr><tr><td>K<sub>D<sub>1</sub></sub>, K<sub>D<sub>2</sub></sub></td><td>Dissociation constants</td><td>0.8</td></tr><tr><td>κ<sub>ij</sub></td><td>Ratio between modified and normal dissociation constants</td><td>1</td></tr><tr><td>δ</td><td>Rate of free virus loss by other means than natural mortality</td><td>0.8 days<sup>-1</sup> ; □</td></tr><tr><td>f<sub>ij</sub></td><td>Fraction of simultaneous binding of antibodies A<sub>1</sub>, A<sub>2</sub></td><td>0, 10<sup>-6</sup>, 1 □; 0, 10<sup>-3</sup>, 1 ⊗ ; □</td></tr></tbody></table>
<table><caption>Table 2: Fraction of simultaneous binding $f_{ij}$, ratio between modified and normal relative activities $\tilde{v}_{ij}$, rate of loss of free viruses by means other than natural mortality $\delta$, natural mortality of infected cells $\gamma_Y$, time offset values $t_V$, $t_{IgG}$ and $t_{IgM}$, and scaling factors $SF_{IgG}$ and $F_{IgM}$, for each individual, used for Figure 5. The values of $f_{ij}$, $\tilde{v}_{ij}$, $\delta$, $\gamma_Y$ and $t_V$ are the results of the Genetic Algorithm.</caption><thead><tr><th>Patient ID</th><th>f<sub>ij</sub></th><th>ṽ<sub>ij</sub></th><th>δ</th><th>γ<sub>Y</sub></th><th>t<sub>V</sub></th><th>SF<sub>IgG</sub></th><th>t<sub>IgG</sub></th><th>SF<sub>IgM</sub></th><th>t<sub>IgM</sub></th></tr></thead><tbody><tr><td>102</td><td>0</td><td>0.47</td><td>2.13</td><td>2.87</td><td>5.2</td><td>1.83 × 10<sup>4</sup></td><td>1</td><td>3.88 × 10<sup>4</sup></td><td>0.5</td></tr><tr><td>105</td><td>8.01 × 10<sup>-19</sup></td><td>1</td><td>0.69</td><td>2.73</td><td>4.77</td><td>2.00 × 10<sup>4</sup></td><td>-0.5</td><td>2.36 × 10<sup>4</sup></td><td>-0.25</td></tr><tr><td>107</td><td>5.65 × 10<sup>-20</sup></td><td>0.43</td><td>1.29</td><td>2.62</td><td>5.04</td><td>1.75 × 10<sup>4</sup></td><td>0</td><td>8.96 × 10<sup>4</sup></td><td>0.5</td></tr><tr><td>115</td><td>1.74 × 10<sup>-16</sup></td><td>0.40</td><td>0.61</td><td>2.75</td><td>5.31</td><td>1.92 × 10<sup>4</sup></td><td>-2</td><td>1.62 × 10<sup>4</sup></td><td>-1.5</td></tr><tr><td>207</td><td>2.68 × 10<sup>-19</sup></td><td>0.61</td><td>0.11</td><td>2.86</td><td>3.35</td><td>3.07 × 10<sup>4</sup></td><td>-0.5</td><td>2.00 × 10<sup>4</sup></td><td>1</td></tr><tr><td>217</td><td>1.28 × 10<sup>-16</sup></td><td>0.83</td><td>0.45</td><td>2.95</td><td>4.22</td><td>2.50 × 10<sup>4</sup></td><td>0</td><td>2.64 × 10<sup>4</sup></td><td>0</td></tr><tr><td>224</td><td>5.08 × 10<sup>-19</sup></td><td>0.46</td><td>0.96</td><td>2.98</td><td>4.66</td><td>2.23 × 10<sup>4</sup></td><td>-1</td><td>2.46 × 10<sup>4</sup></td><td>0</td></tr><tr><td>233</td><td>2.68 × 10<sup>-19</sup></td><td>0.93</td><td>0.09</td><td>2.86</td><td>2.28</td><td>2.39 × 10<sup>4</sup></td><td>-1</td><td>3.46 × 10<sup>4</sup></td><td>-1</td></tr></tbody></table>

<table><caption>Table 3: Summary of the parameter values obtained by the Genetic Algorithm (GA) (see Table 2). In addition to parameter descriptions, we provide median, minimum and maximum values.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Median [min, max]</th></tr></thead><tbody><tr><td>f<sub>ij</sub></td><td>Fraction of simultaneous binding</td><td>3.88 × 10<sup>-19</sup>[0, 1.74 × 10<sup>-16</sup>]</td></tr><tr><td>ṽ<sub>ij</sub></td><td>Ratio between modified and normal relative activities</td><td>0.54 [0.39, 1]</td></tr><tr><td>δ</td><td>Rate of loss of free viruses by means other than natural mortality</td><td>0.65 [0.0882, 2.13] days<sup>-1</sup></td></tr><tr><td>γ<sub>Y</sub></td><td>Natural mortality of infected cells</td><td>2.86 [2.61, 2.98] days<sup>-1</sup></td></tr><tr><td>t<sub>V</sub></td><td>Time between onset of infection and first sampling</td><td>4.71 [2.27, 5.31] days</td></tr></tbody></table>

# A Appendix

## A.1 Proofs of Theorems

*Proof of Theorem 4.1.* Let $J(P_0)$ be the Jacobian matrix associated with the system (8), evaluated at the equilibrium point $P_0$. Then, we have

$$J(P_0) = \begin{pmatrix} -\gamma_{A_1} & 0 & 0 & p'_1(0) & 0 \\ 0 & -\gamma_{A_2} & 0 & p'_2(0) & 0 \\ 0 & 0 & -\gamma_X & 0 & -\beta G^* X^* \\ 0 & 0 & 0 & -\gamma_Y & \beta G^* X^* \\ 0 & 0 & 0 & 0 & -\Phi - (\gamma_V + \delta) \end{pmatrix},$$

with $A_1^* = \frac{p_1(0)}{\gamma_{A_1}}$, $A_2^* = \frac{p_2(0)}{\gamma_{A_2}}$, $X^* = \frac{\Omega}{\gamma_X}$ and $G^* = G(A_1^*, A_2^*)$. Then, the characteristic equation of $J(P_0)$ is

$$\left| \begin{array}{ccccc} \lambda + \gamma_{A_1} & 0 & 0 & -p'_1(0) & 0 \\ 0 & \lambda + \gamma_{A_2} & 0 & -p'_2(0) & 0 \\ 0 & 0 & \lambda + \gamma_X & 0 & \beta G^* X^* \\ 0 & 0 & 0 & \lambda + \gamma_Y & -\beta G^* X^* \\ 0 & 0 & 0 & -\Phi & \lambda + (\gamma_V + \delta) \end{array} \right| = 0.$$

The eigenvalues of $J(P_0)$ are $\lambda = -\gamma_{A_1}$, $\lambda = -\gamma_{A_2}$, $\lambda = -\gamma_X$ and the roots of the polynomial

$$(\lambda + \gamma_Y)(\lambda + \gamma_V + \delta) - \Phi\beta G^* X^* = 0.$$

We obtain the following polynomial

$$\lambda^2 + (\gamma_V + \delta + \gamma_Y)\lambda + \gamma_Y(\gamma_V + \delta) - \Phi\beta G^* X^* = 0. \qquad (12)$$
According to the Routh-Hurwitz criterion, a polynomial of degree two has roots with a negative real part if and only if the coefficients $a_1$ and $a_2$ of the polynomial $\lambda^2 + a_1\lambda + a_2$ are positive. For the polynomial (12), we have

$$ 
\begin{aligned}
a_1 &= \gamma_V + \delta + \gamma_Y > 0, \\
a_2 &= \gamma_Y (\gamma_V + \delta) - \Phi \beta G^* X^*, \\
    &= \gamma_Y (\gamma_V + \delta) (1 - \mathcal{R}_0).
\end{aligned}
$$

Then, $a_2 > 0$ if and only if $\mathcal{R}_0 < 1$. Therefore, the roots of the polynomial (12) have negative real parts if and only if $\mathcal{R}_0 < 1$. We conclude that the disease-free equilibrium $P_0$ is locally asymptotically stable if $\mathcal{R}_0 < 1$, and unstable if $\mathcal{R}_0 > 1$. $\square$

**Proof of Theorem 4.3.** Again, we take System (22). Let first consider the case $\mu = -1$. Then, we have $\frac{\Omega}{\gamma} + \mu Z = X + Y \le \frac{\Omega}{\gamma}$. Therefore, System (22) satisfies the following inequalities on $\mathbb{D} \times \mathbb{R}_+^3$

$$ 
\left\{
\begin{aligned}
\frac{dB_1}{dt} &\le p'_1(0)Y - \gamma_{A_1}B_1, \\
\frac{dB_2}{dt} &\le p'_2(0)Y - \gamma_{A_2}B_2, \\
\frac{dZ}{dt} &= -\gamma Z, \\
\frac{dY}{dt} &\le \beta G \left( \frac{p_1(0)}{\gamma_{A_1}}, \frac{p_2(0)}{\gamma_{A_2}} \right) \frac{\Omega}{\gamma}V - \gamma Y, \\
\frac{dV}{dt} &= \Phi Y - (\gamma_V + \delta)V.
\end{aligned}
\right.
$$

Then, System (22) can be compared to the following linear system

$$ 
\left\{
\begin{aligned}
\frac{dB_1}{dt} &= p'_1(0)Y - \gamma_{A_1}B_1, \\
\frac{dB_2}{dt} &= p'_2(0)Y - \gamma_{A_2}B_2, \\
\frac{dZ}{dt} &= -\gamma Z, \\
\frac{dY}{dt} &= \beta G \left( \frac{p_1(0)}{\gamma_{A_1}}, \frac{p_2(0)}{\gamma_{A_2}} \right) \frac{\Omega}{\gamma}V - \gamma Y, \\
\frac{dV}{dt} &= \Phi Y - (\gamma_V + \delta)V.
\end{aligned}
\right.
\qquad (13)
$$

This last system has the following characteristic equation

$$ 
\left|
\begin{array}{ccccc}
\lambda + \gamma_{A_1} & 0 & 0 & -p'_1(0) & 0 \\
0 & \lambda + \gamma_{A_2} & 0 & -p'_2(0) & 0 \\
0 & 0 & \lambda + \gamma & 0 & 0 \\
0 & 0 & 0 & \lambda + \gamma & S \\
0 & 0 & 0 & -\Phi & \lambda + \gamma_V + \delta
\end{array}
\right| = 0.
$$

with $S := -\beta G (p_1(0)/\gamma_{A_1}, p_2(0)/\gamma_{A_2}) \Omega / \gamma$. It is equivalent to

$$ (\lambda + \gamma_{A_1})(\lambda + \gamma_{A_2})(\lambda + \gamma) \times [\lambda^2 + \lambda(\gamma + \gamma_V + \delta) + \gamma(\gamma_V + \delta)(1 - \mathcal{R}_0)] = 0. \quad (14) $$

The eigenvalues are $\lambda = -\gamma_{A_1} < 0$, $\lambda = -\gamma_{A_2} < 0$, $\lambda = -\gamma < 0$ and the roots of the polynomial

$$ \lambda^2 + (\gamma + \gamma_V + \delta)\lambda + \gamma(\gamma_V + \delta)(1 - \mathcal{R}_0). $$

Using the Routh-Hurwitz criteria, we find that the roots of this last polynomial have negative real parts if and only if $\mathcal{R}_0 < 1$. Thus, we obtain the global asymptotic stability of the trivial equilibrium of the linear system (13). We now use the following comparison result to conclude (Lemma 1 of [28] and [45]).
**Lemma A.1.** Consider two differential systems $x' = h(x)$ and $y' = g(y)$ given on an invariant subset $U$ of $\mathbb{R}^k$, $h, g: U \to \mathbb{R}^k$ are locally Lipschitz functions. Then, the following two conditions are equivalent:

1. For each $x_0, y_0 \in U$ the inequality $x_0 \le y_0$ implies $x(t) \le y(t)$ for all $t \ge 0$, where $x'(t) = h(x(t))$ and $y'(t) = g(y(t))$, $t \ge 0$, with $x(0) = x_0$, $y(0) = y_0$.

2. For all $i = 1, \dots, k$, the inequality

$$h_i(x_1, \dots, x_{i-1}, x_i, x_{i+1}, \dots, x_n) \le g_i(\bar{x}_1, \dots, \bar{x}_{i-1}, \bar{x}_i, \bar{x}_{i+1}, \dots, \bar{x}_n)$$

holds whenever $x_j \le \bar{x}_j$, for all $j \ne i$ and $x_i = \bar{x}_i$.

As our corresponding functions $f$ and $g$ are of class $C^1$, they are locally Lipschitz and we can apply Lemma A.1. Therefore, if we chose the same initial condition on $\mathbb{D} \times \mathbb{R}_+^3$, we obtain a solution of System (22) between $0_{\mathbb{R}^5}$ and the solution of System (13). Hence, the disease-free steady-state $P_0$ is globally asymptotically stable on the subset $\mathbb{D} \times \mathbb{R}_+^3$.

We continue the proof of Theorem 4.3 by now considering the case $\mu = 1$. Using the equation of $Z$, we can show for all $\epsilon > 0$ and all $Z_0 > 0$, the existence of $T_\epsilon \ge 0$, we can take $T_\epsilon = \max\{0, \frac{1}{\gamma} \ln\left(\frac{Z_0}{\epsilon}\right)\}$, such that $0 < Z(t) < \epsilon$, for all $t > T_\epsilon$. We then compare System (22), for $t > T_\epsilon$, to the following linear system

$$\left\{
\begin{aligned}
\frac{dB_1}{dt} &= p'_1(0)Y - \gamma_{A_1}B_1, \\
\frac{dB_2}{dt} &= p'_2(0)Y - \gamma_{A_1}B_2, \\
\frac{dZ}{dt} &= -\gamma Z, \\
\frac{dY}{dt} &= \beta G \left( \frac{p_1(0)}{\gamma_{A_1}}, \frac{p_2(0)}{\gamma_{A_2}} \right) \left( \frac{\Omega}{\gamma} + \epsilon \right) V - \gamma Y, \\
\frac{dV}{dt} &= \Phi Y - (\gamma_V + \delta) V.
\end{aligned}
\right.
\qquad (15)$$

Its characteristic equation is given by

$$0 = (\lambda + \gamma_{A_1})(\lambda + \gamma_{A_2})(\lambda + \gamma) \times \left[ (\lambda^2 + \lambda(\gamma + \gamma_V) + \gamma\gamma_V(1 - \mathcal{R}_\epsilon)) \right], \qquad (16)$$

where

$$\mathcal{R}_\epsilon := \frac{\Phi\beta}{\gamma(\gamma_V + \delta)} \left(\frac{\Omega}{\gamma} + \epsilon\right) G \left(\frac{p_1(0)}{\gamma_{A_1}}, \frac{p_2(0)}{\gamma_{A_2}}\right).$$

Under the condition $\mathcal{R}_0 < 1$ and by choosing $\epsilon > 0$ small enough, we have $\mathcal{R}_\epsilon < 1$. Then, we obtain the global asymptotic stability of the trivial solution of the linear system (15). So, using again Lemma A.1, we conclude the global asymptotic stability of the disease-free steady-state $P_0$ on the set $\mathbb{D} \times \mathbb{R}_+^3$. So far, we have considered the solutions of System (8) that start in the set $\mathbb{D} \times \mathbb{R}_+^3$. Now consider a nonnegative solution such that $A_1(0) < \frac{p_1(0)}{\gamma_{A_1}}$ and/or $A_2(0) < \frac{p_2(0)}{\gamma_{A_2}}$, which means that it does not start in the set $\mathbb{D} \times \mathbb{R}_+^3$. We have two different cases.

Firstly, suppose that we have the existence of $\bar{t}_i > 0$, $i = 1, 2$, such that $A_i(\bar{t}_i) = \frac{p_i(0)}{\gamma_{A_i}}$ and $A_i(t) < \frac{p_i(0)}{\gamma_{A_i}}$ for all $t < \bar{t}_i$. From the equation of $A_i$ in System (8) and under (H4), we can see that $A'_i(t) \ge p_i(Y(t)) - p_i(0) \ge 0$, for $t \le \bar{t}_i$. This means that $A_i$ is an increasing function on $[0, \bar{t}_i]$. Hence, the solution enters the subset $\mathbb{D} \times \mathbb{R}_+^3$ at time $\bar{t} = \max_{i=1,2} A_i(t)$. As $\mathbb{D} \times \mathbb{R}_+^3$ is invariant for the system (8), the solution tends to $P_0$, because $P_0$ is globally asymptotically stable on the set $\mathbb{D} \times \mathbb{R}_+^3$.

Secondly, we suppose that there is no $\bar{t}_i$ as described above (for $i = 1$ or and $i = 2$). We consider that for all $t \ge 0$, $A_i(t) < \frac{p_i(0)}{\gamma_{A_i}}$. Then, $A'_i(t) > 0$. This means that $A_i$ is increasing on $[0, +\infty)$ and admits a limit $l_i := \lim_{t \to +\infty} A_i(t) \le \frac{p_i(0)}{\gamma_{A_i}}$. Furthermore, $\lim_{t \to +\infty} A'_i(t) = 0$. Then, we can see that

$$\lim_{t \to +\infty} p_i(Y(t)) = \gamma_{A_i} l_i \le p_i(0).$$

But, according (H5), $p_i$ ($i = 1$ or $i = 2$) is an increasing function in a neighbourhood of 0. So, we conclude that

$$\lim_{t \to +\infty} Y(t) = 0 \quad \text{and} \quad \lim_{t \to +\infty} A_i(t) = l_i = \frac{p_i(0)}{\gamma_{A_i}}.$$
We have shown that $\mathbb{D} \times \mathbb{R}_+^3$ is globally attractive for System (8). We can now conclude that the disease-free steady-state $P_0$ is globally asymptotically stable. $\square$

*Proof of Theorem 4.4.* An endemic equilibrium $P^* = (A_1^*, A_2^*, X^*, Y^*, V^*)$ satisfies

$$ 
\left\{ 
\begin{array}{ll}
\gamma_{A_1} A_1^* & = p_1(Y^*), \\
\gamma_{A_2} A_2^* & = p_2(Y^*), \\
\gamma_X X^* + \beta G(A_1^*, A_2^*) V^* X^* & = \Omega, \\
\beta G(A_1^*, A_2^*) V^* X^* & = \gamma_Y Y^*, \\
(\gamma_V + \delta) V^* & = \Phi Y^*.
\end{array}
\right.
\qquad (17) $$

Subtracting the third from the fourth equation of the (17) system, we get,

$$ \Omega - \gamma_X X^* - \gamma_Y Y^* = 0, $$

which results in

$$ X^* = \frac{\Omega - \gamma_Y Y^*}{\gamma_X}, \quad \text{with} \quad 0 < Y^* \le \frac{\Omega}{\gamma_Y}. \qquad (18) $$

From the first equation of the system, we obtain,

$$ A_1^* = \frac{1}{\gamma_{A_1}} p_1(Y^*). $$

From the second equation of the system (17), we get

$$ A_2^* = \frac{1}{\gamma_{A_2}} p_2(Y^*). $$

We can see that $(A_1^*, A_2^*) \in \mathbb{D}$, then $P^* \in \mathbb{D} \times \mathbb{R}_+^3$. From the last equation of System (17), we get

$$ V^* = \frac{\Phi Y^*}{\gamma_V + \delta}. \qquad (19) $$

Substituting the equations of $X^*$ and $V^*$ into the fourth equation of the system (17) we get

$$ \beta G(A_1^*, A_2^*) \left( \frac{\Phi Y^*}{\gamma_V + \delta} \right) \frac{\Omega - \gamma_Y Y^*}{\gamma_X} = \gamma_Y Y^*. $$

As we are looking for $Y^* > 0$, the previous equation is equivalent to

$$ \frac{\beta \Phi}{\gamma_Y \gamma_X} (\Omega - \gamma_Y Y^*) G(A_1^*, A_2^*) \left( \frac{1}{\gamma_V + \delta} \right) = 1 $$

We consider the function $I: [0, \frac{\Omega}{\gamma_Y}] \to [0, +\infty)$ defined by

$$ I(Y^*) = \frac{1}{\gamma_Y} \times \frac{\Omega - \gamma_Y Y^*}{\gamma_X} \times \frac{\Phi}{\gamma_V + \delta} \times \beta G \left( \frac{p_1(Y^*)}{\gamma_{A_1}}, \frac{p_2(Y^*)}{\gamma_{A_2}} \right). $$

The problem to solve is the following

$$ I(Y^*) = 1, \quad \text{with} \quad 0 \le Y^* \le \frac{\Omega}{\gamma_Y}. $$

We have

$$ I(0) = \mathcal{R}_0 \quad \text{and} \quad I\left(\frac{\Omega}{\gamma_Y}\right) = 0. $$

Therefore, if $\mathcal{R}_0 > 1$, there is at least one positive solution $Y^*$. $\square$
Proof of Theorem 4.5. Let $A_1(t) = A_1^* + a_1(t)$, $A_2(t) = A_2^* + a_2(t)$, $X(t) = X^* + x(t)$, $Y(t) = Y^* + y(t)$, $V(t) = V^* + v(t)$. The linearization of System (8) around $P^*$ is given by

$$ egin{cases} 
\frac{da_1}{dt} = p'_1(Y^*)y - \gamma_A a_1, \\
\frac{da_2}{dt} = p'_2(Y^*)y - \gamma_A a_2, \\
\frac{dx}{dt} = -\gamma x - \beta G(A_1^*, A_2^*) (V^*x + X^*v) \\
\qquad - \beta G_{A_1}(A_1^*, A_2^*) V^*X^*a_1 - \beta G_{A_2}(A_1^*, A_2^*) V^*X^*a_2, \\
\\
\frac{dy}{dt} = -\gamma y + \beta G(A_1^*, A_2^*) (V^*x + X^*v) \\
\qquad + \beta G_{A_1}(A_1^*, A_2^*) V^*X^*a_1 + \beta G_{A_2}(A_1^*, A_2^*) V^*X^*a_2, \\
\\
\frac{dv}{dt} = \Phi y - (\gamma_V + \delta)v.
\end{cases} 
\quad (20) $$

The characteristic equation associated with the previous system is

$$ \left| \begin{array}{cccccc} \lambda + \gamma_A & 0 & 0 & -\alpha_1 & 0 \\ 0 & \lambda + \gamma_A & 0 & -\alpha_2 & 0 \\ \beta_1 V^* X^* & \beta_2 V^* X^* & \lambda + \gamma + \beta_0 V^* & 0 & \beta_0 X^* \\ -\beta_1 V^* X^* & -\beta_2 V^* X^* & -\beta_0 V^* & \lambda + \gamma & -\beta_0 X^* \\ 0 & 0 & 0 & -\Phi & \lambda + \gamma_V + \delta \end{array} \right| = 0, $$

with

$$ \alpha_i = p'_i(Y^*) \ge 0, \quad \beta_i = \beta G_{A_i}(A_1^*, A_2^*), \quad i = 1, 2 \quad \text{and} \quad \beta_0 = \beta G(A_1^*, A_2^*) > 0. $$

We add the third and the fourth lines that we substitute to line 4 and the characteristic equation becomes

$$ \left| \begin{array}{cccccc} \lambda + \gamma_A & 0 & 0 & -\alpha_1 & 0 \\ 0 & \lambda + \gamma_A & 0 & -\alpha_2 & 0 \\ \beta_1 V^* X^* & \beta_2 V^* X^* & \lambda + \gamma + \beta_0 V^* & 0 & \beta_0 X^* \\ 0 & 0 & \lambda + \gamma & \lambda + \gamma & 0 \\ 0 & 0 & 0 & -\Phi & \lambda + \gamma_V + \delta \end{array} \right| = 0. $$

The subtraction of columns 3 and 4 that we substitute to column 4 gives,

$$ \left| \begin{array}{cccccc} \lambda + \gamma_A & 0 & 0 & -\alpha_1 & 0 \\ 0 & \lambda + \gamma_A & 0 & -\alpha_2 & 0 \\ \beta_1 V^* X^* & \beta_2 V^* X^* & \lambda + \gamma + \beta_0 V^* & -\lambda - \gamma - \beta_0 V^* & \beta_0 X^* \\ 0 & 0 & \lambda + \gamma & 0 & 0 \\ 0 & 0 & 0 & -\Phi & \lambda + \gamma_V + \delta \end{array} \right| = 0. $$

Developing firstly from the fourth line and secondly from the first column we obtain,

$$ -(\lambda + \gamma) ((\lambda + \gamma_A)D_1 + \beta_1 V^* X^* D_2) = 0, $$

with

$$ D_1 = \left| \begin{array}{ccc} \lambda + \gamma_A & -\alpha_2 & 0 \\ \beta_2 V^* X^* & -\lambda - \gamma - \beta_0 V^* & \beta_0 X^* \\ 0 & -\Phi & \lambda + \gamma_V + \delta \end{array} \right|, $$

and

$$ D_2 = \left| \begin{array}{ccc} 0 & -\alpha_1 & 0 \\ \lambda + \gamma_A & -\alpha_2 & 0 \\ 0 & -\Phi & \lambda + \gamma_V + \delta \end{array} \right|. $$

Hence we get,

$$ D_1 = (\lambda + \gamma_A) [(-\lambda - \gamma - \beta_0 V^*)(\lambda + \gamma_V + \delta) + \Phi \beta_0 X^*] + \alpha_2 \beta_2 V^* X^*(\lambda + \gamma_V + \delta), $$
and

$$D_2 = \alpha_1(\lambda + \gamma_A)(\lambda + \gamma_V + \delta).$$

Then, we get the following characteristic equation

$$-(\lambda + \gamma)(\lambda + \gamma_A)(D_1 + \alpha_1\beta_1V^*X^*(\lambda + \gamma_V + \delta)) = 0,$$

$\lambda = -\gamma < 0$ and $\lambda = -\gamma_A < 0$ are two eigenvalues. The others eigenvalues are solution of

$$D_1 + \alpha_1\beta_1V^*X^*(\lambda + \gamma_V + \delta) = 0.$$

This corresponds to

$$(\lambda + \gamma_A)[(\lambda + \gamma + \beta_0V^*)(\lambda + \gamma_V + \delta) - \Phi\beta_0X^*] - \alpha_1\beta_1V^*X^*(\lambda + \gamma_V + \delta) - \alpha_2\beta_2V^*X^*(\lambda + \gamma_V + \delta) = 0.$$

Developing this last equation, we obtain

$$\begin{aligned} & \lambda^3 + \lambda^2[\gamma_A + \gamma + \gamma_V + \delta + \beta_0V^*] \\ & + \lambda[\gamma_A(\gamma + \gamma_V + \delta + \beta_0V^*) + (\gamma_V + \delta)(\gamma + \beta_0V^*) - \Phi\beta_0X^* - \alpha_1\beta_1V^*X^* - \alpha_2\beta_2V^*X^*] \\ & + \gamma_A(\gamma_V + \delta)(\gamma + \beta_0V^*) - \gamma_A\Phi\beta_0X^* - \alpha_1\beta_1V^*X^*(\gamma_V + \delta) - \alpha_2\beta_2V^*X^*(\gamma_V + \delta) = 0. \end{aligned}$$

That we rewrite

$$\lambda^3 + k_1\lambda^2 + k_2\lambda + k_3 = 0, \quad (21)$$

with

$$k_1 = \gamma_A + \gamma_V + \delta + \gamma + \beta_0V^*,$$

$$k_2 = \gamma_A(\gamma_V + \delta + \gamma + \beta_0V^*) + (\gamma_V + \delta)(\gamma + \beta_0V^*) - \Phi\beta_0X^* - (\alpha_1\beta_1 + \alpha_2\beta_2)V^*X^*,$$

$$k_3 = \gamma_A((\gamma_V + \delta)(\gamma + \beta_0V^*) - \Phi\beta_0X^*) - (\alpha_1\beta_1 + \alpha_2\beta_2)(\gamma_V + \delta)V^*X^*.$$

According to the Routh-Hurwitz criteria, the endemic steady-state is locally asymptotically stable if and only if the following conditions on the coefficients of Equation (21) are satisfied

$$k_1 > 0, \quad k_3 > 0 \quad \text{and} \quad k_1k_2 > k_3.$$

It is clear that we have always

$$k_1 = \gamma_A + \gamma_V + \delta + \gamma + \beta_0V^* > 0.$$

We recall the equations that come from the steady state equation (17)

$$\gamma_A A_1^* = p_1(Y^*), \quad \gamma_A A_2^* = p_2(Y^*), \quad \gamma X^* + \beta_0 V^* X^* = \Omega, \quad \beta_0 V^* X^* = \gamma Y^* \quad \text{and} \quad (\gamma_V + \delta) V^* = \Phi Y^*.$$

Then, $k_3$ becomes

$$\begin{aligned} k_3 &= \gamma_A((\gamma_V + \delta)(\gamma + \beta_0V^*) - \Phi\beta_0X^*) - (\alpha_1\beta_1 + \alpha_2\beta_2)(\gamma_V + \delta)V^*X^*, \\ &= \gamma_A\Phi\frac{Y^*}{V^*}\left(\frac{\Omega}{X^*} - \gamma\right) - (\alpha_1\beta_1 + \alpha_2\beta_2)\Phi Y^*X^*. \end{aligned}$$

$Y^* > 0$ and $V^* > 0$, and so $0 < X^* < \Omega/\gamma$. Therefore, under the condition (11), we have $\alpha_1\beta_1 + \alpha_2\beta_2 \le 0$. Then, we directly get $k_3 > 0$. Using the above equations, we can also rewrite $k_1$ and $k_2$ as follows

$$k_1 = \gamma_A + \Phi\frac{Y^*}{V^*} + \frac{\Omega}{X^*},$$

$$k_2 = \gamma_A\left(\Phi\frac{Y^*}{V^*} + \frac{\Omega}{X^*}\right) + \Phi\frac{Y^*}{V^*}\left(\frac{\Omega}{X^*} - \gamma\right) - (\alpha_1\beta_1 + \alpha_2\beta_2)V^*X^*.$$

After some calculations, we obtain

$$\begin{aligned} k_1k_2 - k_3 = & \left(\gamma_A + \Phi\frac{Y^*}{V^*}\right)\left(\gamma_A\left(\Phi\frac{Y^*}{V^*} + \frac{\Omega}{X^*}\right) + \Phi\frac{Y^*}{V^*}\left(\frac{\Omega}{X^*} - \gamma\right)\right) \\ & + \frac{\Omega}{X^*}\left(\gamma_A\frac{\Omega}{X^*} + \Phi\frac{Y^*}{V^*}\left(\frac{\Omega}{X^*} - \gamma\right)\right) + \gamma\gamma_A\Phi\frac{Y^*}{V^*} \\ & - (\alpha_1\beta_1 + \alpha_2\beta_2)(\Omega V^* + \gamma_A V^* X^*). \end{aligned}$$

Then, Condition (11) implies that

$$k_2 > 0 \quad \text{and} \quad k_1k_2 - k_3 > 0.$$

So, the endemic equilibrium $P^* = (A_1^*, A_2^*, X^*, Y^*, V^*)$ is locally asymptotically stable.

□
## A.2 Proofs of Propositions

*Proof of Proposition 1.* The regularity of the functions used in the right-hand side of System (8) guarantees the existence and uniqueness of solutions on an interval $[0, T)$, with $T > 0$. Let us see first that the solution is nonnegative on $[0, T)$. For all $t \in [0, T)$, whenever $(A_1, A_2, X, Y, V) \in \mathbb{R}_+^5$, the derivatives satisfy

$$ \left. \frac{dA_1}{dt} \right|_{A_1=0} \geq p_1(0) \geq 0, \quad \left. \frac{dA_2}{dt} \right|_{A_2=0} \geq p_2(0) \geq 0, \quad \left. \frac{dX}{dt} \right|_{X=0} = \Omega > 0, \quad \left. \frac{dY}{dt} \right|_{Y=0} \geq 0 \quad \text{and} \quad \left. \frac{dV}{dt} \right|_{V=0} \geq 0. $$

Hence, any solution of System (8) that starts nonnegative remains nonnegative (see Theorem 3.4 in [49] and Proposition B.7 in [50] for more details).

Now, we prove that any solution of System (8) is bounded on $[0, T)$. By adding the equations of $X$ and $Y$, we get

$$ \frac{dX}{dt} + \frac{dY}{dt} \leq \Omega - \gamma_m(X + Y), \quad \text{with } \gamma_m = \min\{\gamma_X, \gamma_Y\}. $$

Then,

$$ \begin{aligned} 0 \leq X(t) + Y(t) &\leq e^{-\gamma_m t}(X(0) + Y(0)) + \frac{\Omega}{\gamma_m}(1 - e^{-\gamma_m t}), \\ &\leq X(0) + Y(0) + \frac{\Omega}{\gamma_m}. \end{aligned} $$

So,

$$ \limsup_{t \to +\infty} (X(t) + Y(t)) \leq X(0) + Y(0) + \frac{\Omega}{\gamma_m}. $$

Hence, $X$ and $Y$ are bounded on the interval $[0, T)$. Similarly, we find

$$ \limsup_{t \to +\infty} (A_1(t)) \leq A_1(0) + \frac{p_1(\bar{Y})}{\gamma_{A_1}}, $$

$$ \limsup_{t \to +\infty} (A_2(t)) \leq A_2(0) + \frac{p_2(\bar{Y})}{\gamma_{A_2}}, $$

and

$$ \limsup_{t \to +\infty} (V(t)) \leq V(0) + \frac{\Phi \bar{Y}}{\gamma_V + \delta}, \quad \text{with } \bar{Y} = \sup_{0 \leq s < T} (Y(s)). $$

Therefore the solution is bounded on $[0, T)$. This means that the solution is defined and bounded over the whole interval $[0, +\infty)$. $\square$

*Proof of Proposition 2.* The condition (H3) implies that

$$ I(Y^*) \leq \frac{\Phi \beta \Omega}{\gamma_Y \gamma_X (\gamma_V + \delta)} G \left( \frac{p_1(0)}{\gamma_{A_1}}, \frac{p_2(0)}{\gamma_{A_2}} \right) = \mathcal{R}_0 < 1, $$

for all $Y^* > 0$. Then, it is clear that there is no endemic steady-state. $\square$

*Proof of Proposition 3.* The hypothesis of Proposition 3 implies that the function $I: [0, \Omega/\gamma_Y] \to [0, +\infty)$ is decreasing with $I(0) = \mathcal{R}_0$ and $I(\Omega/\gamma_Y) = 0$. Then, the equation $I(Y^*) = 1$ has a positive solution if and only if $\mathcal{R}_0 > 1$. $\square$

## A.3 Proof of Lemma

*Proof of Lemma 4.2.* Under the condition (H1), $\gamma := \gamma_X = \gamma_Y$, we have

$$ \frac{d}{dt}(X + Y) = \Omega - \gamma(X + Y) = \gamma \left( \frac{\Omega}{\gamma} - (X + Y) \right). $$

Let $B_1 = A_1 - \frac{p_1(0)}{\gamma_{A_1}}$, $B_2 = A_2 - \frac{p_2(0)}{\gamma_{A_2}}$ and

$$ Z = \begin{cases} \frac{\Omega}{\gamma} - (X + Y), & \text{if } X + Y \leq \frac{\Omega}{\gamma}, \\ (X + Y) - \frac{\Omega}{\gamma}, & \text{if } X + Y > \frac{\Omega}{\gamma}. \end{cases} $$
It is clear that $X + Y$ cannot cross $\frac{\Omega}{\gamma}$. Then, with the new functions $B_1, B_2$ and $Z$, System (8) becomes

$$ 
\left\{
\begin{aligned}
\frac{dB_1}{dt} &= p_1(Y) - p_1(0) - \gamma_{A_1} B_1, \\
\frac{dB_2}{dt} &= p_2(Y) - p_2(0) - \gamma_{A_2} B_2, \\
\frac{dZ}{dt} &= -\gamma Z, \\
\frac{dY}{dt} &= \beta G \left(B_1 + \frac{p_1(0)}{\gamma_{A_1}}, B_2 + \frac{p_2(0)}{\gamma_{A_2}}\right) \\
&\quad \times V \left(\frac{\Omega}{\gamma} + \mu Z - Y\right) - \gamma Y, \\
\frac{dV}{dt} &= \Phi Y - (\gamma_V + \delta) V,
\end{aligned}
\right.
\qquad (22) $$

where

$$ 
\left\{
\begin{array}{ll}
\mu = -1 & \text{if } X + Y \le \frac{\Omega}{\gamma}, \\
\mu = 1 & \text{if } X + Y > \frac{\Omega}{\gamma}.
\end{array}
\right.
$$

If $\mu = 1$, then $\frac{\Omega}{\gamma} + \mu Z = \frac{\Omega}{\gamma} + X + Y - \frac{\Omega}{\gamma} = X + Y \ge 0$.

If $\mu = -1$, then $\frac{\Omega}{\gamma} + \mu Z = \frac{\Omega}{\gamma} - \left(\frac{\Omega}{\gamma} - (X + Y)\right) = X + Y \ge 0$.

So, in all cases $\frac{\Omega}{\gamma} + \mu Z \ge 0$. Moreover, the equilibrium $P_0$ becomes $0_{\mathbb{R}^5}$ for System (22) and the set $\mathbb{D} \times \mathbb{R}_+^3$ corresponds now to the non-negative orthant $\mathbb{R}_+^5$. We have already proved in Proposition 1 that $Y$ and $V$ are nonnegative, and as the function $p_i$ is nondecreasing, we also have $B_i'(t) \ge -\gamma_{A_i} B_i(t)$, for all $t \ge 0$. Hence, any solution of System (22) that starts in $\mathbb{R}_+^5$ remains in $\mathbb{R}_+^5$. Then, the subset $\mathbb{D} \times \mathbb{R}_+^3$ is invariant under System (8). $\square$

**Acknowledgements.** CPF thanks grant # 304007/2023-4, National Council for Scientific and Technological Development (CNPq). This work was supported by grants # 2019/22157-5, São Paulo Research Foundation (FAPESP) and Capes 88881.878875/2023-01. MA and CDC thank the Inria international associate team Mo-CoVec and BIO-CIVIP, STIC AmSud 23-STIC-02.

## References

[1] Einav, T. & Bloom, J. When two are better than one: Modeling the mechanisms of antibody mixtures. *PLoS Computational Biology.* **16**, 1-17 (2020)

[2] Santillán, M. On the use of the Hill functions in mathematical models of gene regulatory networks. *Mathematical Modelling Of Natural Phenomena.* **3**, 85-97 (2008)

[3] John, A. & Rathore, A. Adaptive immune responses to primary and secondary dengue virus infections. *Nature Reviews Immunology.* **19** pp. 218-230 (2019)

[4] Clapham, H., Quyen, T., Kien, D., Dorigatti, I., Simmons, C. & Ferguson, N. Modelling virus and antibody dynamics during dengue virus infection suggests a role for antibody in virus clearance, PLoS Computational Biology. (2016)

[5] Estofolete, C., Versiani, A., Dourado, F., Milhim, B., Pacca, C., Silva, G., Zini, N., Santos, B., Gandolfi, F., Mistrão, N., Garcia, P., Rocha, R., Gehrke, L., Bosch, I., Marques, R., Teixeira, M., Fonseca, F., Vasilakis, N. & Nogueira, M. Influence of previous Zika virus infection on acute dengue episode. *PLoS Neglected Tropical Diseases.* **17**, 1-25 (2023)

[6] Reich, N., Shrestha, S., King, A., Rohani, P., Lessler, J., Kalayanarooj, S., Yoon, I., Gibbons, R., Burke, D. & Cummings, D. Interactions between serotypes of dengue highlight epidemiological impact of cross-immunity. *Journal Of The Royal Society Interface.* **10**, 1-9 (2013)
[7] Tricou, V., Minh, N., Farrar, J., Tran, H. & Simmons, C. Kinetics of viremia and NS1 antigenemia are shaped by immune status and virus serotype in adults with dengue. *PLoS Neglected Tropical Diseases*. **5**, 1-10 (2011)

[8] Wen, J., Cheng, Y., Ling, R., Dai, Y., Huang, B., Huang, W., Zhang, S. & Jiang, Y. Antibody-dependent enhancement of coronavirus. *International Journal Of Infectious Diseases*. **100** pp. 483-489 (2020)

[9] Puschnik, A., Lau, L., Cromwell, E., Balmaseda, A., Zompi, S. & Harris, E. Correlation between dengue-specific neutralizing antibodies and serum avidity in primary and secondary dengue virus 3 natural infections in humans. *PLoS Neglected Tropical Diseases*. **7**, 1-8 (2013)

[10] Goncalvez, A., Engle, R., Claire, M., Purcell, R. & Lai, C. 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**, 9422-9427 (2007)

[11] Halstead, S. Dengue Antibody-Dependent Enhancement: knowns and unknowns. *Microbiology Spectrum*. **2** (2014)

[12] Parren, P. & Burton, D. The antiviral activity of antibodies in vitro and in vivo. *Advances In Immunology*. **77** pp. 195-262 (2001)

[13] Sanna, P., Ramiro-Ibañez, F. & Logu, A. Synergistic interactions of antibodies in rate of virus neutralization. *Virology*. **270**, 386-396 (2000)

[14] Anam, V., Guerrero, B., Srivastav, A., Stollenwerk, N. & Aguiar, M. Within-host models unravelling the dynamics of dengue reinfections. *Infectious Disease Modelling*. **9**, 458-473 (2024)

[15] Ansari, H. & Hesaaraki, M. A within-host dengue infection model with immune response and beddington-deangelis incidence rate. *Applied Mathematics*. **3**, 177-184 (2012)

[16] Ben-Shachar, R. & Koelle, K. Minimal within-host dengue models highlight the specific roles of the immune response in primary and secondary dengue infections. *Journal Of The Royal Society Interface*. **12** pp. 103 (2015)

[17] Gómez, M. & Yang, H. A simple mathematical model to describe antibody-dependent enhancement in heterologous secondary infection in dengue. *Mathematical Medicine And Biology: A Journal Of The IMA*. **36**, 411-438 (2018)

[18] Clapham, H., Tricou, V., Chau, N., Simmons, C. & Ferguson, N. Within-host viral dynamics of dengue serotype 1 infection. *Journal Of The Royal Society Interface*. **11** pp. 96 (2014)

[19] Gujarati, T. & Ambika, G. Virus antibody dynamics in primary and secondary dengue infections. *Journal Of Mathematical Biology*. **69**, 1148-1155 (2014)

[20] Nuraini, N., Tasman, H., Soewono, E. & Sidarto, K. A within-host dengue infection model with immune response. *Mathematical And Computer Modelling*. **49**, 1148-1155 (2009)

[21] Pawelek, K., Huynh, G., Quinlivan, M., Cullinane, A., Rong, L. & Perelson, A. Modeling within-host dynamics of influenza virus infection including immune responses. *PLoS Computational Biology*. **8** (2012)

[22] Perelson, A. Modelling viral and immune system dynamics. *Nature Reviews Immunology*. **2** pp. 28-36 (2002)

[23] Rubio, F. & Yang, H. A mathematical model to evaluate the role of memory B and T cells in heterologous secondary dengue infection. *Journal Of Theoretical Biology*. **534** pp. 110961 (2022)

[24] Tang, B., Xiao, Y., Sander, B., Kulkarni, M., Radam-Lac Research Team & Wu, J. Modelling the impact of antibody-dependent enhancement on disease severity of Zika virus and dengue virus sequential and co-infection. *Royal Society Open Science*. **7**, 191749 (2020)

[25] Wodarz, D. & Nowak, M. Mathematical models of HIV pathogenesis and treatment. *BioEssays*. **jb24j/bj1178-1187**. (2002)

[26] Xu, Z., Zhang, H., Yang, D., Wei, D., Demongeot, J. & Zeng, Q. The Mathematical Modeling of the Host-Virus Interaction in Dengue Virus Infection: A Quantitative Study. *Viruses*. **16** (2024)
[27] Adimy, M., Mancera, P., Rodrigues, D., Santos, F. & Ferreira, C. Maternal passive immunity and dengue hemorrhagic fever in infants. *Bulletin Of Mathematical Biology.* **82** pp. 24 (2020)

[28] Camargo, F., Adimy, M., Esteva, L., Métayer, C. & Ferreira, C. Modeling the relationship between antibody-dependent enhancement and disease severity in secondary dengue infection. *Bulletin Of Mathematical Biology.* **83**, 1-28 (2021)

[29] Klein, F., Mouquet, H., Dosenovic, P., Scheid, J., Scharf, L. & Nussenzweig, M. Antibodies in HIV-1 vaccine development and therapy. *Science.* **341**, 1199-1204 (2013)

[30] Scott, A., Wolchok, J. & Old, L. Antibody therapy of cancer. *Nature Reviews Cancer.* **12** pp. 278-287 (2012)

[31] Perelson, A. & Weisbuch, G. Immunology for physicists. *Reviews Of Modern Physics.* **69** pp. 1219-1268 (1997)

[32] Flamand, A., Raux, H., Gaudin, Y. & Ruigrok, R. Mechanisms of Rabies Virus Neutralization. *Virology.* **194**, 302-313 (1993)

[33] Klasse, P. Neutralization of Virus Infectivity by Antibodies: Old Problems in New Perspectives. *Advances In Biology.* **2014** pp. 1-24 (2014)

[34] Pierson, T., Fremont, D., Kuhn, R. & Diamond, M. Structural insights into the mechanisms of antibody-mediated neutralization of flavivirus infection: implications for vaccine development. *Cell Host And Microbe.* **4** pp. 229-238 (2008)

[35] N.D. Durham, A., Waltari, E., Croote, D., Zanini, F., Fouch, M., Davidson, E., Smith, O., Carabajal, E., Pak, J., Doranz, B., Robinson, M., Sanz, A., Albornoz, L., Rosso, F., Einav, S., Quake, S., McCutcheon, K. & Goo, L. Broadly neutralizing human antibodies against dengue virus identified by single B-cell transcriptomics. *ELife.* **8** pp. 1-29 (2019)

[36] Katzelnick, L., Gresh, L., Halloran, M., Mercado, J., Kuan, G., Gordon, A., Balmaseda, A. & Harris, E. Antibody-dependent enhancement of severe dengue disease in humans. *Science.* **358**, 929-932 (2017)

[37] Burton, D. Antiviral neutralizing antibodies: from in vitro to in vivo activity. *Nature Reviews Immunology.* **23**, 720-734 (2023)

[38] Kuiken, T. & Taubenberger, J. Pathology of human influenza revisited. *Vaccine.* **26**, D59-D66 (2008)

[39] Kyle, J., Beatty, P. & Harris, E. Dengue virus infects macrophages and dendritic cells in a mouse model of infection. *The Journal Of Infectious Diseases.* **195**, 1808-1817 (2007)

[40] Dörner, T. & Radbruch, A. Antibodies and B-cell memory in viral immunity. *Immunity.* **27**, 384-392 (2007)

[41] Trouplin, V., Boucherit, N., Gorvel, L., Conti, F., Mottola, G. & Ghigo, E. Bone marrow-derived macrophage production. *Journal Of Visualized Experiments.* **81** pp. 50966 (2013)

[42] Rodenhuis-Zybert, I., Wilschut, J. & Smit, J. Dengue virus life cycle: viral and host factors modulating infectivity. *Cellular And Molecular Life Sciences.* **67**, 2773-2786 (2010)

[43] Samji, T. Influenza A: understanding the viral life cycle. *The Yale Journal Of Biology And Medicine.* **82**, 153-159 (2009)

[44] Driessche, P. & Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical Biosciences.* **180**, 29-48 (2002)

[45] Kirkilionis, M. & Walcher, S. On comparison systems for ordinary differential equations. *Journal Of Mathematical Analysis And Applications.* **299**, 157-173, (2004)

[46] Guzman, M., Gubler, D., Izquierdo, A., Martinez, E. & Halstead, S. Dengue infection. *Nature Reviews Disease Primers.* **2** pp. 1-25 (2016)

[47] Ghosh, I. Within Host Dynamics of SARS-CoV-2 in Humans: Modeling Immune Responses and Antiviral Treatments. *SN Computer Science.* **2**, 482 (2021)

[48] Cambón, M. & Sánchez, Ó. Thermodynamic modelling of transcriptional control: a sensitivity analysis. *Mathematics.* **10**, 1-18 (2022)
[49] Smith, H. An Introduction to Delay Differential Equations with Applications to the Life Sciences. Springer New York, . (2011)

[50] Smith, H. & Waltman, P. The Theory of the Chemostat. Cambridge University Press, Cambridge UK, . (1995)