Abstract
The interplay between the virus, infected cells and the immune responses to SARS-CoV-2 is still under debate. Extending the basic model of viral dynamics we propose here a formal approach to describe the neutralizing versus weakly (or non-)neutralizing scenarios and compare with the possible effects of antibody-dependent enhancement (ADE). The theoretical model is consistent with data available from the literature; we show that weakly neutralizing antibodies or ADE can both give rise to either final virus clearance or disease progression, but the immuno-dynamic is different in each case. Given that a significant part of the world population is already naturally immunized or vaccinated, we also discuss the implications on secondary infections infections following vaccination or in presence of immune system dysfunctions.
1. Background
SARS-CoV-2 is a new virus from the coronavirus family, responsible for the ongoing COVID-19 pandemic. To date, there are more than 300 million cases and over five million deaths worldwide [1]. SARS-CoV-2 is the third betacoronavirus to severely infect humans appearing in the last 20 years, after SARS-CoV-1 and MERS-CoV. This motivates a growing need for efficient drugs and/or vaccines, not only for the time being but also in anticipation of a future coronavirus resurgence.
However, initial promising successes of antiviral treatments raised also the possibility of negative side-effects. On the vaccine front, an auto-immune disease (leading to temporary suspension of clinical trials) appeared during the AstraZeneca vaccine trial (Sept 9th 2020); altogether this context demonstrated the importance of understanding qualitatively and quantitatively the immune response to primary infection and also to challenges (vaccines belong to both categories). In particular, relevant mathematical models of the immune dynamics can be of interest to understand and predict the complicated behavior often observed.
We focus here on adaptive humoral immunity (the antibody-mediated immunity) and refer to future works for an extension to the cellular and/or innate immune system.
For clinical reasons and also for the understanding of those studying vaccines, antibody responses are of paramount importance. However, the neutralization capacities of antibodies is still under discussion, especially since weak or non-neutralizing antibodies can promote infection through a process called antibody-dependent enhancement (hereafter abbreviated ‘ADE’) [2, 3, 4, 5] (see appendix A.1 for details); thus the antibody neutralization capacity and ADE level are important ingredients of the model. Furthermore, such behavior may be accentuated by a challenge (secondary infection or infection following some immune system event or dysfunction), cf. also appendix A.2. Accordingly, we investigated here both primary and secondary COVID-19 infections.
To summarize, we propose a mathematical model of the immune response and virus dynamics that includes the possibility of weakly neutralizing antibodies and / or ADE. and analyze its implications. At the time of writing the final draft (January 2022) a significant part of the world population is either vaccinated or naturally immunized and the consequences of reinfection events are a major source of uncertainty concerning pandemic evolution. This naturally calls for scientific investigation.
2 Methods
2.1 Mathematical model
We present below the viral and immune response model. It is a compartmental model similar to those used to describe the epidemic propagation, see [6, 7, 8, 9] for a general introduction and [10, 11, 12, 13, 14, 15] for COVID-19 specific works.
The viral-host interaction (excluding the immune response) is called the basic model of virus dynamics. It has been extensively validated both theoretically and experimentally, see [16, eq (3.1) page 18], [17, eqns. (2.3)-(2.4) page 26] and references therein. See also [18, 19, 20] for general overviews of mathematical immunology.
The model involves several classes: that of the target cells, denoted T, the infected cells, denoted I, the free virus denoted V and the antibodies denoted A. The model is illustrated in figure 1.
Target cells T, which in our case are the epithelial cells with ACE2 receptors located, for instance in the respiratory tracts including lungs, nasal and trachea/bronchial tissues, are produced at a rate Λ and die at rate μ. The parameters Λ and μ define tissue dynamics in the absence of infection, see also section D.1. When these susceptible cells meet free virus particles V, they become infected at a rate β0. Furthermore, target cells can also become infected via ADE if virus entry is mediated by antibodies. The parameter β1 represents the rate of ADE infection route which is the result of a three-species interaction: T, A and V.
Infected (initially target) cells, denoted I, die at a rate d. Note that this death rate will often be larger than the death rate of uninfected cells because viruses cause cell damage and cell death, [17, 16]. Infected cells produce new virus particles at a rate ω, and the free virus particles which have been released from infected cells decay at a rate c called the clearance rate.
Free virions are neutralized by antibodies A, which can block virus entry into cells but also facilitate phagocytosis, at a rate b. Finally, the antibodies can be stimulated by the free virus with a production rate a while declining at a rate of s (see for instance [17, eq. (9.4) p.126]). Note that alternative proposals for the antibody dynamics exist, see e.g. André and Gandon [21] who assume that immune response, once started, grows at a constant rate while Pawlek et al. [22] design a more complex model that takes into account the macrophage activation. The complete model reads (all constants are positive):
Several hypotheses in this model need to be further documented. The first one in that all infected cells including ADE infected cells support viral replication and can produce virus. However, to date, it is still unclear whether ADE infected cells can support viral replication in vivo, [4], [5]. Here we choose not to distinguish between virus productive and non productive infected cells to keep the model simple (see however the comments in appendix E). For the same reason, we do not discriminate between neutralizing, weakly neutralizing or non-neutralizing antibodies but consider all as members of the same class, the antibodies neutralizing capacity will therefore be the average of the neutralizing power and the average is described by the parameter b; on the other hand the ADE magnitude will be monitored by parameter β1. These parameters are the most important part of the immune response and the object of our study.
3 Stability of equilibria and further considerations
We operate under the assumptions that all parameters are positive and furthermore the following two assumptions hold (see Appendix D for details):
where we define as in [16, eq. (6.2) page 53]):
Note that (7) implies in particular R0 > 1 which is a standard condition for such models. We will further denote
3.1 Stability of the equilibrium without ADE
With these definitions we can give the main theoretical properties of the model depending on the presence or not of the ADE term.
3.2 Stability of the equilibrium with ADE
We investigate now the full model having a non-null ADE term β1 > 0.
The model (1)-(4) has three equilibria:
the trivial equilibrium T = T* = Λ/μ, V = I = A = 0 which is unstable;
the immmunosuppression equilibrium, also unstable, given by :
and a third equilibrium characterized as follows:
the antibody level Af is the unique positive solution of the following second order equation in the unknown A:
the other quantities are :
The following affirmations hold true concerning this third equilibriumwhen β1 is small enough the equilibrium is stable;
when β1 is large enough the equilibrium is stable;
however there exist choices of parameters (in particular values of β1) for which this equilibrium is unstable.
Proof. The proof is presented in appendix D.4. □
3.3 Dynamical aspects
The equilibrium analysis in the previous sections does not yet tell the full story of the evolution of the system (1)-(4). Depending on the parameters, a common behavior is the following: initially A will increase as response to V being above threshold V t; the increase of A will drive both I and V to zero. Such a dynamics is stable over a long period and in practice I and V will keep small values for a time long enough to ensure virus clearance (when V is small enough, due to the random nature of the events, V will disappear).
Taking I and V to be constant equal to zero, the new evolution is:
Note that equations for I and V are missing because if the initial states are V (0) = I(0) = 0 then V (t) = I(t) for all t ≥ 0. This evolution drives T to Λ/μ and A to zero. If however during the slow decay of A a challenge is presented in the form of a virus load V > s/a a new infection will start and V and I will rise again.
In conclusion, the stable equilibrium (12)-(13) is not necessarily reached in practice. The precise dynamics depends crucially on the parameters b and β1, see main text for details.
4 Results
4.1 Theoretical results
We refer the reader to the section 3 for the rigorous statements concerning the theoretical properties of the model (1)-(5). Several situations may occur, but in summary the absence of ADE (i.e., β1 = 0) insures stable equilibrium while intermediate β1 values (neither too small not too large) may provide examples of unstable equilibria; moreover, stochastic events prevent the stable equilibrium state to be reached in practice, cf. section 3.3. The parameters b and β1 are shown to be the most important for the viral-host-antibody dynamics.
4.2 Empirical results: initial infection
Taking into account the available data from the literature and the methodology in appendix B we run a numerical procedure to fit the model parameters to reproduce at best the viral load data in figure 2 (left) and obtained the values in table 1. The numerical simulation for a primary infection corresponding to these parameters is shown in figure 3.
Left: Clinically observed typical variation of SARS-CoV-2 viral load in nasopharyngeal swab normalised using cell quantification. Data taken from [23, figure 3 page 703, patient 4]. Right: : Typical time variations for IgG. Data taken from [24, figure 2 page 1085]. Note that the antibody data is a mean over several days and corresponds to a different patient cohort.
Numerical simulation of the first infection without ADE for model (1)-(4) and parameters in table 1. A good fit for the viral load from figure 2 is obtained (data is truncated below the value 1). On the contrary the fit is not as good for A(t) because data does not correspond to the same patient (joint V (t)/A(t) data was not available).
There is a 20% fall of target cells which either become infected or naturally die. The viral load peaks around 4-5 days after symptoms onset at 1.8 × 106 copies/ml. While SARS-CoV-1 viral load, as MERS-CoV, peaked around 10 days after symptoms onset, most studies agree that SARS-CoV-2 viral load peaks sooner, around 5 days, [23],[25]. Concerning antibodies, they increase sharply until week 2 then slower until a month after infection and start to decrease within 2-3 months [26], [27]. Qualitative agreement is observed with clinically observed variations variations of viral load and partially with antibodies concentration depicted in figure 2 (see references in the figure).
Note that, although we expect agreement between V (t) and the viral load evolution in figure 2 (left) (which corresponds to a precise, real patient) the antibody data from [24, figure 2 page 1085] does not correspond to the same patient (data unavailable) but is a mean value over several days and patients (not always the same). Each individual is likely to have his own immuno-kinetic parameters: the parameters of the individual that may fit the A(t) data from figure 2 (right) are not the same as the parameters that fit the data in left side of the same figure.
The equilibrium state (10) when β1 = 0 (no ADE present) is reached after 2 years for all variables in figure 3. However, viral load and infected cells reach a minimum within several weeks post-infection before increasing and oscillating toward equilibrium state (10) (simulations not shown here). Therefore, if the virus load is very small close to the minimum, all other variables decrease to-wards 0 and the infection has vanished. The equilibrium state (10) is stable but not reached in practice as the patient is cured.
4.3 Empirical results: secondary infection, variants, vaccination
We focus on a scenario where the immuno-kinetic parameters such as the neutralizing efficacy (b) or the ADE parameter (β1) change; the causes can be multiple: a primary infection with a different variant, vaccination, or some immune evolution (aging being an example). In all cases we investigate the infection, called challenge, that takes place with a different set of b of β1 parameters than in table 1.
4.3.1 Variation of the neutralizing capacity b
When there is no ADE, decreasing the neutralizing capacity of antibodies (parameter b) leads on the one hand to a higher viral load peak but on the other hand to higher antibodies concentrations. The less neutralizing the more abundant antibodies are to compensate so that the infection is always solved. The simulations results are presented in figure 4. Infection resolution is obtained with little target cell destruction for larger values of b. On the contrary, low values of b will lead to significant increase of the antibody number and simultaneous decay of target cells, both largely pejorative for the patient.
Numerical simulation of the model (1)-(4) and parameters in table 1. Only the neutralizing capacity b is changed around the nominal value b = 0.52. Lower value of the neutralizing capacity b (solid blue line) leads to lower target cell count and higher antibody levels. See figure 6 for simulation over a longer time span.
In the cases where the viral load reaches low values the infection stops before converging to the theoretical equilibrium.
4.3.2 Presence of ADE (β1> 0)
We investigated in figure 5 the possibility of the ADE mechanism present (β1 > 0), for a range of possible parameter β1 values. We plot all variables upon challenge with the same neutralizing capacity for antibodies. A higher ADE parameter leads to more destroyed target cells, more infected cells, more viral load and more antibodies. However the antibodies concentration is restricted by an upper limit. There is a threshold effect : increasing β1 does not increase significantly the antibody population (see figure 5 and compare with theoretical insights in the proof of point 3b of proposition 2 in Appendix D). Therefore a higher β1 ADE parameter cannot be compensated by more antibodies. For example, unlike β1 = 10−8, if β1 = 10−6 the viral load directly stabilizes to its equilibrium state (13), without reaching a minimum close to 0 while oscillating (simulation not shown here). In this case, the infection wins (leading to respiratory function disruption and possibly patient death). Large values of β1 lead to significant (possibly total) destruction of target cells.
The effect of the ADE parameter β1 (the model (1)-(4)). The secondary (challenge) infection has fixed neutralizing antibodies capacity b but several possible β1 (ADE) parameter values; for all other parameters we use the nominal values given in table 1. See figure 7 for simulation over a longer time span.
5 Discussion
We investigated the immune response to SARS-CoV-2 infection and re-infection through a numerical model; the model can also take into account the possible presence of ADE, either on first infection or to a challenge (secondary or reinfection with a different phenotype, after vaccination, etc.).
As to date there is no clear evidence that ADE occurs in COVID-19 severe patients, we assume that ADE only happens upon challenge.
We started from a classic viral-host dynamics ([17],[28]) that we modified by adding parameter β1 to account for a possible ADE mechanism. In order to keep the model at its lowest complexity, we do not distinguish between ADE triggering antibodies and neutralizing antibodies.
We conducted a theoretical study of our system by computing equilibrium states and stability with and without ADE. We showed that stochastic events may also play a role and prevent the stable equilibrium state to be reached in practice; we identified the parameters b (neutralization capacity) and β1 (ADE presence) to be crucial for the dynamics of our system.
Then, we calibrated our parameters values to match reference viral load from the literature [23, 27] and obtained good results.
We then investigated a secondary infection (or infection following vaccination or other immune event) that can posess different immuno-dynamic parameters. We saw that without ADE, the possible weak antibody neutralizing capacity was systematically compensated with higher concentrations of antibody leading to viral clearance. On the other side, adding ADE was not always associated with viral clearance but possibly high target cell destruction. Simulations and equilibrium analysis showed that antibody concentration had an upper limit which prevented higher ADE to be compensated by an unlimited antibody quantity. Therefore, ADE should be taken in consideration as a serious risk in disease understanding, treatment and vaccine development and scheduling.
On the other hand, we showed that the results are sensitive to the neutralizing antibody capacity (the b parameter); note that a decrease of this parameter can occur in several situations, for instance due to immune function decay, due to the malfunctioning of the antibody immunodominance mechanism that ends up selecting too many weakly neutralizing antibodies or due to miscalibrated therapeutic interventions. Independent of the cause, such a decrease of the neutralizing capacity is susceptible to imply a substantial deterioration of the outcome.
In summary, our results seem to support a picture where ADE presence is primary correlated with important target cell destruction while loss in neutralization capacity is correlated with both higher antibody count (leading to inflammation) and larger target cell destruction.
5.1 Limitations and future work
As any other, our model contains of course several limitations. First, we considered all infected cells to support viral replication, including ADE-infected cells. Concerning SARS-CoV-2, the questions of ADE is still under debate, but for SARS-CoV-1 in vitro ADE evidence suggested abortive viral replication in ADE infected cells. Therefore, if we changed the model (1)-(4) to include this distinction, equilibrium state would change and ADE may be compensated. Similarly, we did not distinguish between former antibodies and novel antibodies secreted upon challenge. This would imply more parameters and change equilibrium levels but without inherently changing variables behavior. The antibody dynamic model can also be changed to include e.g., constant antibody production following a threshold or more specific effects [21, 22]. Regarding parameters, we did not have enough exploitable available data to train our model and fit parameters better. Finally, an unique model can hardly capture the extreme variability of COVID-19 clinical outcomes, see [29]; some studies proposed that some of the variability come from genetics, see e.g., [30] where genetic information from roughly 4,000 people from Italy and Spain was correlated to severity of COVID-19. This may lead to a variability of our model parameters in the form of random variables.
By now (January 2022), billions of people have been vaccinated, at least by one injection, and more than half a billion have been infected. Fortunately, despite the spread of the highly contagious omicron variant, it appears that morbidity and mortality are declining. This implies that, for the time being, the most disastrous consequences of the phenomena included in this model are not being observed. However, it must be emphasized that human polymorphism, viral polymorphism and highly variable environmental conditions, as well as the considerable variety of vaccination protocols, mean that there may be isolates where ADE or the other immune responses we have explored could be significant. It is therefore particularly important to monitor variations in morbidity and mortality around the world so that a rapid response can be implemented if there is any local increase. Finally, the types of vaccines used are very different. For those based on well-established technologies, we do not fore-see any consequences other than those discussed in this work, except perhaps as a function of the vaccination protocols (time lag between primary and booster injections). In contrast, the use of vaccines based on indirect antigen production (adenovirus- or synthetic RNA-based vaccines) requires specific encapsulation of the active ingredient in a variety of capsules or cassettes. These containers can, by themselves, be immunogenic. The consequence would be that after several immunisations, patients would develop a response against the vaccine, rendering it ineffective against the disease. We did not consider this outcome in our work.
The more science will shed light on the full picture of SARS-CoV-2, the more our model can input complex and precise details. In the meantime, the main take-home message is that, with parameters consistent with the available clinical data, the neutralizing capacity and ADE mechanisms may play an important immunological role into the primary and secondary infection outcomes.
Data Availability
All data is available from public sources.
A Additional motivations
A.1 Available evidence detailing the antibody response
For clinical reasons and also for the understanding of those studying vaccines, antibody responses are of paramount importance; SARS-CoV-2 specific antibodies are usually detected during the second week after illness onset (see [31, 32]) and remain active thereafter for an unknown time span (see however [33] for recent information). Antibody responses are mainly directed against the RBD-spike and nucleocapsid proteins. However, the neutralization capacities of these specific antibodies is still under discussion, especially since weak or non-neutralizing antibodies can enhance infection through a process called antibody-dependent enhancement (hereafter abbreviated ‘ADE’) [2, 3, 4, 5]. This has been recently emphasized in the set up of clinical trials (see for example [34]), in a general discussion of the prospects of vaccination [35] and in a perspective accounting for the present situation in terms of SARS-CoV-2 vaccines, therapies and immunity [36, 37]
The present academic interpretation of the ADE is that it occurs through virus-antibody immunocomplexes that facilitate virus internalization in host cells that do not express virus receptor but Fc receptors. ADE is induced when the antibody-virus stoichiometry is below the threshold for neutralization, [2, 3]. As a consequence, neutralizing antibodies may enhance infection when their concentrations fall below a key occupancy threshold, and some poorly neutralizing antibodies may strongly increase infection over a wide dose–response range. ADE has been demonstrated in vitro for many viral infections, including that triggered by SARS-CoV-1 which was reported to infect in vitro human macrophages (see [4]) and human B cell lines via an ADE pathway, (see [5]). Moreover, Qidi Wang et al. reported that a specific spike protein epitope elicited antibodies which could enhance infection via ADE, while other epitopes induced neutralizing antibodies in non-human primates. Furthermore, the authors showed that a SARS-CoV-1 inactivated vaccine could induce ADE and lung pathology in experimental rhesus monkeys [38]. In contrast, Martial Jaume and co-authors showed that vaccine candidates which mediated in vitro ADE infection could still be neutralizing and protective in vivo on rhesus macaques, [5]. Moreover, in most cases, ADE infected cells do not support viral replication, [4, 5]. Instead, ADE may trigger cell apoptosis and promote tissue inflammation and injury with the release of pro-inflammatory cytokines from infected cells, [3, 39]. As a result, whether ADE actually happens in SARS-CoV-1 infected humans and is a factor of disease severity is still a debated research subject since no in vivo human evidence has been demonstrated yet (but this statement is very time-dependent given the present intense research on SARS-CoV-2). Note however, that SARS-CoV-1 infected patients who developed a higher and earlier antibody response were associated with worse clinical outcome. An early antibody response may be weakly neutralizing compared to a later one. As a consequence, a high concentration of those antibodies could lead to ADE and enhancement of infection.
The question of ADE and the link between antibody dynamics and disease evolution is still unclear for COVID-19. J. Zhao et al. reported a strong positive correlation between disease severity and high antibody titers two weeks after illness onset. The antibody level is considered as a risk factor for severe evolution, independently from age, gender and comorbidities [40]. In another study, Wenting Tan et al also came to the same conclusion: higher titers of antiN IgM and anti-N IgG are observed for severe patients [41]. Finally, Baoqing Sun et al observed that severe patients had higher levels N-IgG than S-IgG after the symptoms onset. However, according to the authors, whether N-specific antibodies can block virus infection is still open to question [42]. The secretion of a high level of weakly neutralizing antibody supports the hypothesis of ADE for COVID-19 which can partially explain some clinical complications. In contrast, Mehul S. Suthar et al concluded that the appearance of high titer neutralizing antibody responses early after the infection was promising and may offer some degree of protection against re-infection [43]. This result seems to be confirmed in a recent study in which SARS-CoV-2 infection induced protective immunity against re-exposure in nonhuman primates. However, rhesus macaques do not develop severe clinical complications as reported in human patients, suggesting that if rhesus macaques produce neutralizing antibodies, transposition of this observation to humans is still to be investigated [44]. Finally, a recent study on a recovered cohort of COVID-19 patients showed that elderly patients had significantly higher levels of antibodies than younger patients. However, severe and critical patients were excluded from the study because they received passive antibody treatment before sample collection. As a result, the authors could not directly evaluate the effect of antibodies on virus clearance or disease progression in COVID-19 patients [45]. This suggests that if elderly patients tend to develop higher titers of antibodies, those may not be systematically associated with worse clinical evolution. What should rather be answered is whether disease severity is systematically associated with high antibody levels.
On the other hand, the vaccine community is increasingly aware of this need (see discussion on ADE in [36, 35, 37]) and studies along these lines are required.
Another motivation comes from the fact that the adaptive immune system response starts in about a week; on the other hand in many mild forms infection is resolved in around a week while on the contrary severe forms may at first start as mild and only then become severe; a simplistic view may indicate that the innate immune response is very efficient while the adaptive immune system response may be detrimental. In this case, everyone with a mild first infection (i.e., mostly dealt with by the innate immune system) will, upon re-infection, see an adaptive immune response rising faster (once the memory is in place, its response is faster than the innate immune response) and thus the detrimental effects could be visible for people previously having experienced mild forms, e.g., low age class individuals.
A.2 Evidence on re-infection
The possible unfavorable outcomes of a secondary infection (challenge) following a primary SARS-CoV-2 infection were described in various situations (see for example [46]), but an increasing body of evidence highlights the Kawasaki-like syndrome as a possible negative outcome, see [47, 48, 49, 50, 51].
An italian study [49] indicates that the immune response to SARS-Cov-2 is responsible for the appearance of a pediatric Kawasaki-like syndrome (Kawasaki-like disease or Multisystem Inflammatory Syndrome in Children MIS-C in the US). In this study, 8 to 10 children have been tested positive to IgG, IgM or both (the infection to SARS-CoV-2 preceded the development of the syndrome) and 2 only in PCR (the infection was simultaneous). SARS-CoV-2 infected children who developed the Kawasaki-like syndrome (KLS) were on average older and more severely hit than other children victims of the classical Kawasaki syndrome.
The same phenomena has been observed in the US and UK [47, 50, 51]. Academic studies begin to investivage the interplay between COVID-19 and the MIS-C [52].
The causes of the development of the Kawasaki disease are still unknown. The best accepted hypothesis is that of an abnormal immune response that occurs as a result of the infection provoked by one of several pathogenic agents for the genetically susceptible patients. The triggering pathogens have not yet been identified.
To account for the peaks of the KLS cases following an infection with SARS-CoV-2, two hypothesis may be formulated: the antibodies produced by the children can induce the initiation of an autoimmune disease and syndromes similar to the Kawasaki syndrome. The second hypothesis is an ADE-type mechanism.
B Choice of simulation parameters
Parameters’ order of magnitude were derived from literature, see [53] for μ, [23] for ω, clearance data from [28, 54, 26, 27]. To obtain the precise values, we then fitted the model to the SARS-CoV-2 clinical data available in figure 2 and obtained the values in table 1 (simulation results are shown in figure 3).
C Sensitivity with respect to parameters
We plot here a longer time evolution corresponding to figures 4 and 5. This allows to see the difference between initial dynamics and the long time equilibrium, cf. considerations in section 3.3.
D Mathematical properties of the model
We describe in an incremental way the mathematical properties of the main model (1)-(5). We take advantage of this description to illustrate the hypotheses (6) and (7). The results in sections D.1 and D.3 are known, see e.g., [55, 16, 56, 17] while those in the main text (propositions 1, 2 and their proofs in this appendix) are, to the best of our knowledge, original.
D.1 Model without a virus, nor immune response
In absence of any infection the equations for the target cells are (see [17, 16]):
Since the Jacobian matrix at equilibrium (a 1 × 1 matrix) is the constant −μ therefore the equilibrium is stable, in fact any initial data T (0) will converge to the equilibrium
D.2 Model with virus but no immune response
We employ the basic model of virus dynamics, see [16, eq (3.1) page 18] and also [17, eqns. (2.3)-(2.4) page 26] described by the equations:
The initial conditions are:
which express the fact that the initial state for T is the stable equilibrium seen in section D.1, there are initially no infected cells and the initial viral load is strictly positive.
It is natural to assume that the decay rate of infected cells is at least as large as the decay rate of healthy cells, i.e., assumption (6).
In this model, an infection is only possible if the basic reproduction ratio of the virus in the absence of immune response, defined in (8) is strictly superunitary, that is
Otherwise, that is if R0 ≤ 1, the initial viral load can only decrease. The model has two equilibria:
- trivial equilibrium: T = T* = Λ/μ, V = I = 0. The Jacobian matrix at equilibrium is
. The eigenvalues of this matrix, under condition (22), are all real but not all negative: one of them is λ1 = −μ but the product of the other two is dc− ωβ0T* ≤ 0 thus at least one is positive. Therefore, under assumption (22), this critical point is not a stable equilibrium.- the “immunosuppression” equilibrium (11)
The Jacobian matrix is 
; the characteristic polynomial P (X) = (X +d)(X +c)(X +μ +β0V is) − dc(X +μ) has the following properties: P (− ∞) < 0, P (− δ− c) = dcβ0V is > 0, thus it has a real root which is smaller than − δ − c. The product of all roots is dcβ0V is > 0 and the sum of all roots is − δ− c − μ − β0V is <− δ− c, thus the other two roots have negative real part. Therefore the equilibrium is stable.
It is important to note that the viral load V is is the viral load that the infection will cause in a completely immunodeficient individual. We expect V is to be significantly high, see in section D.3 for details.
D.3 Model: virus and immune response but no enhancement
In this section we consider the model (1)-(4) with no ADE i.e., β(A) = β0 that is β1 = 0. This model is similar to other in the literature (see for instance [17, eq. (2.9) page 29] who consider also the cytotoxic effect of the immune response on the infected cells; however they do not consider virus destruction by antibodies. In particular there virus load is constant. Another similar model is [17, eqns. (8.1)-(8.3)]. With respect to the previous section here the immune response is present. It is triggered by a threshold set at V t (see definition in (9)). It is natural to suppose that the immune response threshold is a very small value and in any case a value smaller than the immunosupression viral load V is in (11). That is we can make the hypothesis that V is > V t i.e. assumption (7) holds.
The Jacobian matrix is:
With these provisions, one can find analytically the critical points (equilibria candidates):
T = T* = Λ/μ, V = I = A = 0, which is the high dimensional analog of equilibrium (17). However, unlike in section D.1, this equilibrium is not stable any more (the determinant of the Jacobian matrix is negative when hypothesis (22) is satisfied.
the immunosuppression equilibrium (11) with A = 0; again this equilibrium is not stable any more because the condition (7) implies that the eigenvalue aV is − s is positive.
the only critical point left is (10) We prove that it is stable in the following. Note that the equilibrium value of the antibody level is positive due to condition (7).
Proof of the Proposition 1.
The equilibrium is stable when the real parts of the eigenvalues of the Jacobian matrix are negative. This is the same as saying that the roots of the polynomials P0(X) = det(X·I − J) have negative real parts (here I is the identity matrix). Such a polynomial is called stable and, if we write P0(X) = γ4X4 + γ3X3 + γ2X2 + γ1X1 + γ0 then, following the Routh-Hurwitz criterion [57], [58, p. 1076], the stability holds true when
and
Unfortunately, checking in general these conditions is very difficult because the expressions involved are highly non-linear in the original parameters of the model (a,b, c, s, etc.). We therefore need to exploit to the full extent the specific setting of our model. To this end we will make the following change of variables:
After replacing all new variables and direct computations, we obtain:
Since all parameters involved are positive we obtain that the condition (24) is satisfied. To check the remaining condition (25) we obtain
where the multi-variable polynomial Q0 is seen, after long but straightforward computations, to be equal to :
Most of the monomials in Q0 have positive coefficients, except the following ones: 
. However, in all cases we can come up with two terms that render the total sum positive. For instance the term − w2a5cβ0μ3sζ2 (term 27 of the polynomial) is negative but, when we combine it with the terms w3a6c2μ3sζ/2 (half of the third term) and 
(half of the term 79), both appearing with positive coefficients, we obtain a positive number 
.
The interested reader can check that in the same way that:
- the 36-th monomial compensate with monomials 7 and 92;
- the 38-th monomial compensate with monomials 3 and 104;
- the 61-th monomial compensate with monomials 14 and 128;
- the 73-th monomial compensate with monomials 13 and 145;
- the 75-th monomial compensate with monomials 14 and 146;
- the 87-th monomial compensate with monomials 20 and 154;
- the 89-th monomial compensate with monomials 22 and 155.
This allows to state that Q0 > 0 which concludes the proof. □
D.4 Full model: virus, immune system and ADE
Proof of the Proposition 2.
We consider the model (1)-(4) with β(A) = β0+β1A (β1 > 0). The analysis of this dynamics is more involved. The first two equilibria, having A = 0 are the complete analogues of the equilibria seen in previous sections and have no dynamical interest. Since A = 0 the parameter β1 that multiplies A has no impact and the proof of the instability of the trivial equilibrium and immunosupression equilibrium follow exactly the same arguments as before.
To find the third equilibrium, note that after immediate computations we find that the antibody level is solution of the second order equation (12). Such an equation has two solutions but exactly one is positive because the product of roots is negative; thus only a single point is an admissible equilibrium, namely the positive solution of (12) (with respect to the unknown A); setting to zero all derivatives we obtain the other values as in (13).
To prove the properties of this equilibrium we start with the point 3c of the proposition; consider the values a = s = c = b = ω = 1, μ = 1.e − 3, δ = 2, Λ = 4, β0 = 0.0011 and β1 = 0.01188; all hypotheses are satisfied and the numerical values of the equilibrium are T = 333.33, I = 1.83, V = 1, A = 0.83 while the eigenvalues are − 3.45, 0.50, 0.01 and − 0.90. Since some eigenvalues are real and positive the equilibrium is not stable for this set of parameters. This completes the proof for this point. In practice the evolution oscillates indefinitely between a state with high T value and one with very low T value.
Note that the point 3a of the conclusion is just a consequence of the continuity and the proposition 1, because both the equilibrium and the coefficients of the polynomial P (X) = det(X ×Id−J) evaluated at the equilibrium depend smoothly on β1. Since we proved that (24) and (25) are true for β1 = 0 by continuity the terms in the two conditions will remain strictly positive for β1 small enough and by the Routh-Hurwitz criterion the equilibrium will be stable.
The only point remaining to be proved is 3b. Note that when β1 → ∞ the positive root Af of the equation (12) converges to some quantity A∞, and β(Af) → ∞; moreover, we obtain from the definition of Tf that limβ1→∞ Tf = 0 and 
. Consider the Jacobian matrix:
Let us compute P (X) = det(XId − J) = det(J − XId):
Thus the polynomial P (X) = det(XId − J) can be written, to first order in β1, as P (X) = R(X) + β(A)V (X + d)(X2 + X(c + bA) + aAbV), where R(X) is a fourth order polynomial with leading term X4 and coefficients independent of β1. Note that (X + d)(X2 + X(c + bA) + aAbV) is a stable polynomial. To finish the proof we invoke Lemma 1 below for ψ = β(A)V.
□
Let Z3 = ϕ3X3 + ϕ2X2 + ϕ1X + ϕ0 be a stable polynomial of order 3 with ϕ3 > 0 and Z4 = φ4X4 + φ3X3 + φ2X2 + φ1X + φ0 a polynomial of order four with φ4 > 0. Then, for ψ large enough the polynomial Z4(X) + ψZ3(X) is stable.
Proof. Since ϕ3 > 0 using the reciprocal of the Routh-Hurwitz criterion all coefficients ϕk are strictly positive and ϕ1ϕ2 > ϕ0ϕ3. For ψ large enough this allows to check the Routh-Hurwitz criterion for the fourth order polynomial Z4(X) + ψZ3(X): the coefficients will be positive and the last remaining condition is 
, which is satisfied for ψ large enough (leading term (ϕ1ϕ2− ϕ0ϕ3)ϕ3 is positive). □
E Extended model including a latent phase
We present here a version of the main model (1)-(5) extended to take into account a latent phase of the cells. The interest of such a model is to give a finer description of all states of the attacked cells; this comes however at the price of requiring several more parameters (including the transition rate η > 0 from the latent to infected, virus-producing, cells). In practice the choice of the model depends on the outcomes of interest and available data to fit. In our case the data to fit was relatively scarce thus we kept the restricted model (1)-(5) for the numerical simulations. Denoting L the number of latent infected cells (i.e., cells already infected but not yet producing viruses) we can write this model as:
The model is illustrated in figure 8.
Simulation in figure 4 for a longer time span.
Simulation in 5 for a longer time span.
Acknowledgements
Ghozlane Yahiaoui is supported by the Engineering and Physical Sciences Research Council (EPSRC): CDT Grant Ref. EP/L015811/1.
Footnotes
↵* The order of the authors is alphabetic
Antoine.Danchin{at}normalesup.org
Oriane.Pagani-Azizi{at}espci.fr
Gabriel.Turinici{at}dauphine.fr
Ghozlane.Yahiaoui{at}maths.ox.ac.uk
Viral data fit is improved, theoretical mathematical proofs of the stability of the equilibria added.
References
- [1].↵
- [2].↵
- [3].↵
- [4].↵
- [5].↵
- [6].↵
- [7].↵
- [8].↵
- [9].↵
- [10].↵
- [11].↵
- [12].↵
- [13].↵
- [14].↵
- [15].↵
- [16].↵
- [17].↵
- [18].↵
- [19].↵
- [20].↵
- [21].↵
- [22].↵
- [23].↵
- [24].↵
- [25].↵
- [26].↵
- [27].↵
- [28].↵
- [29].↵
- [30].↵
- [31].↵
- [32].↵
- [33].↵
- [34].↵
- [35].↵
- [36].↵
- [37].↵
- [38].↵
- [39].↵
- [40].↵
- [41].↵
- [42].↵
- [43].↵
- [44].↵
- [45].↵
- [46].↵
- [47].↵
- [48].↵
- [49].↵
- [50].↵
- [51].↵
- [52].↵
- [53].↵
- [54].↵
- [55].↵
- [56].↵
- [57].↵
- [58].↵
