https://orcid.org/0000-0002-6936-5388
Figures
Abstract
The impact of variants of concern (VoC) on SARS-CoV-2 viral dynamics remains poorly understood and essentially relies on observational studies subject to various sorts of biases. In contrast, experimental models of infection constitute a powerful model to perform controlled comparisons of the viral dynamics observed with VoC and better quantify how VoC escape from the immune response. Here we used molecular and infectious viral load of 78 cynomolgus macaques to characterize in detail the effects of VoC on viral dynamics. We first developed a mathematical model that recapitulate the observed dynamics, and we found that the best model describing the data assumed a rapid antigen-dependent stimulation of the immune response leading to a rapid reduction of viral infectivity. When compared with the historical variant, all VoC except beta were associated with an escape from this immune response, and this effect was particularly sensitive for delta and omicron variant (p<10−6 for both). Interestingly, delta variant was associated with a 1.8-fold increased viral production rate (p = 0.046), while conversely omicron variant was associated with a 14-fold reduction in viral production rate (p<10−6). During a natural infection, our models predict that delta variant is associated with a higher peak viral RNA than omicron variant (7.6 log10 copies/mL 95% CI 6.8–8 for delta; 5.6 log10 copies/mL 95% CI 4.8–6.3 for omicron) while having similar peak infectious titers (3.7 log10 PFU/mL 95% CI 2.4–4.6 for delta; 2.8 log10 PFU/mL 95% CI 1.9–3.8 for omicron). These results provide a detailed picture of the effects of VoC on total and infectious viral load and may help understand some differences observed in the patterns of viral transmission of these viruses.
Author summary
SARS-CoV-2 has been characterized by the successive emergence of Variants of Concern (VoC) that have caused large epidemic rebounds. However, as VoC emerged in very different contexts of pre-existing immunity, the comparison of their intrinsic effect of viral dynamics and infectiousness remains poorly understood. Here we analysed data from 78 non-human primates infected by different VoC (Beta, Delta, Gamma and Omicron BA.1) and we used a mathematical model to quantity the impact of VoCs on viral load and infectivity. Compared with the historical variant, Omicron and Delta variants were associated with a longer and larger excretion of infectious particles, which was attributed in the model to an enhanced capability to escape the immune response. While Delta variant was associated with a larger peak viral load than the historical variant, no such effect was observed for Omicron variant. This suggests that the increased transmissibility of Omicron variant does not stem from higher viral load levels but rather from its ability to maintain high levels of infectious particles over time. Altogether, these results illustrate the importance of quantifying both viral load and infectiousness to better understand some differences observed in the patterns of viral transmission of VoCs.
Citation: Marc A, Marlin R, Donati F, Prague M, Kerioui M, Hérate C, et al. (2023) Impact of variants of concern on SARS-CoV-2 viral dynamics in non-human primates. PLoS Comput Biol 19(8): e1010721. https://doi.org/10.1371/journal.pcbi.1010721
Editor: Rustom Antia, Emory University, UNITED STATES
Received: November 8, 2022; Accepted: July 12, 2023; Published: August 9, 2023
Copyright: © 2023 Marc et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All animal data, statistical analysis code and modelling files are available on Zenodo at link https://doi.org/10.5281/zenodo.7304183.
Funding: This work was funded by the Bill & Melinda Gates Foundation through INV-017335 (JG). The NHP experiment was part of BIOVAR and PRODEVA programs funded by the ANRS-MIE project EMERGEN (ANRS0151). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The sever acute respiratory coronavirus 2 (SARS-CoV-2) is the causative agent of the Coronavirus-induced disease 2019 (COVID-19) cumulating more than 500 million cases and over 18 million death as measured by excess mortality as the end of 2022 [1,2]. Repeatedly, several variants have emerged and although most of them vanished quickly, some of them, called Variants of Concern (VoC), in particular alpha, beta, gamma, delta and omicron have caused dramatic epidemic rebounds [3–5]. These variants have acquired specific mutations enhancing their infectious capacities and escaping the immune response, leading to a dramatic loss of efficacy of monoclonal antibodies [6]. They have also caused a large drop in vaccine efficacy against disease acquisition even though until now vaccine remain largely effective against severe disease [7–9].
While several millions of individuals have been infected by these VoC, we still do not have a precise understanding on the effects of VoC on viral load. Even though some effects on larger levels of viral excretion have been reported [10–13], these studies often lack of robustness, and may be biased by many confounding factors that complicate comparisons, in particular reporting biases, heterogeneity in the incubation period and vaccination coverage.
In that context where human clinical data are difficult to interpret, the non-human primate (NHP) experimental model offers a unique opportunity to describe infection with SARS-CoV-2 in detail in a fully controlled environment. Since 2020, our group has conducted many studies to evaluate the effects of antiviral drugs or vaccines in this model [14,15], and showed its large predictive value [16]. Here, we analysed retrospectively viral load data obtained in 78 animals that were included as control arms of these studies and that were infected with different strains of SARS-CoV-2 (historical, beta, gamma, delta and omicron (BA.1)). In addition, we performed longitudinal measures of viral culture to evaluate a potential effect of VoC on viral infectivity. Using the techniques of mathematical modelling, we characterize the viral kinetics in these animals and we discuss their biological insights.
Results
Variant of concern viral kinetics
Several biomarkers were measured, both genomic RNA and subgenomic RNA were quantified at regular interval over all the study period and infectious titers at 2 times points. All macaques developed a rapid infection with genomic viral load peaking between 2- and 3-day post-infection (dpi) for the historical and beta variant, 3.5 dpi for variant delta and 4 dpi for variants gamma and omicron (BA.1). Genomic viral load was cleared at 8 dpi for the historical variant, 10 dpi for the beta variant, at 12 dpi for variants delta and omicron (BA.1) and at 14 dpi for variant gamma (Fig 1 and S1 Table). In addition to viral RNA, infectious titers were measured for 41 animals. Infectious titers were measured by Tissue Culture Infectious Dose (TCID50) from nasopharyngeal swab sampled at 2 time points per animal (day 2, 3 or 4 plus at day 5 or 7 post-infection). As we included several control animals from different studies, infected with either TCID50 or Plaque Forming Units (PFU), all TCID50 were converted to PFU assuming 1 PFU = 0.7 TCID50 [17]. All infectious titers quickly dropped to undetectable levels for the historical variant at 5 dpi, where for the other variants the infectious titers remained consistent over the course of the infection (Fig 1).
Both limit of quantification and detection are depicted as empty dots, the latter being lower. Upper limit of detection is depicted as filled squares.
Viral dynamic model
To account for the quick drop in infectious titers observed in the historical variant, (Figs 1 and S1) several models incorporating an action of an antigen-mediated immune response were tested (Fig 2). All models, except a model targeting the viral production parameter, provided an improvement of BIC compared to a target cell limited model (Table 1). We found that a model targeting the infectious ratio best described our data (Model 1 in Table 1). In the following, we discuss the parameter values of the final constructed model accounting for both an effect of the immune effector and variant specific effect on the parameters (see below). For the historical variant, we estimated the infectivity rate parameter β at 1.86×10−5 copies-1.d-1 (95% confidence interval (CI) 1×10−5–3.39×10−5) and the loss rate of infected cells δ at 1.38 d-1 (95% CI 1.22–1.55), corresponding to a half-life of 12 hours. We estimated the viral load production parameter p at 9.44×105 copies.cells-1.day-1 (95% CI 2.1×105–1.68×106). This corresponds to a within-host basic reproductive number R0 (i.e., the number of newly infected cells by one infected cell at the beginning of the infection) of 3.1 (95% CI 2–4.3) and a burst size (i.e the total number of infectious virus produced by one cell over its lifespan at the beginning of the infection) of 136 (95% CI 121−153).
The basic model is a target cell limited model without any immune response. The parameters are: β the infectivity rate, k the transfer rate between non-productive and productive infected cells, δ the loss rate of productive infected cells, p the viral production rate, μ the ratio of infectious virus, g the transfer rate between the compartments of the immune response and c the loss rate of both infectious and non-infectious virus.
VoC specific effect on viral dynamic parameters
Once an effect of the immune response was selected, a covariate search algorithm was used to find the most likely VoC associated effects (see methods) and considered the historical variant as the reference. Several variant-specific covariates were found on viral kinetics parameters that we detail below (Fig 3 and S2 Table). First, beta variant was characterized with a reduced infected cells death rate (δ) by a factor of 0.7 (95% CI 0.6−0.9) compared with the historical variant (p-value < 0.01). This led to an infected cell half-life of 17 hours and resulted in a longer period of viral load shedding as infected cells produced viruses for longer period of time. Gamma variant had an effect on the parameter θ (p-value < 0.001), the amount of immune effector F20 required to reduce by half the infectious ratio, increasing it by a factor of 9508 (95% CI 387−50 041) resulting in higher peak viral load and a longer duration of infectious virus shedding (Fig 4). Variant delta is characterized by an effect on both θ (p-value < 0.001) and the viral production parameter p (p-value < 0.05), increasing those parameters by factors 336 (95% CI 49−1191) and 1.78 (95% CI 1−3) respectively. Finally, omicron variant (BA.1) affected the parameters of the immune system θ (p-value < 0.001), the viral production rate parameter p (p-value < 0.001) and the infectious ratio μ (p-value < 0.001) modifying them by factors 229 (95% CI 27−884), 0.07 (95% CI 0.02−0.2) and 18 (95% CI 4−51) respectively (Fig 4). The model well reproduced the viral load of all animals in the individuals fits (S2 Fig). To challenge our model assumption on the immune response compartment, we performed a sensitivity analysis on the number of transitioning compartments (j) and the mean time spent in those compartments (τ). We also performed the covariate selection on all models. We found that largely similar VoC-specific covariates were selected regardless of the delay in the immune response (S3 and S4 Figs). In addition, we found that the model assuming a time spent of 3 days in the transitioning compartment yielded the best results (S3 Table).
We represent the mean value and 95% confidence interval of populations parameters for each variant. We represent only parameters having at least one variant-specific effect. Full table for population parameters is in S2 Table. The dashed black line represents the historical value.
Using simulations, we sampled parameters considering both the uncertainty in the estimation and the inter-individual variability (see methods). We represent the mean viral load of all variants and its 95% confidence interval. Dotted lines are the limits of detections.
Predicted impact of variants in a natural infection setting
The main limitation of translating these results to humans is the fact that infection in animals is done with a large inoculum dose (105−106 PFU), while human infections are presumably initiated with much lower virus dose [18]. Human experimental infections were performed with 10 TCID50 [19] in the nose, i.e., 10,000–100,000 times less virus than in the animal model. Using simulations with a lower inoculum of 1, 10 and 100 PFU (see methods), we are able to derive metrics of interest for each variant. The results obtained between 1 PFU are identical to those observed with 10 PFU, while 100 PFU predicts a more rapid kinetics (S5 Fig). We present the results with 10 PFU in the following.
The historical variant is characterized by a mean time to peak of 4.3 dpi (95% CI 3.7−4.8) and of 3.5 dpi (95% CI 3−3.9) for genomic RNA and infectious titers respectively. We found a mean peak viral load of 6.3 log10 copies/mL (95% CI 5.5−7) and of 2.1 PFU/mL (95% CI 1.2−2.9) for genomic RNA and infectious titers, respectively.
The reduced infected cell clearance rate of the beta variant resulted in a longer period of viral load shedding. The duration of the acute infection stage was consequently increased from 10.9 days (95% CI 9.5–13.1) for the historical variant to 13.4 days (95% CI 11.1–15.7) for the beta variant.
All variants except beta have shown an effect on the antigen-mediated response, greatly reducing its impact on viral kinetics. As the effect of the antigen-mediated response was reduced, the infectious ratio was increased leading to more infectious particles produced over longer periods of time. This led to the increase of the infectious titers clearance stage duration from 1.5 days for the historical variant (95% CI 0.6–1.9) to 6 days (95% CI 4.4–7.5), 3.8 days (95% CI 3.1–4.6) and 3.7 days (95% CI 2.8–4.5) for the gamma, delta and omicron variants respectively (Fig 5). This is in line with numbers of studies showing the immune escape capabilities of those variants [20–22].
An effect increasing the viral production parameter (p), as observed for the delta variant, results in largely higher peak viral load of 7.6 log10 copies/mL (95% CI 6.8–8.2) and peak infectious titers of 3.7 PFU/mL (95% CI 2.4–4.6). Conversely, an effect reducing the viral production parameter, as observed for the omicron variant, results in lower peak viral load compared to the historical variant of 5.6 log10 copies/mL (4.8–6.3) but very similar peak infectious titers at 2.8 PFU/mL (95% CI 1.9–3.8). This is due to an effect of omicron on the infectious ratio, increasing the proportion of infectious virus produced.
We represent the mean and 95% confidence interval for each variant. The dashed black line represents the historical mean value.
Discussion
Here, we used mathematical models to characterize in detail the SARS-CoV-2 viral dynamics with the main variants of concern, using data obtained in an experimental model of non-human primates. We evaluated the impact of an antigen-mediated immune response on the viral dynamics and found that an effect reducing the infectious ratio best described our data. Gamma, delta and omicron (BA.1) variants showed a significant ability to escape this response, leading to a higher ratio of infectious virus than the historical variant, especially at later timepoints. Interestingly, delta variant was associated with an increased viral production rate and therefore higher viral loads, whereas omicron variant was associated with a lower viral production rate. Altogether, our model predicts that omicron infections are associated with lower peak viral RNA and reduced duration of viral RNA clearance compared to delta variant, while their kinetics of infectious titers are similar. Accordingly, it suggests that the increased transmissibility of omicron variant is not caused by increased viral load, but rather to its ability to maintain high levels of infectious particles, and may suggest longer duration of infectiousness for both delta and omicron variants. They illustrate that the quantification of infectious titers over time is crucial to inform further public health policies and adjust the isolation period accordingly [23].
Our study has important limitations that need to be acknowledged. First, in this experimental model the inoculated dose is extremely high compared to a natural infection, which is typically initiated with 1–10 infectious particles (19). This is inherent limitation of NHP studies leads to a rapid saturation of target cells making it difficult to accurately estimate the early phase of infection. Here, however, our prediction obtained assuming a lower inoculum is in line with what was observed in humans, with a peak viral load around 5 days post infection [24,25]. In the future, studies evaluating viral dynamics with lower inoculum could be helpful to tease out a potential impact of the inoculum size on viral dynamics. Second, as our model selection process was performed using all data available, it assumes that the same model applies to all variants, which may not be true. Nonetheless, data limitation did not allow to perform model selection on each variant separately. Third, we developed an extension of the target cell limited model considering the effect of an antigen-mediated immune response targeting the infectious ratio parameter μ. This effect on the infectious ratio could be attributable to an effect of the innate response, in particular the interferon (IFN) response, that can confer antiviral, antiproliferative and immunomodulatory functions [26]. Our model also identified that the effects of the immune response was variant-dependent (captured by the covariates on the parameter θ in the model) consistent with the observations that SARS-CoV-2 variants exhibits mutations in nucleocapsid, membrane and non-structural proteins that antagonize IFN signalling in cells [27,28]. Unfortunately, this hypothesis remains speculative as IFN was not measured in this study, and the cytokines that were measured (e.g. MCP1, IL15 and IL1RA) did not correlate with the viral kinetics (S7 Fig). Of note, these cytokines peaked at day 2, i.e., earlier than predicted by our model which assumes a peak of the immune response by day 6 (S6 Fig). This could be due to the fact that these cytokines were measured in plasma, and not in the site of infection, but we nevertheless tested the sensitivity of this prediction to our hypotheses. For that purpose, we conducted a sensitivity analysis on the kinetics of the immune response by assuming different number of compartments (j) and mean time spent in those compartments (τ), allowing to modulate the mean delay of immune response and its distribution. Interestingly, in all simulations, the model predicted a peak of the immune response later than day 4. Using the final model, we nevertheless show that the impact of the immune response in reducing the infectious ratio is already significant at day 3 post infection (S6B Fig).
Regarding virus infectivity, several aspects need to be discussed regarding both their measurements and their interpretation. First, infectious titers are only a measure of in vitro infectivity, and to what extent they translate into infectiousness is unknown. Second, recent results obtained with other experimental cell line suggest that Vero E6 cells may underestimate viral production of Omicron-variants as represented by the parameter θ [29]. Related to that question, one should acknowledge that the amount of infectious virus remains limited, both in terms of number of points and in the capability to measure it precisely. Although the data available here allowed us to estimate all parameters with good precision, including the changes in infectivity over time, some aspects of the model (such as the effect of the immune system on infectivity, noted θ) remain phenomenological. Finally, in a context where more than half of the world population has received at least one dose of COVID-19 vaccine [30], and probably even more have been infected by one or different variants, our results will need to be complemented by studies to evaluate the impact of vaccination, previous infection or both on viral dynamics.
Materials and methods
Experimental procedure
Data comes from studies performed on cynomolgus macaques to evaluate the viral dynamics of SARS-CoV-2 variants. Our study includes 78 cynomolgus macaques (Macaca fascicularis) coming from control arms of several studies and have received no pharmacological interventions besides placebo. All animals were infected with doses ranging from 7×104 to 106 PFU of different SARS-CoV-2 strains. Animals are infected via both nasopharyngeal and intratracheal route with 10% of the initial volume administered in the nose and 90% in the trachea. The study is composed of 5 groups, each infected with a different SARS-CoV-2 strains: 44 Historical (hCoV-19/France/lDF0372/2020 strain; GISAID EpiCoV platform under accession number EPI_ISL_406596), 9 Bêta (B.1.351—hCoV-19/USA/MD-HP01542/2021, BEI NR-55283), 5 Gamma (P.1 - hCoV-19/Japan/TY7-503/2021, BEI NR-54984), 11 Delta (B.1.617.2—hCoV-19/USA/MD-HP05647/2021, BEI NR-55674) and 9 Omicron (B.1.1.529 –hCoV-19/USA/MD-HP20874/2021, BEI NR-56462). For each group both genomic RNA and subgenomic RNA swab samples were quantified using real time PCR in both the nasopharynx and in the trachea. For 41 animals (13 Historical, 3 Beta, 5 Gamma, 11 Delta and 7 Omicron (BA.1)) infectious titers were measured at 2 time points, early (2, 3 or 4 days post infection) and late (5 or 7 days post infection) using Tissue Culture Infectious Dose (TCID50) from nasopharyngeal swab samples [16]. As we included animals from different studies that were inoculated with different methods (PFU or TCID50), we normalized all measures of infectious titers by converting all TCID50 measurements to Plaque Forming Units (PFU) using the Formula 1 PFU = 0.7 TCID50 [17]. As no infectious titers were measured in the trachea samples, we focused the main analysis on the nasopharyngeal compartment. The results mainly focus on the genomic viral load as the subgenomic is a directly proportional to the latter.
Basic viral dynamic model
We used a previously described model of SARS-COV-2 viral dynamics to reconstruct the nasopharyngeal viral load of infected animals. In this model, target cells (T) become infected cells (I1) at a rate β. Infected cells transition into productive infected cells (I2) at a rate k and produce infectious virus (VI) at a rate pμ and non-infectious virus (VNI) at a rate p(1−μ). Productive infected cells are cleared at a rate δ and both infectious and non-infectious virus are cleared at a rate c. The basic within-host reproductive number, representing the number of newly infected cells by one infected cell, is and the burst-size, representing the number of infectious virus produced by one infected cell over its lifespan, is
. The model is described with the following set of ordinary differential equations:
(1)
(2)
(3)
(4)
(5)
Assumption on parameter values
Some parameters of the model were fixed to ensure identifiability. The transfer rate parameter between infected cells and productive infected cells was fixed to k = 4 day-1 (corresponding to a mean duration of the eclipse phase, i.e. the time for infected cells to start producing viruses, of hours) [31]. The viral clearance c was set to 10 day-1 based on previous work [14,16,24]. As only the product pT0 is identifiable, we choose to fix the initial number of target cell to T0 = 12 500 cells following the same assumptions as in [16]. Animals are infected using both an intranasal and intratracheal route using 90% of the initial inoculum in the trachea and 10% in the nose [13]. We take this into account by adding a 0.1 factor in our initial conditions. We introduced a parameter h representing the proportion of the inoculum that arrive on the site of infection. Because this parameter cannot be reliably estimated or experimentally measured, a profile likelihood approach was used, testing values ranging from 10% to 90%. Therefore, our initial conditions were set to:
(6)
(7)
(8)
(9)
(10)
Where VI(0)i is the number of infectious virus of subject i (obtained assuming that 1 PFU = 1 infectious particle), VNI(0)i is the number of non-infectious virus (obtained by the difference between the total number of RNA copies and the number of infectious titers in the inoculum) and hi the proportion of the inoculum actively initiating the infection.
Models incorporating antigen-mediated immune response
To account for the quick drop in infectious titers observed for the historical variant (Figs 1 and S1), we introduced a model incorporating an action of an antigen-mediated immune response. To allow some variability in this response we assumed a certain production phase before the immune response takes effect. We introduced this delayed effector compartment using the Linear Chain Trick (LCT) [32] by modelling a successive number of transitioning compartments. This assumes an Erlang distribution with parameters j and τ representing the number of transitioning compartments and the mean time spent in those compartments respectively. We fixed those parameters to j = 20 compartments and τ =3 days to account for the setup of the immune response. We then performed a sensitivity analysis on those parameters varying both the number of compartments (from 5 to 30) and the mean time spent in those compartments (from 1 day to 6).
The equations for the transfer compartments are written as follows:
(11)
(12)
(13)
In the following only the compartment F20 will serve as the effector for the action of the immune system. The transfer rate parameter g is then written as and fixed to 6.67 d-1 and the loss rate of the final effector dF is fixed to 0.4 d-1 [24]. All immune response compartments were set to 0 at t = 0. Several modes of action of the response system were tested:
Model 1: Immune effector decreases the infectious ratio μ
In this model, the immune effector directly decreases the infectious ratio parameter μ using an Emax function type expression:
(14)
(15)
With θ being the amount of immune effector F20 needed to reduce by half the infectious ratio.
Model 2: Immune effector increases infected productive cells death rate δ
The death rate of infected cells is increased in proportion to the amount of immune effector F20.
Where φ is the strength of the immune system.
Model 3: Immune effector reduces the infectivity rate β
In this model, the immune effector blocks virus entry in the cells by reducing the infectivity parameter β.
Model 4: Immune effector reduces the production rate p
In the same way as model 1, the viral load production parameter is reduced by the immune effector with an Emax type function:
(19)
(20)
All models were compared based on the Bayesian Information Criterion (BIC). We selected the model that yielded the lowest BIC and the best individual fits.
Model 5: Immune effector induce refraction to infection.
The immune response renders target cells refractory to infection [24] as follows:
(21)
Models were estimated using all data (historical and variants) and compared based on BIC, precision of the estimation and goodness of fit. A covariate search algorithm was then used on the selected model (see below) to identify VoC-specific effects.
Statistical model
Parameter estimation was performed using non-linear mixed effect modelling. The statistical models describing the genomic RNA, subgenomic RNA and the infectious titers are:
(22)
(23)
(24)
Where the superscript 1, 2 and 3 refers to the genomic RNA, subgenomic RNA and infectious titers, respectively. We denote yij is the jth observation of subject i at time tij, with i ∈ 1, …, N and j ∈ 1, …, ni with N the number of subject and ni the number of observations for subject i. The function describing the total viral load kinetics V(tij, Ψi) predicted by the model at time tij defined as: VI(tij, Ψi)+VNI(tij, Ψi) predicted by the model at time tij. The vector of individual parameters of subject i is noted Ψi and ϵij is the additive residual Gaussian error of constant standard deviation σ. The vector of individual parameters depends on a fixed effects vector and on an individual random effects vector, which follows a normal centered distribution with diagonal variance-covariance matrix Ω. All parameters follow a log-normal distribution to ensure positivity except both parameters μ and h which follows logit-normal distribution and are bounded between 0 and 1. We assumed random effect on all parameters and removed them using backward procedure, if they were < 0.1 or their RSE > 50%. All biomarkers (i.e. genomic RNA, subgenomic RNA and infectious titers) were fitted simultaneously.
Selection of variant-specific effect on the viral dynamic parameters
Using the best model selected at the previous step, we sought to identify VoC-specific effect on the parameters of the model (β, δ, p, μ and θ). We first performed a backward selection of the random effects removing non-significant ones (i.e. relative standard error > 50%) if the BIC wasn’t degraded by more than 2 points. We then used the Conditional Sampling use for Stepwise Approach on Correlation tests (COSSAC) to identify variant specific effect [33]. Then a backward procedure was used to remove any non-significant covariate effect with a Wald test (i.e. the covariate was removed if its coefficient effect relative standard error was > 50%). This procedure was repeated until all nonsignificant covariate effects had been eliminated. Additionally, we performed a sensitivity analysis on our best structural model. We tested for several delays in the establishment of the antigen-mediated effector (from 1 to 6 days) and on the number of transitions compartments (from 5 to 30) and then performed the covariate search on all model combinations.
Simulation of natural human infection
Finally, we used our final model to assess the impact of variants of concern on viral load and viral infectivity in a natural infection setting. We used a starting inoculum of 10 infectious virus, as described in an experimental challenge conducted in England [19] to simulate a human infection. The initial conditions are then written as:
We provided confidence interval on the mean predicted viral load, considering both the uncertainty in the estimation and the inter-individual variability. We first sampled M = 100 population parameters in their estimation distribution and then, for each variant, sampled N = 30 individual parameters from each sets of population parameters (leading to 3000 individual parameters per variant). We calculated the predicted viral load of all individuals and derived the mean viral load over the simulated individuals at all times with its 95% inter quantile range. Additionally, we provided the distribution of several viral dynamic metrics, namely:
- the area under viral load curve,
- the peak and time to peak viral load
- the duration of the clearance stage, calculated as the time interval between the peak viral load and the time to undetectable viral load
- the duration of the acute phase, calculated as the time between the first and the last detectable viral load [34].
Additional simulations were done with different inoculum (1 and 100) to assess sensibility to initial conditions.
Supporting information
S1 Fig. Relationship between genomic RNA and infectious titers.
We represent the longitudinal values of genomic RNA for each individual and if the associated PFU sampe is detectable or not.
https://doi.org/10.1371/journal.pcbi.1010721.s001
(PDF)
S2 Fig. Individual fit of genomic RNA, subgenomic RNA and infectious titers in all animals.
Undetectable values are represented as empty dots. Values above the upper limit of quantification are represented as squares.
https://doi.org/10.1371/journal.pcbi.1010721.s002
(PDF)
S3 Fig. Sensitivity analysis on the covariate selection algorithm.
We performed a sensitivity analysis on our best model. The model IDs are represented on top, as described in S3 Table. The scale represents the magnitude of the covariate effect rescaled for each row with 0 being the minimum value and 1 the maximum. Empty tiles indicate that no covariates were selected for this variant-parameter relationship.
https://doi.org/10.1371/journal.pcbi.1010721.s003
(PDF)
S4 Fig. Consistency of the covariate selection algorithm.
We represent the number of times a covariate was found on a variant-parameter relationship across all 24 models. Empty tiles indicate that no covariates were found for this variant-parameter relationship.
https://doi.org/10.1371/journal.pcbi.1010721.s004
(PDF)
S5 Fig. Impact of different inoculum on viral dynamic simulations.
Using simulations, we sampled parameters considering both the uncertainty in the estimation and the inter-individual variability (see methods). Only the mean viral load was shown for clarity.
https://doi.org/10.1371/journal.pcbi.1010721.s005
(PDF)
S6 Fig. Dynamics of the immune response and its effect.
A) Median trajectory of the last compartment of our immune response B) Median trajectory of the infectious ratio parameter μ over time. We used the population parameters of our best model to simulate the median trajectory of each variant.
https://doi.org/10.1371/journal.pcbi.1010721.s006
(PDF)
S7 Fig. Dynamics of cytokines and correlation with viral load metrics.
A) Median concentration of measured cytokines. B) Correlation between AUC of viral load predicted by our model and the cytokine AUC. C) Correlation between peak viral load predicted by our model and peak cytokine concentration.
https://doi.org/10.1371/journal.pcbi.1010721.s007
(PDF)
S1 Table. Characteristics of the 78 animals analysed.
Descriptive statistics of the animals calculated on the raw data.
https://doi.org/10.1371/journal.pcbi.1010721.s008
(DOCX)
S2 Table. Estimates of the population parameter and covariate effects for the best model.
**The standard error for the R0 parameters were calculated using the delta method.
https://doi.org/10.1371/journal.pcbi.1010721.s009
(DOCX)
S3 Table. Sensitivity analysis on the delayed immune response.
Using the best structural model (i.e. Model 1 including an effect on the infectious ratio) we tested several delays for the immune response to take place and performed the covariate search algorithm on all models.
https://doi.org/10.1371/journal.pcbi.1010721.s010
(DOCX)
Acknowledgments
We would like to thank everyone in the CEA and at Pasteur Institute that have helped for data collection. We thank Alan Perelson for helpful discussions.
References
- 1.
WHO Coronavirus (COVID-19) Dashboard | WHO Coronavirus (COVID-19) Dashboard With Vaccination Data. [cited 25 Apr 2022]. Available: https://covid19.who.int/data
- 2. Wang H, Paulson KR, Pease SA, Watson S, Comfort H, Zheng P, et al. Estimating excess mortality due to the COVID-19 pandemic: a systematic analysis of COVID-19-related mortality, 2020–21. The Lancet. 2022;399: 1513–1536. pmid:35279232
- 3. Campbell F, Archer B, Laurenson-Schafer H, Jinnai Y, Konings F, Batra N, et al. Increased transmissibility and global spread of SARS-CoV-2 variants of concern as at June 2021. Eurosurveillance. 2021;26: 2100509. pmid:34142653
- 4.
World Health Organization. COVID-19 weekly epidemiological update, edition 42, 1 June 2021. World Health Organization; 2021 Jun. Available: https://apps.who.int/iris/handle/10665/341622
- 5.
World Health Organization. COVID-19 weekly epidemiological update, edition 110, 21 September 2022. World Health Organization; 2022 Sep. Available: https://apps.who.int/iris/handle/10665/363125
- 6. Choudhary MC, Chew KW, Deo R, Flynn JP, Regan J, Crain CR, et al. Emergence of SARS-CoV-2 Resistance with Monoclonal Antibody Therapy. medRxiv. 2021; 2021.09.03.21263105. pmid:34545376
- 7.
Effect of Covid-19 Vaccination on Transmission of Alpha and Delta Variants | NEJM. [cited 27 Sep 2022]. Available: https://www.nejm.org/doi/full/10.1056/nejmoa2116597
- 8.
Virological characteristics of SARS-CoV-2 vaccine breakthrough infections in health care workers | NCRC. In: 2019 Novel Coronavirus Research Compendium (NCRC) [Internet]. 21 Aug 2021 [cited 27 Sep 2022]. Available: https://ncrc.jhsph.edu/research/virological-characteristics-of-sars-cov-2-vaccine-breakthrough-infections-in-health-care-workers/
- 9. Andrews N, Stowe J, Kirsebom F, Toffa S, Rickeard T, Gallagher E, et al. Covid-19 Vaccine Effectiveness against the Omicron (B.1.1.529) Variant. New England Journal of Medicine. 2022;386: 1532–1546. pmid:35249272
- 10. Elie B, Roquebert B, Sofonea MT, Trombert-Paolantoni S, Foulongne V, Guedj J, et al. Variant-specific SARS-CoV-2 within-host kinetics. Journal of Medical Virology. n/a. pmid:35373851
- 11. Blanquart F, Abad C, Ambroise J, Bernard M, Cosentino G, Giannoli J-M, et al. Characterisation of vaccine breakthrough infections of SARS-CoV-2 Delta and Alpha variants and within-host viral load dynamics in the community, France, June to July 2021. Eurosurveillance. 2021;26: 2100824. pmid:34533119
- 12. Naveca FG, Nascimento V, de Souza VC, Corado A de L, Nascimento F, Silva G, et al. COVID-19 in Amazonas, Brazil, was driven by the persistence of endemic lineages and P.1 emergence. Nat Med. 2021; 1–9. pmid:34035535
- 13. Li B, Deng A, Li K, Hu Y, Li Z, Shi Y, et al. Viral infection and transmission in a large, well-traced outbreak caused by the SARS-CoV-2 Delta variant. Nat Commun. 2022;13: 460. pmid:35075154
- 14. Maisonnasse P, Guedj J, Contreras V, Behillil S, Solas C, Marlin R, et al. Hydroxychloroquine use against SARS-CoV-2 infection in non-human primates. Nature. 2020 [cited 5 Aug 2020]. pmid:32698191
- 15. Maisonnasse P, Aldon Y, Marc A, Marlin R, Dereuddre-Bosquet N, Kuzmina NA, et al. COVA1-18 neutralizing antibody protects against SARS-CoV-2 in three preclinical models. Nat Commun. 2021;12: 6097. pmid:34671037
- 16. Antonio G, Pauline M, Flora D, Mélanie A, Sylvie B, Vanessa C, et al. SARS-CoV-2 viral dynamics in non-human primates. PLOS Computational Biology. 17: e1008785. pmid:33730053
- 17.
Davis BD, Dulbecco R, Eisen HN, Ginsberg HS, Wood WBBD Davis. Nature of viruses. New York: Harper and Row. 1972.
- 18. Tang S, Mao Y, Jones RM, Tan Q, Ji JS, Li N, et al. Aerosol transmission of SARS-CoV-2? Evidence, prevention and control. Environ Int. 2020;144: 106039. pmid:32822927
- 19. Killingley B, Mann AJ, Kalinova M, Boyers A, Goonawardane N, Zhou J, et al. Safety, tolerability and viral kinetics during SARS-CoV-2 human challenge in young adults. Nat Med. 2022; 1–11. pmid:35361992
- 20. Chakraborty C, Sharma AR, Bhattacharya M, Lee S-S. A Detailed Overview of Immune Escape, Antibody Escape, Partial Vaccine Escape of SARS-CoV-2 and Their Emerging Variants With Escape Mutations. Frontiers in Immunology. 2022;13. Available: https://www.frontiersin.org/articles/10.3389/fimmu.2022.801522 pmid:35222380
- 21. Mlcochova P, Kemp SA, Dhar MS, Papa G, Meng B, Ferreira IATM, et al. SARS-CoV-2 B.1.617.2 Delta variant replication and immune evasion. Nature. 2021;599: 114–119. pmid:34488225
- 22. Willett BJ, Grove J, MacLean OA, Wilkie C, De Lorenzo G, Furnon W, et al. SARS-CoV-2 Omicron is an immune escape variant with an altered cell entry pathway. Nat Microbiol. 2022;7: 1161–1179. pmid:35798890
- 23. Qiu X, Nergiz AI, Maraolo AE, Bogoch II, Low N, Cevik M. The role of asymptomatic and pre-symptomatic infection in SARS-CoV-2 transmission—a living systematic review. Clinical Microbiology and Infection. 2021;27: 511–519. pmid:33484843
- 24. Néant N, Lingas G, Hingrat QL, Ghosn J, Engelmann I, Lepiller Q, et al. Modeling SARS-CoV-2 viral kinetics and association with mortality in hospitalized patients from the French COVID cohort. PNAS. 2021;118. pmid:33536313
- 25. Gonçalves A, Bertrand J, Ke R, Comets E, de Lamballerie X, Malvy D, et al. Timing of antiviral treatment initiation is critical to reduce SARS-CoV-2 viral load. CPT Pharmacometrics Syst Pharmacol. 2020; psp4.12543. pmid:32558354
- 26. Thiel V, Weber F. Interferon and cytokine responses to SARS-coronavirus infection. Cytokine & Growth Factor Reviews. 2008;19: 121–132. pmid:18321765
- 27. Shalamova L, Felgenhauer U, Wilhelm J, Schaubmar AR, Büttner K, Schoen A, et al. Omicron variant of SARS-CoV-2 exhibits an increased resilience to the antiviral type I interferon response. PNAS Nexus. 2022;1: pgac067. pmid:36713328
- 28. Guo K, Barrett BS, Morrison JH, Mickens KL, Vladar EK, Hasenkrug KJ, et al. Interferon resistance of emerging SARS-CoV-2 variants. Proceedings of the National Academy of Sciences. 2022;119: e2203760119. pmid:35867811
- 29. Planas D, Bruel T, Staropoli I, Guivel-Benhassine F, Porrot F, Maes P, et al. Resistance of Omicron subvariants BA.2.75.2, BA.4.6, and BQ.1.1 to neutralizing antibodies. Nat Commun. 2023;14: 824. pmid:36788246
- 30. Ritchie H, Mathieu E, Rodés-Guirao L, Appel C, Giattino C, Ortiz-Ospina E, et al. Coronavirus Pandemic (COVID-19). Our World in Data. 2020 [cited 26 Sep 2022]. Available: https://ourworldindata.org/covid-vaccinations
- 31. Ke R, Zitzmann C, Ho DD, Ribeiro RM, Perelson AS. In vivo kinetics of SARS-CoV-2 infection and its relationship with a person’s infectiousness. Proceedings of the National Academy of Sciences. 2021;118: e2111477118. pmid:34857628
- 32. Hurtado PJ, Richards C. Building mean field ODE models using the generalized linear chain trick & Markov chain theory. Journal of Biological Dynamics. 2021;15: S248–S272. pmid:33847236
- 33. Ayral G, Si Abdallah J-F, Magnard C, Chauvin J. A novel method based on unbiased correlations tests for covariate selection in nonlinear mixed effects models: The COSSAC approach. CPT: Pharmacometrics & Systems Pharmacology. 2021;10: 318–329. pmid:33755345
- 34. Kissler SM, Fauver JR, Mack C, Olesen SW, Tai C, Shiue KY, et al. Viral dynamics of acute SARS-CoV-2 infection and applications to diagnostic and public health strategies. PLOS Biology. 2021;19: e3001333. pmid:34252080
- 35. Comets E, Lavenu A, Lavielle M. Parameter Estimation in Nonlinear Mixed Effect Models Using saemix, an R Implementation of the SAEM Algorithm. Journal of Statistical Software. 2017;80: 1–41.
- 36.
Monolix. In: Lixoft [Internet]. [cited 24 May 2020]. Available: http://lixoft.com/products/monolix/