Viral dynamic modeling of SARS-CoV-2 in nonhuman primates

Non-human primates infected with SARS-CoV-2 exhibit mild clinical signs. Here we used a mathematical model to characterize in detail the viral dynamics in 31 cynomolgus macaques infected with 10 6 pfu of SARS-CoV-2 for which nasopharyngeal and tracheal viral load were frequently assessed. We identied that infected cells had a large daily viral production (>10 4 virus) and a within-host reproductive basic number of 6 and 4 in nasopharyngeal and tracheal compartment, respectively. After peak viral load, infected cells were rapidly cleared with a half-life of 9 hours, with no signicant association between cytokine elevation and clearance. Translating our model to the context of human-to-human infection, human mild infection may be characterized by a peak occurring 4 days after infection, a viral shedding of ~11 days and a generation time of 4 days. These results improve the understanding of SARS-CoV-2 viral replication and better understand the infection to SARS-CoV-2 in humans.


Introduction
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) which originated in Wuhan, China, at the end of December 2019, has spread rapidly around the world, resulting at the end of July 2020 in more than 500,000 deaths [1]. Fortunately, the majority of infections do not lead to hospitalizations, and about 80% of subjects infected with SARS-CoV-2 will experience asymptomatic or pauci-symptomatic infection characterized by speci c (anosmia) or general symptoms (fever, fatigue) [1][2][3] . However, in 20% of cases evolution toward more severe symptoms can take place leading to pneumonia or acute respiratory distress syndrome 4 . In other acute or chronic viral diseases (HIV, HCV, in uenza), the characterization of viral load kinetics has played an important role to understand the pathogenesis of the virus and design better antiviral drugs [5][6][7] . In the case of SARS-CoV-2, viral kinetics remain poorly characterized and its association with disease evolution is controversial. This is due to the fact that many studies rely on large transversal analyses with few patients having serial data points or, in contrary, on detailed small series of patients, often with mild diseases [8][9][10] . In that perspective, the analysis of data generated in non-human primates is a unique opportunity to characterize in detail the viral dynamics during natural infection, and to study the effects of antiviral therapy [11][12][13] .
Here, we used data generated on nonhuman primates (NHP) infected by SARS-CoV-2. In this study, 31 cynomolgus macaques were infected with 10 6 pfu of a SARS-CoV-2 isolate and treated with hydroxychloroquine (HCQ) ± azithromycine (AZM) in either pre-or post-exposure prophylaxis 14 . Because our analysis did not nd any support for any antiviral effect of HCQ against SARS-COV-2, we here analyzed together all animals and we neglected treatment with HCQ ± AZM as a rst approximation. We develop a mathematical model of SARS-CoV-2 infection in NHP to provide estimation of key parameters driving viral dynamics.

SARS-CoV-2 viral kinetics
A total of 31 animals were infected using a combined intra-nasal and intra-tracheal infection with 10 6 pfu of a primary SARS-CoV-2 isolate (BetaCoV/France/IDF/0372/2020). Animals developed a rapid infection, with viral loads peaking 2 days post infection (dpi) in both nasopharyngeal and tracheal compartments.After the peak, both nasopharyngeal and tracheal viral loads rapidly declined exponentially, with a median rate of 0.8 log 10 copies/mL every day (Fig. 1).The slope in viral decline was more rapid in the trachea than in the nasopharynx, leading to a rst measurement below the limit of quanti cation (LOQ=8514 copies/mL) 7 and 9 dpi, respectively. Overall, the area under the viral load curve (AUC) was larger in the nasopharynx than in the trachea (45 vs 38 log 10 copies.day/mL, p<10 -4 ).
Importantly, there was no evidence of an antiviral effect in animals receiving various doses of HCQ compared to those receiving vehicles, even after adjustement on HCQ exposure 14 . Thus, the effect of treatment was neglected as a rst approximation. We later challenged this hypothesis but found no effect of HCQ in reducing viral production (see supplementary information le 1).

Viral kinetic model
Nasopharynx and trachea were considered as two distinct compartments of the upper respiratory tract (URT) where each one is described by a target cell limited model 11,15,16 . The model given by ordinary differential equations (1) to (5) describes theinterplay between susceptible cells (T), non-productive infected cell (I 1 ), productive cells (I 2 ) and infectious and non-infectious viruses (V I and V NI , respectively).
The subscript X denotes the compartment of interest either nasopharynx or trachea. A schematic representation of the model is given in Fig. 2.
Through respiration and movements of ciliary cells, viruses can migrate from one compartment to the other. Thus, we tested the possibility for viruses to exchange from nasopharyngeal to tracheal compartment and vice versa by linking both with a bidirectional rate constant g.However, due to the sparseness of the data, the parameter g could not be precisely estimated (CI 95% included 0) and was set to 0. Then, using a backward selection procedure we found that the infectivity rate β and the viral production p were different between nasopharyngeal and tracheal compartments (Table 1). For β, we found estimates of 1.8×10 -4 and 2.5×10 -5 mL/virion/day(p <10 -4 ) and p was estimated to 1.9×10 4 and 3.6×10 4 virions/cell/day (p<0.05) in nasopharynx and trachea, respectively. However, viral production depends on the assumed initial number of target cells, T 0 . For nasopharyngeal and tracheal compartments, we assumed a total number of epithelial cells of 2.25×10 7 and 1.25×10 8 cells. Assuming that only 0.1% of cells express ACE2 receptors at their surface, we have 2.25×10 4 and 1.25×10 5 susceptible cells in nasopharynx and trachea, respectively. This leads to 4.4×10 8 and 4.5×10 9 virions/mL/day in nasopharynx and trachea respectively. The loss rate of infected cells, δ, was not found different between the two compartments and estimated to 1.88 d -1 corresponding to an infected cell half-life of 9 hours.
These parameter estimates allow us to derive the basic reproduction number R 0 corresponding to the number of infected cells generated by a single infected cell at the beginning of the infection. We found R 0,N equal to 5.9 and R 0,T to 4.0 . Together with the high inoculum in the trachea (see methods), the viral load scarcely increased in the trachea while clearly increased until 3 dpi in the nasopharynx (Fig. 3). One can also derive the viral burst size N corresponding to the number of viruses produced by an infected cell over its lifespan. We found N N = 18,300 and N T = 9,600 .

Sensitivity analysis
Although values at which the viral clearance c and the eclipse rate k were xed based on the current literature (10 d -1 and 3 d -1 respectively, see material and methods). We tested models with different values of these parameters. We considered 9 models resulting from the combination of models with c = 5; 10 or 20 d -1 and k = 1; 3 or 5 d -1 . Models assuming larger eclipse phase duration degraded the t to the data with a substantial increase in the Bayesian Information Criterion (BIC) while other models were broadly undistinguishable from the reference model and within ~3 BIC points. Thus, we kept c = 10 d -1 and k = 3 d -1 for further analysis.

Immune markers during SARS-CoV-2 infection
Among the 30 cytokines tested, 6 greatly varied during the infection and peaked at 2 dpi, namely CCL11, CCL2, IFN-, IL-15, IL-1RA and IL-2 ( Fig. 4). Spearman correlation tests were performed between the area under the cytokine curve and the predicted log 10 AUC viral load, but there was no trends suggesting an association between cytokine and viral loads (Fig. 5).
We further considered potential association between IFN-and viral dynamics. For that purpose, we tested variations of models (Model 1 to 4, see supplementary information le 2) where IFN-levels over time, noted F(t), could modulate viral kinetic parameters. We considered 4 potential scenarios 11,15,16 where IFN could i) reduce infectivity by a factor 1/(1+φF(t)), ii) preserve φF(t) target cells from infection, iii) increase the loss of infected cells by a factor (1+φF(t)) and iv) lower viral production by a factor 1/(1+φF(t)). Overall, models assuming an increased infectivity or the formation of preserved cells showed a moderate gain in BIC (< 5 points) but did not provide better individual ts of the data. Thus, we did not consider those models for further analysis. Models, parameters estimates and individual ts of these models are provided in supplementary information le 2.

Simulation of a human infection
Given the physiopathology of a human infection by SARS-CoV-2, viral replication will start after the rst cell is infected in both nasopharynx and trachea. By taking into account this aspect of the infection and scaling the compartments to human size (see methods), we predict that the viral load peak in human is 1.5×10 7 ( ) and 9.7×10 6 ( copies/mL in nasopharynx and trachea, respectively. In addition, the viral loads were predicted to be detectable during 10.4 and days in nasopharynx and trachea, respectively (Fig. 6). Lastly, we calculated the generation time T g . Assuming that the infectiousness of the virus is proportional to the infectious viral loads, the generation time corresponds to the mean duration between infection of the index case and the subsequent infection of a contact. We found T g,T = days and T g,N = days suggesting that infectiousness was maximal within the rst week of infection, but could still occur 2-3 weeks after infection. in this experimental model (10 6 pfu), explains that tracheal viral loads barely increased post-infection and that nasopharyngeal viral loads rapidly peaked at 3 dpi. After peak viral load, the rapid loss rate of infected cells was su cient to explain that clearance of the virus occurred around day 7 in both compartments, and we did not nd evidence in this model for a role of an immune response mediated by cytokines in accelerating the viral clearance.

Discussion
Although the number of animals and the very detailed kinetic data allowed a precise estimation of parameters, some hypotheses were made that will need to be con rmed. First, the estimation of R 0 depends on the estimation of the viral infectivity rate, β. In the experimental design, animals were inoculated with a total number of 10 10 RNA copies and 10 6 pfu. Although only a small proportion of established the infection in the considered compartments (total V 0 equal to 8.10 7 , see Table 1), this nonetheless implies that a large number of cells are probably readily infected upon inoculation, making it di cult to precisely estimate the rate of infection by circulating virus. This led to uncertainty in the estimate of R 0 with a 95% con dence interval ranging between 1 and 18 in both compartments. Third, the number of target cells was xed to 0.1% of the total susceptible cells. This hypothesis is supported by the fact that in humans 0.1% of alveolar type II cells expressing the ACE2 receptor, gate for SARS-CoV-2 to enter host cells 17,18 . Such estimate is unknown to our knowledge in cynomolgus macaques. Although the estimate of T 0 should therefore be taken with caution, we nonetheless note that it is consistent with the burst size of the virus. Indeed, we estimated the burst size of infectious virus to be equal to 9.6 and 18.3 in the tracheal and nasopharyngeal compartments, respectively, which is consistent with the estimate of secondary infection, R 0, estimated in the two compartments. Fourth, we relied only on measures in both compartments, which may not re ect the kinetics in the lower respiratory tract. It is in particular possible that the kinetics of both the virus and the immune response may be different in the lung, and that both cytokine responses and the lesions as observed by CT scans may be associated with viral loads in the lungs. In our experiments the rst viral load measurements in bronchoalveolar lavages (BAL) was made at 6 dpi, and were all below the limit of detection at the next available data point at day 14, precluding a more detailed analysis of the kinetics in the lungs.
We also evaluated the potential role of immunity in viral resolution by investigating several immune response models. We considered 4 models of the immune response involving IFN-concentrations 11,19 but none of them showed any strong improvement of data tting. This is not surprising since none of the cytokines measured here were related to viral dynamics (Fig. 5). As pointed out above, this may be due to the fact that we only relied on data in the URT but may also re ect the fact that the infection in NHPs was mild. This is consistent with data obtained in patients with an asymptomatic infection, in which the immune response and the cytokine response remained low throughout the infection period 20 .
Finally, we evaluated whether the parameters found in macaques could be relevant for human infection.
To do this, we used the same parameter values and we only modi ed two parameters to be re ective of human infection. First we used the same number of target cells as previously estimated in human URT 21,22 ; second, we assumed that the infection was initiated by only one infectious virus (see more discussions on the number of virus in droplets during small speeches 23 ). Interestingly, our model could well recapitulate a typical viral kinetic curve. Indeed, we found that the viral load measured on nasopharyngeal swabs peaked at day 4 (CI 95% = [2.2-8.7]) post infection, which is similar to results found in prospective cohorts of asymptomatic or mild infections 24 . Interestingly, the predicted duration of viral shedding was 10.4 days (CI 95% = [8.9-16.6]), remarkably in line with estimates of 9 days 25 and 12 days 9 , but lower than other estimates of about 20 days 20,24 . Moreover, we could use the model to predict the course of infectious virus. We predicted that infectious virus was below the limit of detection at day 7, which is also earlier than the estimates of Van

Experimental procedure
Our study includes 31 cynomolgus macaques (16 male, 15 female) infected with 10 6 pfu of a primary isolate of SARS-CoV-2. Each animal received 5 mL of the total inoculum: 4.5 mL were injected by the intra-tracheal route and 0.5 mL by the intra-nasal route. Nasopharyngeal and tracheal swabs were collected longitudinally at days 1, 2, 3, 4, 5, 7, 9, 11, 13, 16 and 20 post-infection (pi). In parallel, bronchoalveolar lavages were collected at day 6, 13 and 20. Viral RNA levels were assessed in each sample using a real-time PCR, with 8540 and 180 copies/mL as quanti cation and detection limits, respectively.
The original study included 6 groups treated either by a high dosing regimen (Hi) of HCQ (90 mg/kg loading dose and 45 mg/kg maintenance dose) ± AZM (36 mg/kg of at 1 dpi, followed by a daily maintenance dose of 18 mg/kg), a low dosing regimen (Lo) of HCQ (30 mg/kg loading dose and 15 mg/kg maintenance dose) or a vehicle (water) as a placebo. Among the 23 treated animals, 14 were treated at 1 dpi with a Lo (n = 4), Hi (n = 5) dose of HCQ or HCQ + AZTH (n = 5). A late treatment initiation was also investigated in 4 animals receiving a Lo dosing regimen. Finally, 5 animals received HCQ at the Hi dose starting at day 7 prior to infection as a pre-exposure prophylaxis (PrEP). Given that HCQ was not associated to an antiviral effect 14 , we pooled the data from either treated and untreated groups. Thus, we consider all animals as controls in our primary analysis but we investigated a potential effect of HCQ in reducing the viral production p (see supplementary information le 1).

Viral dynamic model of SARS-CoV-2
Model describing nasopharyngeal and tracheal swabs We developed a viral kinetic model that simultaneously characterizes nasopharyngeal and tracheal SARS-CoV-2 infection kinetics (Fig. 2). In this model, both nasopharynx and trachea were considered as distinct compartments, each one described by a target cell limited models 15,16 . In the following, only a generic description of the model is presented and subscripts were omitted. Target cells (T) are infected by infectious viruses ( ) at an infectivity rate β. After infection, infected cells enter an eclipse phase (I 1 ) where they do not produce virions and enter into productively infected cells (I 2 ) at a rate k. I 2 cells produce p virions daily and are lost at a per capita rate δ. A proportion µ of the virions are infectious ( ) and the remaining (1-µ) are noninfectious viruses ( ). Finally, both and are cleared at a rate c.

Fixed parameters
Because not all parameters could be properly estimated based on total RNA data only, some parameters had to be xed based on the experimental setting or the current literature. First, the term pT 0 is the only identi able quantity in our model, were T 0 is the initial amount of susceptible cells at the beginning of the infection. Given the estimated surface of the trachea (S T = 9 cm 2 ) and nasal cavities (S N =50 cm 2 ) in cynomolgus macaques (data not shown) and the apical surface of one epithelial cell s=4 10 -7 cm 2 /cell 21 , one can calculate the number of target cells exposed to the virus in the nasopharynx and the trachea and , respectively. However, only a small fraction of these cells may express the ACE2 receptors needed for SARS-CoV-2 to enter a cell. Hence, assuming a proportion of 0.1% of cells expressing the ACE2 receptor, we set and cells in nasopharynx and trachea respectively. Second, we supposed the proportion of infectious viruses µ remained constant over time and equal to 0.001. Indeed, infectious titers measured by cell culture on tracheal swabs taken at 3 dpi (see details in supplementary information le 3) ranged from 1.2 to 3.2 log 10 TCID 50 /mL i.e. 1,000 to 100,000 times lower than total RNA (Fig. 7).Third, we supposed that the eclipse phase duration was equal in the nasopharynx and in the trachea, and we set k=3 d -1 as the viral loads have been shown to increase 8h after infection 29 . Fourth, the viral clearance c was assumed equal whatever the infectious nature of the virus and its location. Thus, we xed c to 10 d -1 as the clearance is assumed to be fast in other respiratory infections 16 where γ indicates the xed effects and the random effects, which are supposed to follow a normal distribution of mean zero and standard deviation ω. The residual error e ijk is assumed to follow a normal distribution of mean zero and constant standard deviation .
Standard errors were calculated by drawing parameters in the asymptotic distribution of the parameter estimators. For each parameter, we calculated the 2.5 and 97.5% percentile to derive the 95% con dence interval.

Modelling strategy
In order to reduce the number of parameters, we tested whether the transition rate between nasopharynx and trachea g could be set to 0. Then, we tested the possibility for parameters β, δ, p and V 0 to be equal in both nasopharynx and trachea and tested if their variances could be set to 0. To do this, we used a backward selection procedure and stopped once the BIC did not decrease anymore. Lastly, based on the Empirical Bayes Estimates (EBE), we screened the random effects for correlations. Only correlations with a Pearson's coe cient >0.8 were implemented in the model.

Plasma cytokine analysis
In all 31 macaques, the concentration of 30 cytokines were measured at 0, 2, 5, 7 and 9 dpi. Among them, CCL11, CCL2, IFN-, Il-15, Il-1Ra and Il-2 were of particular interest as their kinetic changed during the infection (Fig. 4). To identify the cytokines to be incorporated in the model, we correlated the are under the cytokine curve (AUC, calculated by the linear trapezoidal method) with the AUC of log 10 viral loads predicted by the model and calculated Spearman correlations (Fig. 5). Eventually, Spearman correlations with p-value<0.1 were further investigated and implemented in the viral dynamic model. In addition, we tested models involving IFN-, a central element of the immune response against SARS-CoV-2 in humans 27 .

Parameter estimation
Parameters were estimated with the SAEM algorithm implemented in MONOLIX® software version 2018R2 allowing to handle the left censored data 30 . Likelihood was estimated using the importance sampling method and the Fisher Information Matrix (FIM) was calculated by stochastic approximation. Graphical and statistical analyses were performed using R version 3.4.3.

Simulation of the infection time course in humans
Assuming that the number of viruses entering a host is low in humans, the infection would start after a virus had infected its rst cell. Thus, we scaled the compartments to human size considering a total URT cell count of 4x10 8 and 30 mL of volume 16,22 and assumed that the infection started after infection of one productive infected cell in the nasopharynx and in the trachea (i.e. 1 cell). We then computed the median and associated CI 95% of peak viral load and viral shedding, as well as the generation time T g corresponding to mean duration between a primary case and its associated secondary infections assuming that the infectiousness is proportional to the infectious viral load levels .  Figure 1 Nasopharyngeal and tracheal SARS-CoV-2 viral loads in infected cynomolgus macaques treated by hydroxychloroquine ± azithromycin. Dashed lines represent the lower limit of detection (LOD) and lower limit of quanti cation (LOQ).

Figure 2
Schematic representation of the model of nasopharyngeal and tracheal viral loads. Parameters c, k and µ were assumed identical in nasopharynx and trachea while we tested wheter β, p and δ could be different in both compartments (see methods) Page 17/21 Figure 3 Nasopharyngeal (blue) and tracheal (red) individual predicted dynamics by the model described in equations (1) to (5). Full dots are the quanti able viral loads and crosses the observation below the limit of quanti cation. The dotted line represents the limit of quanti cation and the dashed line the limit of detection.
Page 18/21   Relationship between infectious titers/total viral load ratio and total viral loads. Each dot corresponds to a tracheal swab at 3 dpi used to measure the infectious titers (see supplementary information le 3) and the total viral loads. Green, red and purple dashed lines correspond to ratios of 10-5, 10-4 and 10-3 respectively.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.