SARS-CoV-2 viral dynamics in non-human 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 for which nasopharyngeal and tracheal viral load were frequently assessed. We identified that infected cells had a large burst size (>104 virus) and a within-host reproductive basic number of approximately 6 and 4 in nasopharyngeal and tracheal compartment, respectively. After peak viral load, infected cells were rapidly lost with a half-life of 9 hours, with no significant association between cytokine elevation and clearance, leading to a median time to viral clearance of 10 days, consistent with observations in mild human infections. Given these parameter estimates, we predict that a prophylactic treatment blocking 90% of viral production or viral infection could prevent viral growth. In conclusion, our results provide estimates of SARS-CoV-2 viral kinetic parameters in an experimental model of mild infection and they provide means to assess the efficacy of future antiviral treatments.

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 December 2020 in more than 1,600,000 deaths [1]. Fortunately, the majority of infections do not lead to hospitalizations [2], and the vast majority of subjects infected with SARS-CoV-2 will experience asymptomatic or pauci-symptomatic infection characterized by specific (anosmia) or general symptoms (fever, fatigue) [3-5]. In other acute or chronic viral diseases (HIV, HCV, influenza), the characterization of viral load kinetics has played an important role to understand the pathogenesis of the virus and design better antiviral drugs [6-8]. In the case of SARS-CoV-2, viral kinetics has been found to be associated with mortality in hospitalized patients [9] but the association in less severe patients is unclear. This is due to the fact that many studies rely on large cross-sectional analyses with few patients having serial data points or, in contrary, on detailed small series of patients [10-12]. 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 [13-15].

Here, we used data generated on nonhuman primates (NHP) and treated with hydroxychloroquine (HCQ) ± azithromycine (AZTH) in either pre- or post-exposure prophylaxis [16] to develop a mathematical model of SARS-CoV-2 infection. Although our analysis did not reveal any significant antiviral efficacy of HCQ, the large data generated (31 cynomolgus macaques with frequent measurements of both nasopharyngeal and tracheal viral loads) allowed us to characterize in detail the key parameters driving the viral dynamics and to evaluate putative immune response mechanisms during the infection.

Methods

Ethics statement

Cynomolgus macaques aged 37–40 months and originating from Mauritian AAALAC certified breeding centers were used in this study. Animals were housed under BSL-2 and BSL-3 containment when necessary (Animal facility authorization #D92-032-02, Prefecture des Hauts de Seine, France) and in compliance with European Directive 2010/63/EU, the French regulations and the Standards for Human Care and Use of Laboratory Animals, of the Office for Laboratory Animal Welfare (OLAW, assurance number #A5826-01, US).

The protocols were approved by the institutional ethical committee “Comité d’Ethique en Expérimentation Animale du Commissariat à l’Energie Atomique et aux Energies Alternatives” (CEtEA #44) under statement number A20-011. The study was authorized by the “Research, Innovation and Education Ministry” under registration number APAFIS#24434-2020030216532863v1.

Experimental procedure

Our study includes 31 cynomolgus macaques (16 male, 15 female) infected with 106 pfu (corresponding to 1010 total RNA copies) of a primary isolate of SARS-CoV-2 (BetaCoV/France/IDF/0372/2020). Each animal received 5 mL of the total inoculum: 90% (4.5 mL) were injected by the intra-tracheal route and 10% (0.5 mL) by the intra-nasal route [17,18]. Nasopharyngeal and tracheal swabs were collected longitudinally at days 1, 2, 3, 4, 5, 7, 9, 11, 13, 16 and 20 post-infection (pi) and eluted in 3 mL of Universal transport medium (Copan) or Viral Transport Medium (CDC, DSR-052-01). Viral RNA levels were assessed in each sample using a real-time PCR, with 8540 and 180 copies/mL as quantification and detection limits, respectively ()

Nasopharyngeal and tracheal SARS-CoV-2 viral loads in infected cynomolgus macaques treated by a placebo or hydroxychloroquine ± azithromycin.

Dashed lines represent the lower limit of detection (LOD) and lower limit of quantification (LOQ).

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) ± AZTH (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) (see )

Viral dynamic model of SARS-CoV-2

Model describing nasopharyngeal and tracheal swabs

Given that HCQ was not associated to any antiviral effect [16], we pooled the data from either treated and untreated groups. Thus, we consider all animals as controls in our primary analysis but we also investigated a potential effect of HCQ (see details in ). 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 [13,19,20]. In this model, susceptible cells (T) are infected by an infectious virus (VI) at a rate β and generate non-productive infected cells (I1). After an eclipse phase of duration 1/k on average, infected cells start to be productive (I2) and produce both infectious and non infectious viruses at a rate μp and (1-μ)p respectively. Productive infected cells die at a rate δ and viruses are cleared at a rate c (see ). The model is given by ordinary differential Eqs (1) to (5) where subscript X denotes the compartment of interest either nasopharynx or trachea. The basic reproduction number is and the burst-size is .

Fixed parameters

Because not all parameters can be estimated when only total viral RNA are measured, some parameters had to be fixed based on the experimental setting or the current literature. As only the product p×T is identifiable, we chose to fix T the initial concentrations of susceptible cells, as follows. We measured the surfaces (S) and volumes (V) of the nasopharynx and the trachea in one euthanized animal and obtained V = 6.3 mL; V = 1.2 mL; S = 50 cm2; S = 9 cm2. Assuming an apical surface of one epithelial cell of s = 4 10−7 cm2/cell [21], the number of target cells exposed to the virus in the nasopharynx and the trachea are and , respectively. Only a fraction of these cells expresses both angiotensin-converting enzyme 2 (ACE2) and the type II transmembrane serine protease (TMPRSS2), and we fixed this proportion to 0.1% i.e. T = 1.98×104 and T = 1.88×104 cells in nasopharynx and trachea, respectively. Also, we supposed that the eclipse phase duration was equal in the nasopharynx and in the trachea, and we set k = 3 d-1 based on in vitro studies showing that virus release occurs 8 hours after infection [22]. Third, we supposed the proportion of infectious virus μ remained constant over time. TCID50 is a proxy for the infectious viruses. shows that the ratio TCID50/total viral loads is 10−5 (range 10−6–10−4). However, not all the infectious viruses present in the samples may be detected in Vero-E6 cells experiments. Hence, we fixed μ to 10−4, the upper bound of the observed ratio of infectious virus. Finally the viral clearance c was fixed to 10 d-1, consistent with the rapid viral clearance rate of influenza virus [20]. Sensitivity analyses exploring the consequences of those different assumptions were performed, with c in a range 5–20 d-1, μ in a range 10−5–10−3 and the proportion of target cells being equal to 1% (see ). Thus, the estimated parameters in each compartment were V0, p, β and δ.

Relationship between TCID50 and viral loads. Each symbol corresponds to a tracheal swab obtained at 3 dpi.

A) Viral load levels in samples with no detectable infectious virus (circles) and those with detectable infectious virus (triangles). B) Correlation between the viral loads and the infectious virus, as measured by TCID50. Red, green and purple dashed lines represent putative ratios of infectious virus of 10−4; 10−5 and 10−6 respectively.

Statistical model

The structural model used to describe the observed log10 viral loads Y of the ith animal at the jth time point in the kth compartment (k = 1 for nasopharyngeal or k = 2 for tracheal) is

Where f is the function describing the total viral loads dynamics over time (V(t)+V(t)); θik is the vector of parameters of subject i in the k compartment and e is the additive residual error. Individual parameters θik are supposed to follow a log-normal distribution with a median value that depends on the compartment: where γ indicates the fixed effects and η the individual random effects, which are supposed to follow a normal distribution of mean zero and standard deviation ω, and exp(α) is the vector of the ratio values between the nasal and tracheal compartments. The residual error eijk 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% confidence interval.

Modelling strategy

As a first step, we considered nasopharynx and trachea parameters as two distinct compartments and we then tested whether the virus could migrate from one compartment to the other at a constant first order migration rate g. Parameter g was set to 0 if it did not improve the fit to the data or could not be precisely estimated. In order to reduce the number of remaining parameters, we fixed the ratio of effective inoculum, V0,/V0,, to 9. Then we tested the possibility for estimated parameters to be equal in both nasopharynx and trachea (e.g., α = 0) and tested then if their variances could be set to 0 (e.g., ω = 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 coefficient >0.8 were implemented in the model. All models tested are presented and compared in .

Viral titers determination on Vero-E6 cells

Vero-E6 cells (mycoplasma-free) were seeded in 24-well plates (2×105cells/well) and cultured in DMEM (Thermo Fisher Scientific) containing 1% PS (Penicillin 10,000 U/mL; Streptomycin 10,000 μg/mL) supplemented with 5% FBS (Foetal Bovine Serum) and incubated at 37°C in the presence of 5% CO2. The next day, cells were inoculated in triplicate with 100μL per well of the tracheal swab dilutions (1:2, 1:10, 1:50,1:250) in DMEM, 1% PS, 0.1% TPCK trypsin and incubated for 1 hour at 37°C. After removal of the inoculum, 1mL of DMEM, 1% PS, 0.1% TPCK trypsin was added in each well before incubation at 37°C in the presence of 5% CO2 for 72 hrs. The presence of a cytopathic effect (CPE) was visualized under the microscope and the TCID50 i.e. the tissue culture infective dose leading to 50% of the maximal cytopathic effect,was calculated using the method of Reed and Münch. All experiments were conducted under strict BSL3 conditions.

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. 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 log10 viral loads predicted by the model.

Models assuming a compartment for the innate immune response

We considered additional models incorporating a compartment for an antigen-dependent immune response, F, given by . In these models, F could either i) reduce viral infectivity ii) and iii) put target cells into a refractory state, iv) reduce viral production and v) increase the loss of infected cells (see ).

Simulations of a prophylactic treatment

We used the median estimated parameter values of the model to simulate the effects of an antiviral treatment, initiated before infection, on the viral kinetics. We explored the effects of 10 to 100 fold lower viral inoculum (corresponding to doses of 105 and 104 pfu) as well as the effects of drugs acting on viral production, viral entry or the proportion of infectious virus. For each scenario, we also calculated the 95% confidence interval of the median time to viral clearance, using the method described above.

Parameter estimation

Parameters were estimated with the SAEM algorithm implemented in MONOLIX software version 2018R2 allowing to handle the left censored data [23]. 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.

Results

SARS-CoV-2 viral kinetics

In our experiment, cynomolgus macaques developed a rapid infection, with viral loads peaking 2 days post infection (dpi) in both nasopharyngeal and tracheal compartments. Afterwards, both nasopharyngeal and tracheal viral loads rapidly declined exponentially, with a similar median rate of 1.9 d-1, corresponding to a daily reduction of 0.8 log10 copies/mL (). Because the viral load peaked later and higher in in the nasopharynx than in the trachea (7.9 and 7.2 log10, respectively), the median time to unquantifiable viral load nonetheless occurred later in the nasopharynx than in the trachea (9 and 7 dpi, respectively).

Overall, the area under the viral load curve (AUC) was larger in the nasopharynx than in the trachea (45 vs 38 log10 copies.day/mL, p<10−4). In addition to viral RNA, we also measured virus titers in Vero-E6 cells using tracheal swabs obtained at 3 dpi (see Methods). All samples contained more than 4 log10 copies/mL however viral growth was observed only in those having more than 6 log10 copies/mL (). In those for which viral culture could be obtained, the ratio of TCID50 (median tissue culture infective dose) to the total number of RNA copies ranged between 10−4 to 10−6 ().

Viral dynamic model

Importantly, there was no antiviral effect in animals receiving various doses of HCQ compared to those receiving vehicles, even after adjustment on the mean HCQ exposure [16]. Thus, the effect of treatment was neglected as a first approximation. We later challenged this hypothesis but found no significant effect of HCQ in reducing viral production or viral infectivity (see ).

Physiologically, viruses can migrate from the nasopharynx to the trachea through respiration and movements of ciliary cells at the mucosal surface. Thus, we tested the possibility for viruses to move from nasopharyngeal to tracheal compartment and vice versa by linking both with a bidirectional rate constant g. However, possibly due to data paucity, the parameter g was not significantly different from 0 (CI95% included 0) and was therefore set to 0 in the following. Then, using a backward selection procedure we found that the infectivity rate β and the viral production p were different between nasopharyngeal and tracheal compartments (). The final model well fitted the data and allowed the estimation of several parameters of the infection (). We estimated the effective initial viral load to 1.7×106 and 8.0×107 copies/mL in tracheal and nasal compartments, respectively, which corresponds to a total inoculum of approximately 108 total RNA, i.e., about 1% of the total injected dose (see Methods). We found estimates of β of 1.2×10−3 (CI95% 0.5×10−3−2.8×10−3) and 1.8×10−3 (CI95% 0.6×10−3−5.6×10−3) mL/virion/day (p<10−4) in nasopharynx and trachea, respectively, while p was estimated to 4.8×104 (CI95% 1−8.8×104) and 2.2×104 (CI95% 0.7−4.1×104) virions/cell/day (p<0.05). Consequently the product p×T0 was equal to 9.5×108 (CI95% 2.0–17×108) and 3.9×108 (CI95% 1.2–17.7×108) virions/mL/day in nasopharynx and trachea, respectively.

Nasopharyngeal (blue) and tracheal (red) individual predicted dynamics by the model described in Eqs (1) to (5). Full dots are the quantifiable viral loads and crosses the observation below the limit of quantification. The dotted line represents the limit of quantification and the dashed line the limit of detection.

The loss rate of infected cells, δ, was not found different between the two compartments and estimated to 1.9 (CI95% 1.6–2.3) d-1 corresponding to an infected cell half-life of 9 (CI95% 7–13) hours. These parameter estimates allow us to derive the basic reproduction number R corresponding to the number of infected cells generated by a single infected cell at the beginning of the infection. We found R equal to 5.6 (CI95% 1.3–21) and R equal to 3.8 (CI95% 0.7–18) in the nasopharynx and the trachea, respectively. 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 = 25,000 (CI95% 5,900–72,000) virions and N = 11,000 (CI95% 4,000–19,000) virions.

Sensitivity analysis

In our main analysis, the viral clearance c and the proportion of infectious virus μ were fixed. We tested the robustness of our results to different values of these parameters (see ) and we used model averaging (MA) to compute the averaged parameters values following a methodology presented in [24]. Overall, models were broadly undistinguishable as they provided a Bayesian Information Criterion within a 4 points range. Model averaged parameter estimates of δ, R0,N and R0,T were equal to 1.9 (CI95% 1.5−2.3) d-1, 6.8 (CI95% 1.3−29) and 4.4 (CI95% 1−26), respectively.

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) but there was no association between cytokine and viral loads (). A model assuming an effect of the immune response in cell protection resulted in a reduced BIC of 6 points. However, the gain in fitting criterion was entirely due to 3 individuals (see ) and led to more uncertainty in parameter estimates due to increased complexity (see Table B in ). Moreover none of the cytokines measured during the experiments, including IFN-a, showed a correlation with viral dynamics (). Thus, overall, there was no clear evidence in these data that a more complex model improved the understanding of viral dynamics over a simple target cell model.

Correlation between the individual predicted AUC log10 viral load and cytokine AUC.

None of the 6 cytokines that increased during the infection was significantly correlated with the predicted viral load AUC in the nasopharynx or the trachea.

Expected profiles with prophylaxis treatments

We used estimated parameter values of the model to simulate the effects of an antiviral treatment, initiated before infection, on the viral kinetics. We explored different viral inoculum (ranging from 104 to 106 pfu), drugs mechanisms of action (blocking viral production, viral entry, or infectious virus production), with different levels of antiviral efficacy. As the conclusions were not sensitive to the mechanism of action, we present below the results for a prophylactic drug blocking viral entry (). Other mechanisms of action are presented in

Median viral kinetic profiles in the nasopharynx (blue) and the trachea (red) according to the inoculum size and the level of an antiviral initiated in prophylaxis and blocking viral entry β. Treatment efficacy of 0% (no treatment, solid line), 90% (dashed line) and 99% (dotted lines) were considered. Point-dotted line represents the threshold below which no virus could be cultured in vitro ().

In both compartments, the median time to viral clearance increased with lower doses of inoculum, and ranged from 10.1 (CI95% 8.4–14.7) days with 106 pfu to 11.7 (CI95% 9.5–19.0) days with 104 pfu in the nasopharynx. A 90% effective antiviral treatment administered upon infection would dramatically reduce peak viral load in all scenario and maintain virus below the threshold level of infectivity of 6 log10 copies/mL (). A 99% effective antiviral treatment could in addition abrogate viral load viral load, with time to viral clearance ranging from 3.5 (CI95% 1.9–6.2) days with 106 pfu to 0.5 (CI95% 0.2–2.5) days with 104 pfu in the nasopharynx.

Discussion

We here developed a mathematical model for SARS-CoV-2 viral dynamics using nasopharyngeal and tracheal swabs obtained in 31 infected macaques. Using this model we could estimate key parameters of virus pathogenesis, in particular the viral infectivity rate (equal to 1.2×10−3 and 1.8×10−3 in tracheal and nasopharyngeal compartments, respectively), and the loss rate of infected cells, estimated to 1.9 d-1 in both compartments (e.g., a half-life of 9 hours). Consequently, we estimated the number of secondary cell infection resulting from one infected cells, R, to approximately 4 and 6 in tracheal and nasopharyngeal compartments, respectively. This value of R together with the large viral inoculum used in this experimental model (106 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 sufficient to explain that clearance of the virus occurred around day 10 in both compartments, and we did not find evidence in this model for a role of an immune response mediated by cytokines in accelerating the viral clearance.

Although the number of animals and the very detailed kinetic data allowed the precise estimation of several parameters, some limitations exist. First, animals were infected with a large inoculum (106 pfu, i.e. 1010 RNA copies) which rapidly saturate target cells making the estimation of β less robust and difficult to distinguish the processes of clearance of the inoculum from those of de novo viral infections. In the future, the analysis of subgenomic RNA, that quantifies intracellular viral transcription, will provide important information to distinguish these two processes and provide a more precise estimate of R0 [25,26]. Second, a number of unknown parameters were fixed to ensure identifiability, in particular the proportion of infectious virus, μ, and the number of target cells available in each compartment, T0. In this model, p×T0 is the only reliable quantity since the estimate of p depends on the number of susceptible cells. We here estimated the number of alveolar type II cells by analysis of the size of the tissues of one euthanized animal and we assumed that 0.1% of these cells could be target for the virus. Although this value led to more coherent parameter estimates (see ), this is higher than the value of 1% of cells expressing the ACE2 receptor and serine protease TMPRSS2 found in humans [27,28]. Specific experiments will be needed to estimate the proportion of target cells in the URT of macaques. Third, we relied only on measures in both compartments of the upper respiratory tract, which may not reflect 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 first 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 aimed to evaluate whether more complex models including an antigen-dependent immune response could improve the data fitting. A model assuming that cells could be put in an antiviral state improved the BIC (see ), however this improvement was due to 3 individuals, and this came to the expense of a deteriorated precision of parameter estimates. Further, none of the observed cytokines were found associated with viral dynamics (), suggesting that this improvement in data fitting was not supported by our data. In addition, none of the animals had detectable antibodies until day 7, and only 25% had detectable antibodies by day 14, suggesting that the humoral response played a minor role in viral clearance. The role of the immune response in this experimental model of mild infection is unclear, but our findings are 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 [29]. Accordingly our estimate of the duration of viral shedding was between 10 and 12 days in the nasopharynx, depending on the initial inoculum, very close to the values of 9 days (Diamond Princess [30]) and 12 days [11] estimated in mild or asymptomatic individuals. From an experimental setting, the analysis of viral dynamics in animals infected with different size of viral inoculum could also bring insights on the need to use more complex models [31].

Finally, we used the model to inform on the use of prophylactic drugs in this macaque model. Given our estimate of R0, a 90% effective treatment should be able to prevent virus growth and would maintain the viral load levels below 106 copies/mL at all times, making the chance of detecting infectious virus very limited, consistent with our modeling predictions in humans ([32]).

Effects of hydroxychloroquine.

(DOCX)

Click here for additional data file.

Model file for Monolix.

(TXT)

Click here for additional data file.

Sensitivity analysis.

(DOCX)

Click here for additional data file.

Model building.

(DOCX)

Click here for additional data file.

Immune response models.

(DOCX)

Click here for additional data file.

Simulations.

(DOCX)

Click here for additional data file.

Nasaopharyngeal and tracheal swabs data file.

(TXT)

Click here for additional data file.

22 Nov 2020

Dear Mr. Gonçalves,

Thank you very much for submitting your manuscript "SARS-CoV-2 viral dynamics in non-human primates" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

The reviewers agreed that your work would be a great contribution to the current literature on SARS-CoV-2 and COVID-19. However, the reviewers raised many, in our opinion valid and valuable, points of criticism. We encourage you to address these in your revision as we believe that this strengthen your paper and increase its impact. Especially the point raised by Reviewer #3 on the very high challenge dose in the experiments and the complications that introduces for the estimation of R0 would be important to address. Reviewer #1 asks you to explore additional mechanisms that might be able to give you a better model fits (point 1). Reviewer #4 ask you to further justify some model assumptions (point 5). In our view, the paper would gain if you justified which mechanisms you chose to include and which potential mechanisms you omitted.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Roland R Regoes

Associate Editor

PLOS Computational Biology

Rob De Boer

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: To characterize the SARS-CoV-2 nasopharyngeal and tracheal viral load dynamics in 31 cynomolgus macaques, the authors developed a target-cell limited mathematical model for each tissue compartment. The model includes an eclipse phase and the possibility of production of non-infectious viruses. The model is used to estimate basic reproductive ratios for each tissue compartment and other kinetics parameters. The authors also performed correlation tests between immune markers with viral load and explored competing models that include association of viral load and IFN-α levels. They finally made projections of the effect of >90% effective prophylactic interventions. The manuscript is well written and the model is based on a well-known model of virus dynamics. My biggest concern is related to the accuracy of the model fits in capturing the observed behavior, and therefore, the accuracy of parameter estimates and conclusions based on the model.

Major comments

1. It looks like the viral load observations have different kinetic patterns beyond a viral expansion and clearance, and therefore, the model is not able to capture them. This might be either due to high measurement error or to an additional mechanism that is not in the model (or its extensions using the IFN-α-related models). For example, for the nasopharyngeal observations (blue in Fig 4) there is the possibility of more than one episode/peak (MF6, 13,17,18,19, 21,29) or a plateau/sustained phase where the virus is maintained for several days (MF1, 6, 12 and some of the others mentioned before). In that sense, it is not clear if this phrase in page 5 is completely accurate: “After the peak, both nasopharyngeal and tracheal viral loads rapidly declined exponentially “. It also looks like the model doesn’t capture the nasal dynamics of animals MF4, 5, 10, 26. I suspect that adaptations of this model structure that allows dynamics with a sustained viral load, that may add a couple of parameters could give a significantly lower BIC. I think a larger exploration of models has to be performed to confirm the parameter estimation and conclusions based on the model fits.

2. The paper discusses that the viral peak happens about 2 days post-infection. But there is not a discussion of why the model cannot capture the viral peak and that 2-day time.

3. There is not a discussion of the reasons to select the values of V0 presented in Table 1. There is only a citation to reference 20 (Long Q-X et al, Nat Med) without giving any reason. I wonder how the model would behave using different values of V0. It is possible that allowing for a lower V0 (or other model structures) can capture the dynamics of some of the nasal viral loads that had low levels or its behavior seems delayed (MF4, 5, 10, 24, 26).

Minor comments

4. It was not completely clear if the %RSE in Table 1 was related to the fixed or random effects of the statistical model.

5. There is a missing “)” at the end of page 6 when describing the estimates of βn.

6. The words “we tested” are repeated in the 4th line in the modeling strategy section in page 17.

Reviewer #2: Goncalves et al. previously reported on their results of a study in non human primates and created a dataset of longitudinal virus load data in the upper respiratory tract (URT) for 31 cynomolgus macaques in response to mild SARS-CoV-2 infections. Here, the authors perform mathematical modeling and, using their data, estimate viral kinetics parameters and hypothesize that potent prophylactic treatment can drive the virus to extinction without prior viral growth. The relatively simple model used in their manuscript was compared against more complicated models including virus migration between nasopharyngeal and tracheal compartment, and cytokine mediated immune responses. These models were not found to explain the data significantly better hinting at a minor role of migration and immune response in the URT during mild SARS-CoV-2 infections.

Longitudinal virus load data, in particular data including the pre-symptomatic phase, are extremely rare and hence, despite being limited to the URT and mild infections, the authors study provides an important dataset and estimates for viral kinetics parameters.

However, there are several issues that require attention.

1) You report a virus threshold below which viral culture could not be obtained which is a very interesting result. In Fig. 2B, you compare the ratio of TCID50 and viral load with the viral load (“In those for which viral culture could [be] obtained, the ratio of TCID50 to the number of RNA copies ranged between 10-4 to 10-6 (Fig. 2B). “). It is not clear to me what the importance/interpretation of this graph is. Can you please elaborate on this. Are only the ratios important or is the comparison to the viral load important (e.g. correlated/uncorrelated)? Why does this ratio represent an underestimate of the total infectious virus? (p.16 (“the ratio of titers to RNA copies represent an underestimate of the total infectious virus ”))

2) In the model you split the virus into infectious and uninfectious virus. For clarification it may be beneficial to label viral load as total viral load (V^I + V^NI) in your figures. Furthermore, there are some other issues related to the viral load(s) that require attention (Table 1).

- Why is the unit of initial conditions V_T0 and V_N0 copies and not cp/ml? As is, with the given units for beta and T_T(t=0), T_N(t=0), the units in eq. (1) don’t work out.

- If V_T0 and V_N0 are initial conditions for infectious virus only, then label as V^I_T0, V^I_N0. Furthermore, I do not understand how you obtain the values for V_T0 and V_N0, given that (i) the viral inoculation is 10^10 copies (discussion), (ii) 90% (4.5/5ml) of the inoculum are applied to the tracheal and 10% (0.5/5ml) to the nasopharyngeal compartment and (iii) a fraction mu = 10^(-3) of virions is assumed to be infectious.

Other minor issues

3) In the experiments, is there a reason for splitting the virus inoculation 90% tracheal vs 10% nasopharyngeal? Is this supported by literature of how virus enters the URT?

4) (p.5) “After the peak, both nasopharyngeal and tracheal viral loads rapidly declined exponentially, with a median rate of 0.8 log10 copies/mL every day (Fig. 1). The slope in viral decline was more rapid in the trachea than in the nasopharynx, leading to a first measurement below the limit of quantification (LOQ=8514 copies/mL) 7 and 9 dpi, respectively. “

- Is the 0.8 log10 copies/ml the combined median value for nasopharyngeal and tracheal viral loads?

- Given figure 1 it seems to me that the first measurements below LOQ are much earlier than the reported 7 and 9 dpi. Are these average/median values? Please clarify.

5) (p.6) Equation (2): Typo (VI_X instead of V^I_X)

6) (p.7) “Assuming that only 0.1% of cells express ACE2 receptors at their surface (17) “

- According to figure 1b) in ref. (17) approximately 1% of epithelial cells express ACE2. To my understanding figure 1a) from where I assume you extracted the 0.1% shows all cell types. But since you only model epithelial cells I believe you should work with 1%. Since this will change your T_0 estimates, it will probably also have an impact on p and p*T_0.

7) ”We found R0,N equal to 5.9 (2.0 – 16) and R0,T equal to 4.0 (1 – 13). 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. 4). 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 NN = 22,000 (8,100 – 36,000) virions and NT = 10,000 (4,500 – 19,000) virions.”

- Please refer to the methods section where you give the formulas for R_0 and N.

8) (p.10, Discussion) “we could estimate key parameters of virus pathogenesis, in particular the production rate from infected cells (equal to 1.9 × 104 vs 3.6 × 104 virions/cell/day in tracheal and nasopharyngeal compartments, respectively)”

- You say twice (in the results and the methods) that only p*T0 can be identified. For clarification you should state in the discussion that these estimates are obtained for fixed T0. This seems particularly important to me since you claim a certain degree of uncertainty in your T0 estimate on p.11.

9) (p.16, Statistical model) How is the function f defined? (How) does it relate to model equations (1)-(5)? Please explain.

10) (Supplementary Information file 2) “As the viral clearance c and the eclipse phase rate µ could not be estimated from the data”.

- Mu is not the eclipse phase but the fraction of infectious virus.

11) The manuscript should be checked for typos.

Reviewer #3: In this study Goncalves et al analyze the kinetics of viral load during the entire course of SARS-CoV-2 infection in non-human primates. They describe the dynamics with a mathematical model and estimate model parameters by fitting to data. They also collect longitudinal levels of a panel of cytokines involved in the innate immune response and test for association between these markers and viral load kinetics.

There are many very strong and unique aspects of this study, which I think will make it quite valuable for the field. Longitudinal within-host viral load kinetics are still very poorly characterized for COVID-19, since it is very difficult to catch human subjects early in infection (before symptom onset) and sample them frequently. For this reason, the full timecourse of viral load, its variation between individuals, and its association with the disease course are still not clear. Thus, animal models are still very useful. This study represents a very large cohort of animals with very frequent sampling and a large number of biomarkers measured, so it helps fill in some of these gaps. The authors use a mechanistic mathematical model to describe the data, as opposed to just fitting a simple curve, which has the benefit that it then be used to examine how kinetics might change under different parameter perturbations. They use a rigorous statistical method for group-level fitting, which make it easy to then use the model to simulate expected inter-patient variance in kinetics. The availability of so many immunological markers measured alongside viral loads is unique, and the authors use some very nice methods to test for the impact of these markers on viral kinetics. The Discussion section does a very thorough job of describing the many caveats of their work.

I think there are two major limitations to this work. I do not think the authors can or should do anything differently to address them, I am simply pointing them out to help with the editorial decision.

The first is about the study design in the animal model. It is not clear how much this is helping us understand human infection, because the animals seem to be inoculated with an extraordinarily high dose of virus, and in nearly all animals the entire upslope of the viral load is missing (just like it generally is in human data collected post-symptom onset). We know for many viral infections inoculum size influences disease course, so its possible this high-dose infection model might be qualitatively different than a more realistic infection scenario. And, with the upslope missing, it is nearly impossible to accurately estimate R0 - for infection models within-host or at the population level, it is the early exponential growth phase that provides the information that allows R0 to be an identifiable parameter combination (assuming you separately know/can estimate the lifespan of infected cells). Also, I worry that this initial “overdose” of virus is obscuring some inter-patient variation that exists in natural infection. In humans, there seems to be considerable variation in the time to symptom on set (and therefore likely also in viral load peak), whereas in this model viral load peaks extremely quickly and at a similar time in all animals, so it is unclear how relevant this model is.

The second limitation is that there is really not much about this data that motivates the need to apply a mathematical model, and even more specifically the quite complicated model the authors used. I myself am a modeler and do not need to be convinced of the value of models in general … it’s just that here it sort of seems like a model was used just for the sake of using a model. In the paper, Figure 1 is shown, and then the authors launch into a complex model description. But why? When I look at Figure 1, nothing suggests to me you would need or want to apply any sort of mathematical model to this data. It just looks like decay over time. If you are going to write a whole paper about applying a viral dynamics model to data, it would be helpful to provide the reader with some motivation for that! And the model seems relatively complex - again, you need to motivate this. What is the biology you think is going on and why did you decide to capture it this way? How are the nasopharyngeal and tracheal compartments related, physiologically? How would you possibly have enough information to identify migration of cells? Why did you feel the need to include infectious and non-infectious virus? What is the biology behind this and why do you think it’s necessary to explain your data? Why do you think this process is target cell limited? Where did these #s of cells come from? Why are non-productive cells needed? Then looking at Figure 4, honestly the model does not seem to be doing a particularly amazing job explaining any of this data. Mostly, the data seems to be all over the place and barely following a trend. The model seems to be missing the peak in most animals, and not accurately capturing the early kinetics in the few animals for which this info is available.

Despite these concerns, I do still feel like the results of the paper are useful, because we are still very early on in our understanding of this virus and understandably animal models are not yet perfected, and the model the authors use is at least reasonable so can provide a starting point for other work that probes these mechanisms more carefully, and for now can at minimum serve as a calibrated simulation tool for recreating within-host dynamics.

A few more specific minor comments, mainly about typos/grammar:

* Author summary: Typo in these sentences

* "We found that viral load rapidly declined after peak viral load [despite we found] no association between model parameters and immune markers.”

* The within-host reproductive basic reproduction number was estimated to BE 6 and 4 in nasopharynx and trachea suggesting that a prophylactic therapy blocking viral entry or production with 90 efficacy could be sufficient to prevent viral growth and peak viral load.

* Introduction:

* "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,..”. I think instead of “transversal” you mean “cross-sectional”

* Methods:

* The paper really does not make sense as it is written, because you jump into Results before explaining the methods. This doesn’t work for a modeling paper. Please include all the model motivation and details that are now at the very end of the paper BEFORE the results section. It is impossible to interpret the results without understanding why you chose the model you did, what all the parameters mean, if you did group level fitting or individual level, which parameters were fixed vs fit, and why were some values fixed and what source did you use to estimate them?

* I did not understand why it was necessary to include non-infectious virus in the model (I don’t think you are fitting to data that breaks down virus this way), and I did not understand the sentence “Since the ratio of titers to RNA copies represent an underestimate of the total infectious virus and that this ratio ranged from 10-6 to 10-4 (Fig. 2), we fixed μ to 10-3 "

* Results:

* "In those for which viral culture could obtained, the ratio of TCID50 to the number of RNA copies ranged between 10-4 to 10-6 (Fig. 2B).” What is TCID50 ? This term was never defined

* Figure 1: This figure is pretty useless - it just looks like a pile of points all on top of each other and you can’t follow one animal over time. I think the authors should be able to come up with a better way to show this data. For example, different color lines/points for each animal. And maybe also showing the group mean/median over time.

* Figure 2B: don’t understand what is being shown here. It seems like the same measure (viral load) is being used in both the x and y axis. And again, what is TCID50?

* Figure 3: This figure caption is short and non-informative. It should describe the model and what each of the parameters mean.

* Figure 6: It is obvious from looking at these plots that there is there is no correlation here. I am not sure it is necessary to report all these statistics

* Table 1: What does RSE% mean? Abbreviations should be explained. Table should have another column stating the description of each parameter. The meanings of most of these parameters aren’t explained anywhere - not in the text, or figure caption.

* Re inoculum size - why is it predicted to take longer for viral clearance with lower innoculum? Would be nice to offer an intuitive explanation, since I found it surprising! Is this because of your assumption that target cell limitation is only mechanism of control? If so, it might be a good experimental test of target cell limited model - try smaller innoculums experimentally and see what happens!

* Discussion: Sentences with grammar problems

* “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”

* "Such estimate is unknown to our knowledge in cynomolgus macaques”

* And I don’t understand what this means: "…we nonetheless note that it leads to coherent parameter estimates verifying the condition of >R0 (i.e., the burst size of the infectious virus is larger than the number of secondary cell infection) "

Final note: When submitting a paper for review, the authors should include figures in the main text near where they are first cited, with the figure captions directly below the images. It is very inconsiderate to reviewers to expect them to jump back and forth between different sections of the paper to try to interpret figures. Also, some of the images are bad quality, which makes it very difficult to interpret the results

Reviewer #4: The authors use a simple within-host model to describe the kinetics of SARS-CoV-2. Although the model itself is not new, the authors make use of both in vivo and in vitro data to estimate important parameters (e.g. R0 and the lifespan of infected cells) and distinguish likely infectious virus from non-infectious virus. Overall the manuscript is clear and well-written though lacking in detail in a few places.

Major points:

1. The authors assume HCQ reduces viral production, but are there other possible mechanisms (e.g. reducing the transmission rate, Supplementary File 1, Fig S1), and if so, does accounting for these support the assumption of no effect of HCQ on viral dynamics?

2. Results page 10: The authors focus on drugs with 90 and 99% efficacy and identify 90% as a target efficacy for future drug development that would prevent infection. But were lower efficacies investigated? For example, was 90% the lowest efficacy at which infection was prevented? Even if so, drugs with lower efficacies may still limit infection to a significant enough degree that it would be beneficial to include them as potential targets.

3. Discussion: It would informative to discuss how the parameter estimates compare to those of other within-host models for SARS-CoV-2 (at least in human studies, there have been a wide range of reported R0 and infected cell lifespan estimates).

4. Methods: To improve reproducibility, more information on the fitting procedure is needed, for example, on the backward selection procedure, and the implementation in Monolix (e.g. initial estimates). These could be included as supplementary material. I would encourage the authors to share their code for complete reproducibility.

5. Given that infection is mild here it seems perhaps unlikely that infection is limited by depletion of all susceptible cells which - despite their assumed low abundance — would likely still involve considerable damage and inflammation. It seems at least equally plausible that the infection is limited by innate immune mechanisms or NAbs. Could you justify the choice of a target cell limited model a priori? Does the choice of model influence the key conclusions? Some discussion of this issue would strengthen the paper.

Minor points:

1. p7 - p x T0 is referred to as total viral production — this seems a little misleading as it’s never the case that all T0 targets are simultaneously infected. p15 - you say “First, the term p×T0 is the only identifiable quantity in our model, “… this is a typo - I presume you’re saying that only this combination is identifiable, not p and T0 individually.

2. Results of prophylaxis treatment, page 10: The authors investigate different levels of viral inoculum, but it is not clear how these levels, quoted in PFU, were translated to V(0), quoted in copies (e.g. Table 1), to perform the simulations?

3. Model comparisons: In model comparison tables (e.g. Supplementary FiIe 3, Table S1) it would be easier for the viewer to distinguish between model support by quoting the difference in BIC between each model and the reference model (rather than the absolute BIC value, which is meaningless)

4. Fig 4: Some macaques, nasal virus load persists for some time and there is even a suggestion of a second peak (e.g. MF13, MF17, MF29). Can the authors speculate as to why this may be (in some cases these later measurements appear to be above the threshold for obtaining infectious virus)?

5. Fig 7: It would be useful to the reader to add a horizontal line to each panel marking the threshold below which infectious virus could not be obtained in the in vitro experiments.

6. Fig 7 bottom row panels (tracheal viral loads): As I understand things, the solid curves represent dynamics for a drug with 0% efficacy (i.e. no drug) and should be the same regardless of mode of action. In other words, for each viral inoculum panel in Fig 7, the solid curve should be the same as that in the corresponding viral inoculum panel in Figures S1 and S2 (Supplementary File 4). However, this does not seem to be the case (e.g. the solid curves in Fig 7 appear to peak later than those in Figs S1 and S2)?

More generally, the authors should show somewhere in the main text or supplementary material how each mode of action is incorporated into the model equations.

7. General comment: Did the authors find any differences in dynamics between male and female macaques? As differences have been reported with respect to human infection, it would be interesting to know if any differences are apparent in the macaque data.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: None

Reviewer #3: No: The authors just say "data available upon request" but they should provide the data in spreadsheets with the SI or on a public repository. Ideally same for Monolix code

Reviewer #4: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Reviewer #4: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see

21 Jan 2021

Submitted filename: 19012021_response_to_reviewers.docx

Click here for additional data file.

11 Feb 2021

Dear Mr. Gonçalves,

We are pleased to inform you that your manuscript 'SARS-CoV-2 viral dynamics in non-human primates' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Roland R Regoes

Associate Editor

PLOS Computational Biology

Rob De Boer

Deputy Editor

PLOS Computational Biology

***********************************************************

13 Mar 2021

PCOMPBIOL-D-20-01686R1

SARS-CoV-2 viral dynamics in non-human primates

Dear Dr Gonçalves,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Alice Ellingham

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol

Table 1

Population parameter estimates of the final model described by Eq (1) to (5).

Numbers in parenthesis are the relative standard error expressed in percentage (RSE%) associated either to the fixed or the standard deviation (SD) of random effects.

Parameters (units)	Fixed effects (RSE%)	SD of random effects (RSE%)
βT (mL.copie/d)	1.8×10−3 (42)	0.3 (33)
βN (mL.copie/d)	1.2×10−3 (12)
pT (copies/cell/d)	2.2×104 (49)	1.0 (27)
pN (copies/cell/d)	4.8×104 (40)
VT,0 (copies/mL)	8.0×107 (9)	-
VN,0 (copies/mL)	1.7×106 (9)	-
δ (1/d)	1.9 (9)	0.2 (37)
c (1/d)	10 (fixed)	-
μ (unitless)	10−3 (fixed)	-
k (1/d)	3 (fixed)	-
TT,0 (cells/mL)	1.88×104 (fixed)	-
TN,0 (cells/mL)	1.98×104 (fixed)
σT	1.06 (6)	-
σN	1.19 (6)	-