## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Zika plasma viral dynamics in nonhuman primates provides insights into early infection and antiviral strategies

Edited by Bruce R. Levin, Emory University, Atlanta, GA, and approved July 3, 2017 (received for review March 9, 2017)

## Significance

In light of the recent outbreak of Zika virus (ZIKV), an understanding of the within-host viral dynamics is essential to assess persistence in vivo, transmission risk, and antiviral therapeutics. Using mathematical modeling we find that, in nonhuman primates, the viral dynamics of ZIKV are characterized by a short lifetime of infected cells (<10 h) during which enough viral particles are produced to infect ∼11 other cells. Higher disease burden is associated with changes in natural killer cell subset concentrations and with elevated expression of the cytokine MCP-1, although the mechanisms behind these associations remain unclear. In order for an antiviral treatment to effectively reduce the time to plasma viral clearance therapy should be initiated at the time of infection or given prophylactically.

## Abstract

The recent outbreak of Zika virus (ZIKV) has been associated with fetal abnormalities and neurological complications, prompting global concern. Here we present a mathematical analysis of the within-host dynamics of plasma ZIKV burden in a nonhuman primate model, allowing for characterization of the growth and clearance of ZIKV within individual macaques. We estimate that the eclipse phase for ZIKV, the time between cell infection and viral production, is most likely short (∼4 h), the median within-host basic reproductive number *R*_{0} is 10.7, the rate of viral production is rapid (>25,000 virions d^{−1}), and the lifetime of an infected cell while producing virus is ∼5 h. We also estimate that the minimum number of virions produced by an infected cell over its lifetime is ∼5,500. We assess the potential effect of an antiviral treatment that blocks viral replication, showing that the median time to undetectable plasma viral load (VL) can be reduced from ∼5 d to ∼3 d with a drug concentration ∼15 times the drug’s EC_{50} when treatment is given prophylactically starting at the time of infection. In the case of favipiravir, a polymerase inhibitor with activity against ZIKV, we predict a dose of 150 mg/kg given twice a day initiated at the time of infection can reduce the peak median VL by ∼3 logs and shorten the time to undetectable median VL by ∼2 d, whereas treatment given 2 d postinfection is mostly ineffective in accelerating plasma VL loss in macaques.

Zika virus (ZIKV) is a flavivirus with a ∼10.7-kb positive sense single-stranded RNA genome (1) that is primarily transmitted among humans via a mosquito vector. First identified in Uganda in 1947 (2), it has come to prominence due to an outbreak in 2007 in Micronesia (3) and emergence in Brazil in 2015 (4). Human infection with ZIKV is usually accompanied by a relatively mild, self-limiting fever (3, 5) but can be associated with more severe effects such as Guillain–Barré syndrome (6) and fetal microcephaly (7). Recent reports of non-mosquito-borne transmission (reviewed in ref. 8) including sexual transmission (9, 10) and maternal transmission (11, 12) raise the concern that ZIKV may spread and persist in large parts of the world. Vaccines represent the most effective way to counter the epidemic, with promising results in preclinical studies (13, 14). Despite the fact that human trials have begun, approved vaccines may not become available for years (15). Treatment with antivirals is needed, especially in some populations, such as severely infected patients or individuals for whom viremia remains detectable in the blood or in other compartments for extended periods, as has been observed in pregnancy (16) and in semen (17, 18). So far no antiviral drug has shown a significant effect against ZIKV in vivo. The WHO has highlighted a number of potential drugs that could be relevant against ZIKV (www.who.int/blueprint/priority-diseases/key-action/zika-rd-pipeline.pdf?ua=1) and a high throughput in vitro screen (19) has identified a number of potential small molecules with activity against ZIKV but that will require further testing in animal models.

With that perspective, the establishment of a relevant nonhuman primate (NHP) model (14, 20⇓–22) represents crucial progress toward a better understanding of virus pathogenesis and the development of more efficacious therapeutic strategies. Zika infection in NHPs has recapitulated many of the key clinical findings in human infection, including rapid control of acute viremia, early invasion of the central nervous system, and prolonged viral shedding and fetal pathology in pregnant animals (14, 20, 21, 23, 24). Thus, the study of ZIKV in NHP models can provide insight into viral dynamics as well as providing an essential tool for testing potential antiviral drugs and vaccines (14).

Mathematical modeling of viral dynamics has provided a quantitative understanding of viral dynamics for a number of viruses, including HIV (25, 26), influenza (27), and the flaviviruses hepatitis C virus (28), West Nile virus (WNV) (29), and dengue (30⇓–32). One of the main insights of these models was to provide estimates of parameters that are not easily measurable, such as the half-life (t_{1/2}) of productively infected cells in vivo and the basic reproductive ratio, *R*_{0}, the number of secondary cells infected by viral production from one infected cell in a population of target cells. In the context of viruses for which antivirals are available, such as HIV and hepatitis C virus, mathematical modeling has been used to optimize and evaluate in silico the impact of antiviral strategies (33⇓⇓–36).

Here we use mathematical modeling to characterize ZIKV dynamics in vivo in NHPs, to explore correlations between viral dynamic parameters and the immunologic response, and to simulate the effects of a potential antiviral agent.

## Results

### Zika Viral Kinetics.

After s.c. infection of nine rhesus macaques with 10^{6} pfu of a Thai isolate of ZIKV (*Materials and Methods*), the plasma viral load (VL; data given in *SI Appendix*, Table S1) was high by day 1 postinfection (>10^{5} RNA copies per mL in all monkeys) and peaked at day 2 (Fig. 1). The peak VL was above 10^{6} copies per mL for all monkeys and above 10^{7} copies per mL for six of the nine monkeys. Following peak viremia, there was a rapid decline in VL, and the VL became undetectable (<200 copies per mL) in eight of the nine monkeys by day 5. The VL in the remaining monkey was undetectable by day 10.

### Viral Kinetic Model and Parameter Estimates.

Using nonlinear mixed effects modeling we fit eight different models to the plasma VL data, including ones with immune effects via measured IFN-α, natural killer (NK) cell, or monocyte chemoattractant protein 1 (MCP-1) dynamics. Using model selection theory (37) we were unable to distinguish models incorporating explicit immune effects from a target cell-limited model, because all models had similar Bayesian information criterion (BIC) values (details are given in *SI Appendix*). We chose to use a simple model, a target cell-limited model with an eclipse phase as given by Eq. **1** in *Materials and Methods*, because it allowed us to test the effects of antiviral therapy. We considered fits of this model with different fixed values of *V*(0) (effective initial plasma VL concentration), *k* (rate of transition to productive infection), and *c* (rate of viral clearance). The model fit with *V*(0) = 10^{4}/mL and *k* = 6 d^{−1} was chosen based on the BIC across all values of *c* tested (*SI Appendix*, Figs. S1*A* and S2*A*). With these parameter values, the quality of model fit with *c* = 5 d^{−1}, 10 d^{−1}, or 15 d^{−1} was broadly indistinguishable (within ∼2 BIC points, *SI Appendix*, Fig. S1*A*), and we selected the model with *c* = 10 d^{−1} to analyze further. We also considered the effect of the initial target cell density, *T*(0) or *T _{0}*, selecting

*T*= 10

_{0}^{5}/mL (

*SI Appendix*and

*SI Appendix*, Fig. S1

*C*and

*D*) and ω, the SD of the random effect term, selecting ω = 0.4 (

*SI Appendix*and

*SI Appendix*, Fig. S1

*B*). With these values, the fitting process provided an estimate of the median basic reproductive ratio

*R*

_{0}of 10.7 (

*SI Appendix*, Fig. S2

*B*) with a relative SE (r.s.e) of 18% and the productively infected cell death rate, δ, was estimated as 4.5 d

^{−1}(

*SI Appendix*, Fig. S2

*C*, r.s.e. 15%). The estimated value of

*p*is dependent on

*T*, the initial target cell density, and only the product

_{0}*pT*can be reliably estimated (38), with the fitting process estimating

_{0}*pT*as 1.2 × 10

_{0}^{9}d

^{−1}⋅mL

^{−1}(

*SI Appendix*, Fig. S2

*D*). The individual parameters derived from this population fit provide estimated VLs for each monkey as shown in Fig. 1, and the predicted target cell dynamics for each individual monkey (

*SI Appendix*, Fig. S3) show that target cells are largely but not completely eliminated by infection. Estimates of the value of

*R*

_{0}for each monkey range between 8.5 and 16.5 (Fig. 1). The goodness of fit and estimated parameters for each combination of

*V*(

*0*) and

*k*are shown in

*SI Appendix*, Fig. S2, and under the model averaging approach (

*Materials and Methods*) the distribution of

*R*

_{0}in the population had a median of 10.5 with range 5.3–19.5 (5th–95th percentiles). We also estimated parameters

*R*

_{0}, δ, and

*pT*

_{0}by fitting the model to each individual monkey using a nonlinear least-squares approach (

*SI Appendix*, Fig. S4). This approach gave values of

*R*

_{0}that ranged between 5.3 and 25.9 and that have a very similar median (10.6) to that estimated by population modeling.

### Correlation Between Plasma VL and VL in Other Compartments.

There was a statistically significant positive correlation, after Bonferroni correction for multiple testing, between the area under the log_{10} VL curve (AUC) between days 0 and 14 in plasma calculated by the model and in saliva (*SI Appendix*, Fig. S5*A*, Spearman correlation, *P* = 0.05). There were no other statistically significant correlations at the *P* = 0.05 level between modeled AUC in plasma and AUC in other compartments (*SI Appendix*, Fig. S5 *B*–*E*) or between peak log_{10} VL in plasma as predicted by the model and peak log_{10} VL in other compartments (*SI Appendix*, Fig. S6).

### Immunological Markers During ZIKV Infection.

We tested for correlations between immune cell subset or cytokine concentrations during ZIKV infection and the modeled plasma VL. None of the Spearman correlations were statistically significant after correction for multiple testing, but we considered those cell subsets and cytokines with the strongest correlations with VL to look for suggestions of an immune response to ZIKV infection that we could further examine via model fitting.

Of the cell subsets, NK CD16^{+} CD69^{+} cells showed the strongest positive correlation (*r* = 0.83 between subset concentration at day 0 and peak log_{10} VL, *SI Appendix*, Fig. S7 *A* and *B*), and NK CD16^{−} CD56^{+} cells (*SI Appendix*, Table S4) show the strongest negative correlation (*r* = −0.83 between subset concentration at day 2 and log_{10} VL AUC between days 0 and 14, *SI Appendix*, Fig. S7 *C* and *D*). Biologically we might expect NK cells to aid in clearance of infected cells, but model fitting incorporating the NK CD16^{−} CD56^{+} subset did not provide strong support for or against these cells controlling ZIKV VL (Model G, *SI Appendix*, Table S5).

Of the cytokines, MCP-1 (*SI Appendix*, Table S2) showed the strongest positive correlation (*r* = 0.88 between MCP-1 concentration at day 1 and estimated infected cell half-life, *SI Appendix*, Fig. S7 *E* and *F*). A mechanism by which MCP-1 recruits target cells might explain the observed positive correlation, but model fitting (Model H, *SI Appendix*, Table S5) does not provide strong evidence for or against this hypothesis. Additionally, we hypothesized that IFN-α (*SI Appendix*, Fig. S8 and Table S3) might be involved in immune control of VL, despite the absence of a correlation with VL (*SI Appendix*, Fig. S9). Models incorporating an IFN-α response (Models D–F, *SI Appendix*, Table S5) were all indistinguishable from our basic model, Eq. **1**, using model selection theory.

### Evaluating the Effect of an Antiviral Treatment.

For an antiviral that reduces the viral production rate, the effect on the VL dynamics can be assessed using a target cell-limited model if the EC_{50} against ZIKV and the free drug concentration are known. Because we do not know how lowered VL due to antiviral therapy will influence immune responses, models incorporating the measured immune markers in these data are not suitable for studying the effect of an antiviral. The model, Eq. **1**, predicts that when the drug concentration *C* is low relative to its EC_{50}, *R*_{0} declines almost linearly with increasing drug concentration (Fig. 2*A*). However, once the drug concentration is well above the EC_{50} further increases in *C* have little effect. To reduce the median value of *R*_{0} to below 1, meaning no establishment of infection, a constant free drug concentration 9.7 times the EC_{50} (i.e., 9.7 EC_{50}) is required. Similarly, the peak VL decreases monotonically with increasing *C* (Fig. 2*B*).

Despite the monotonic reduction in *R*_{0} with increasing *C*, the time until VL becomes undetectable initially increases with *C* before declining at higher ratios of *C* to EC_{50} (Fig. 2*C*). This initial increase at low drug concentrations is due to the effect of the drug in reducing the number of viral particles produced and therefore the rate at which target cells become productively infected. Consequently, the time until the target cell population is substantially depleted is increased in the model, leading to a prolonged viral infection, albeit at lower levels than without treatment. The median time until the VL becomes undetectable is maximal, at 16.1 d, when *C* is 7.6 EC_{50} and becomes lower than the time with no drug (4.7 d) when *C* is 12.8 EC_{50}, reaching 3 d when *C* is ∼15 EC_{50}. The VL AUC exhibits similar behavior as constant *C* is increased (Fig. 2*D*). Whether infection will be prolonged at low drug concentrations in reality is unclear because an adaptive immune response most likely will arise and limit the infection (23).

To model the effects of an antiviral drug with activity against ZIKV under more realistic conditions of changing drug concentration we examined favipiravir, an antiviral drug, approved for human use, with activity against many RNA viruses (39) and with known pharmacokinetics (PK) in monkeys (40). We used a previously developed PK model that shows differences between male and female monkeys (*SI Appendix*), assumed a 150 mg/kg twice daily dosing regime initiated at the time of infection, and chose PK (for male monkeys) and viral kinetic parameters randomly from their estimated distributions to simulate VL dynamics for 5,000 monkeys. The peak median VL was substantially reduced by ∼3 logs for male monkeys and the time to first undetectable median VL was shortened by ∼2 d (Fig. 3*A*, red). Initiating treatment 2 d after infection had essentially no effect on peak log_{10} VL and only shortened the time to undetectable VL by ∼0.5 d (Fig. 3*A*, blue). Assessing the reduction in VL burden by the AUC calculated between 0 and 14 d postinfection (dpi) (Fig. 3*B*) shows a reduction in median log_{10} VL AUC from 14.5 d copies per mL with no drug (Fig. 3*B*, black) to 4.7 d copies per mL (with large variation between simulations) when treatment is initiated at the time of infection (Fig. 3*B*, red). When treatment is initiated 2 d after infection, the effect on median log_{10} VL AUC in less substantial (Fig. 3*B*, blue). The equivalent analysis for female monkeys shows a less substantial effect on VL (*SI Appendix*, Fig. S10).

## Discussion

This study presents an analysis of the within host dynamics of ZIKV in NHPs. After s.c. infection with 10^{6} pfu of a Thai isolate, ZIKV is found at high concentrations in the plasma at day 1 and peaks at day 2 in these animals, followed by a rapid decline. These rapid viral kinetics resemble those observed in influenza A infections in humans, where VL tends to peak 2–3 dpi and declines to undetectable levels by 6–7 dpi in many individuals (27). This similarity motivated the use of a target cell-limited viral infection model developed by Baccam et al. (27) for influenza to fit the ZIKV plasma VL data, and the use of similar methods of analysis as in ref. 41.

The early peak VL (2 dpi) makes it is unlikely that an adaptive immune response has been mounted and is responsible for the postpeak viral decline. However, innate immune responses may play a role, as has been suggested for dengue (31), and cytokines such as type I IFN could place some uninfected cells into an antiviral state, protecting them from infection. Fitting models that include IFN-α, using its measured dynamics, did not improve the fit to the data sufficiently to justify including these innate immune responses in the model, consistent with the observation that IFN-α is not correlated with VL (*SI Appendix*, Fig. S9) and that ZIKV disrupts IFN signaling via STAT2 degradation (42⇓–44). Similarly, the modeling did not strongly support or reject including an NK cell-mediated innate immune response. It may be that these data do not contain enough information to allow inference about the presence of an innate immune response, and additional VL measurements at earlier time points or after infection with lower doses may allow more detailed assessment of the contribution of the innate immune system to control of ZIKV in these animals. It may be that inclusion of immune data could reduce the unexplained individual variability of the parameters, but these data do not contain enough information to allow us to estimate this. Given the positive relationship between the expression of MCP-1, a cytokine implicated in recruitment of cells to the site of infection (45), and the predicted infected cell half-life (*SI Appendix*, Fig. S7*F*), we also considered a viral dynamics model in which the target cell population was increased in an MCP-1–dependent manner (model H, *SI Appendix*, Table S5). In fitting this model to the VL data from each monkey, the model BIC suggested that the available data cannot support or reject the hypothesis that MCP-1 affects the VL via this mechanism.

We highlight the point that these data do not contain enough information to allow for either acceptance or rejection of the hypothesis that an immune response is responsible for control of plasma ZIKV in these animals. Our modeling approach is to select the simplest model from those that provide equivalently good fits to the data, and as such we have performed analysis without explicit inclusion of an immune effect.

The initial viral infection is characterized by a basic reproductive ratio (*R*_{0}) of 10.7, suggesting that early after infection each infected cell on average infects ∼11 other cells. To attain rapid spread, we estimate the mean eclipse phase for ZIKV (i.e., the time between infection and viral production) is ∼4 h. This estimated eclipse phase length is shorter than that for WNV in mice (29) and that estimated for influenza A in humans (27). In vitro ZIKV viral production from infected human skin cells was detected at 6 h, the earliest time examined (46). We also estimate that once cells begin producing virus they live an average of 5.3 h to yield a total lifespan, which includes the eclipse phase, of slightly less than 10 h.

Using the model, Eq. **1**, we can only estimate the product *pT*_{0} (38), where *T*_{0} is the initial target cell density and *p* is the rate of viral production from a productively infected cell. Because the target cells for ZIKV are not currently fully characterized, we have no a priori estimate of *T*_{0}. However, we can constrain *p* because the amount of infectious virus made over an infected cell’s productive lifetime of 5.3 h needs to be greater than *R*_{0}. For a French Polynesian isolate of ZIKV it has been estimated that 500–1,000 virions correspond to 1 pfu (20). Thus, with *R*_{0} ∼11, an infected cell needs to produce at a minimum 500 × 11= 5,500 virions over 5.3 h. This implies a daily production rate *p* ≥ 25,000 virions d^{−1}. From our model fitting, we estimated that *pT*_{0} = 1.2 × 10^{9} d^{−1}⋅mL^{−1}, so that with *p* ≥ 25,000 virions d^{−1} we estimate *T*_{0} ≤ 4.8 × 10^{4}/mL. Our finding of a low initial target cell density combined with the broad range of cell types that ZIKV has been shown to be able to productively infect in vitro (44, 46⇓⇓⇓–50) might suggest a role of an early innate response, not captured in these data, placing potential target cells into an antiviral state, thus reducing the effective target cell population.

The early peak viremia and rapid viral decay in this NHP model, along with the limited longitudinal sampling, means that the chosen value of *V*(*0*) and the estimate of *R*_{0} should be taken cautiously. Plasma samples taken during the first few hours postinfection would allow for a more accurate assessment of the likely value of *V*(*0*) and would give insight into how rapidly a s.c. infection is trafficked to the blood. Also, studies with different inoculum sizes would help in determining a relationship between dose and *V*(*0*). The number of data points available also limits our ability to estimate the distribution of parameters found in the population, because these data are unable to support fitting the variance of the estimated parameter distributions. Additional animals as well as more frequent sampling would improve the estimates of these distributions. The lack of VL AUC correlation between the plasma and semen compartments (*SI Appendix*, Fig. S5*D*) might be due to the small number of data points available in this study, or it might indicate that different mechanisms affect the VL in the semen, which would restrict our ability to assess sexual transmission risk from plasma samples. Overall, the persistence of virus in saliva and semen (17) and the prolonged presence in the plasma of pregnant macaques and humans (16, 20) as well as in a child born to ZIKV-infected parents (51) remains to be explained.

In VL dynamics from a study with different ZIKV strains and with lower inoculum doses (23) it was seen that peak VL tended to be later than in these data, and a model incorporating immune control from IFN-α (reducing viral production) and antibodies (reducing viral infectivity) was able to improve the fit of the model. The estimated value of δ was found to be very similar to this study, whereas the estimated value of *R*_{0} was much lower (∼3). This discrepancy between the estimated values of *R*_{0} can be attributed to a later peak with lower inoculum doses as well as a higher estimate of *V*(*0*) with equivalent inoculum doses, highlighting the need for early samples to assess the true initial plasma viral concentration.

Because antiviral drugs against ZIKV are being tested in rhesus macaques (52) we investigated the effect of an antiviral drug at concentration *C* that reduces the production of virus by a factor of 1 – *C*/(EC_{50} + *C*). Increasing the drug concentration reduces *R*_{0}, but we estimate that a drug concentration of at least 12.8 EC_{50} is required to reduce the time to undetectable VL (Fig. 2), where these results are based on our target cell limitation model assumptions. This approach can allow for rapid screening of potential repurposed drug candidates: if a drug is unable to reach this concentration with a safe dose, it is unlikely to have a substantial positive impact of the VL dynamics, although both differences between EC_{50} in vitro and in vivo and differences between the VL dynamics in the NHP model and in humans need to be accounted for. We identified favipiravir, a small-molecule antiviral, as capable of reducing the median viral burden in NHP infections in simulations and found that substantial reduction can be achieved for the majority, but not all, of simulations when treatment is started at the time of infection. The PK of favipiravir is highly variable in NHPs (40), resulting in the broad distribution of VL AUCs seen in our simulations (Fig. 3*B*), although the sustained infections seen in our simulations when the PK is unfavorable, giving low drug concentrations, would likely be mitigated by later immune responses. However, if an antiviral strategy is to have any substantial effect on peak VL, time to undetectable viremia or VL AUC it needs to be initiated before peak viremia. Because peak plasma VL occurs very early in this NHP model, this effectively means that treatment needs to be given prophylactically. In humans where lower viremia levels are often observed (e.g., refs. 18, 53, and 54), as well as when lower inoculum doses are used in NHPs (14) the viral dynamics will be altered and the peak of plasma viremia will likely be delayed, affecting the assessment of antiviral impact and potentially allowing for later antiviral treatment initiation.

The mathematical modeling approach in this study describes the characteristics of the viral dynamics of ZIKV as observed in the plasma of an NHP model and provides a tool by which to assess, in silico, the effect of candidate antivirals.

## Materials and Methods

### ZIKV-Infected Rhesus Macaques.

Ten Indian-origin rhesus macaques (five male and five female) were infected s.c. with 10^{6} pfu of a Thai isolate of ZIKV as described in ref. 21. Samples were collected from blood (plasma), urine, saliva, cerebrospinal fluid, semen, and vaginal secretions for 4 wk postinfection. Samples were taken regularly, and viral RNA levels in samples were assessed using real-time PCR, with a limit of detection of 200 copies per mL. Additionally, peripheral blood mononuclear cells (PBMC) were analyzed from the same blood samples. PBMC samples were stained with fluorescent antibody mixes and flow cytometry used to quantify the concentration of T-cell, B-cell, and NK subsets. The set of immune cells measured includes total lymphocytes, eight T-cell subsets, three B-cell subsets, and four NK cell subsets, and both total concentration and activated (as indicated by expression of CD69) concentration were measured (21). Cytokine expression in plasma samples was also measured by multiplex bead immunoassays at days 0, 1, 3 and 7. All animals were housed at Bioqual Inc. or the Wisconsin National Primate Research Center (WNPRC). The described studies were approved by the Institutional Animal Care and Use Committee of Bioqual Inc. or WNPRC, respectively.

### Viral Dynamic Model.

Plasma viral dynamics in nine of the 10 animals were characterized using a standard target cell-limited model with an eclipse phase originally developed to analyze influenza infection kinetics (27). One animal in which there were only two plasma VLs above the detection limit of the assay was excluded from this analysis. The model with an eclipse phase was used in preference to a model without an eclipse phase (*SI Appendix*). The model is given by the following system of ordinary differential equations:

In this model, target cells (*T*) are infected by free virus (*V*) at a rate proportional to both virus and target cell concentrations, as given by the law of mass action term *βVT*. After infection, the cell enters an eclipse phase, denoted *I*_{1}, where it is infected but does not produce virus. It then transitions to a productively infected state (*I*_{2}) at rate *k*, that is, 1/*k* is the average duration of the eclipse phase. Productively infected cells are cleared at per capita rate *δ*, incorporating the observed cytopathic effects of ZIKV (55, 56), and produce virus at rate *p* per cell. Free virus is cleared at rate *c* per virion. Because of the short-term nature of the infection as measured in blood with a peak VL at day 2 in all nine animals, target cell replenishment was ignored in the model.

The model predicts an initial decline of virus until the end of the eclipse phase of the first infected cells, followed by an exponential increase until target cells become limiting. At the end of this expansion phase, the target cells are largely but not completely depleted (*SI Appendix*, Fig. S3), the number of new infections is low, and the loss of infected cells largely outnumbers the number of new cell infections. Consequently, the number of infected cells, and hence the total viral production, declines. Mathematical analysis shows that this decline is exponential with a rate proportional to the infected cell death rate, δ, when it is much smaller than *c* and *k* (57). In this model all parameters were assumed to be constant, but in principle they could vary over time. Thus, for instance, an adaptive immune response, if present, could be associated with an increasing value of δ over time. One can derive from this model the basic reproductive number, *R*_{0}, representing the number of secondary infected cells that would be created by a single infected cell introduced into a population of *T*_{0} target cells, given by *R*_{0} *= βpT*_{0} */*δ*c*. The model was reparametrized as a function of *R*_{0}, which we constrained to be larger than 1 to ensure productive infection. More details on viral dynamics models during acute infection can be found in refs. 57⇓–59. We additionally tested a number of other potential models, including a model without an eclipse phase and models incorporating innate immune responses (*SI Appendix*). Model selection was based on BIC, with more complex models only being accepted if an improvement in BIC of >2 points was observed (60).

### Fixed Parameters.

Because not all model parameters can be estimated using only VL data (38) a number of parameters had to be fixed (details are given in *SI Appendix*). First the number of nonproductively and productively infected cells at time of infection is assumed to be 0 [*I*_{1}(*0*) *= I*_{2}(*0*) *= 0*], and the number of target cells at the time of infection, *T*(*0*), is assumed to be 10^{5} per mL. As shown in ref. 38, for a standard viral dynamic model the viral production rate *p* and the initial target cell density *T*(*0*) are dependent, and with measurements of VL only one is able to be estimated. Here the assumed value of *T*(*0*) was chosen to give biologically reasonable estimates of *p* (*SI Appendix*). Note that although the value of *T*(*0*) was chosen to be plausible it does not affect our main findings regarding the basic reproductive number *R*_{0}, the half-life of infected cells, and the effect of antiviral treatment. Additionally, we do not assume that all target cells and infection events are in blood, but we model the effective impact of target cell infection on the plasma VL as has been done, for example, in models of hepatitis viruses (28), whose target cells are mainly or exclusively in the liver. After initial model fitting to determine the effect of the rate of clearance of free virus, *c*, on the model fit (details are given in *SI Appendix*) we set *c* to a fixed value of 10 d^{−1}.

### Estimated Parameters.

The macaques were inoculated s.c. with 10^{6} pfu ZIKV. Assuming that 1 pfu represents between 500 and 1,000 copies of viral RNA (20) and that a typical macaque plasma volume is about 300 mL, the inoculum would at most generate between 1.7 × 10^{6} and 3.3 × 10^{6} viral RNA copies per mL of plasma. Because it is likely that much of the early infection occurs in the skin and that only a small proportion of the inoculum is transported to the blood immediately after infection, we used values of *V*(*0*) between 10^{3} and 10^{5} viral RNA copies per mL in our model fitting. Here *V*(*0*) should be interpreted as an effective initial viral concentration.

The eclipse phase duration for ZIKV is unknown. We let *k* be 4, 6, or 8 d^{−1} (representing a mean eclipse phase of 6, 4, or 3 h). We restricted our range of possible eclipse phase lengths to short times because the Zika plasma VL can be as high as 10^{7} RNA copies per mL by 1 d postinfection, suggesting rapid viral replication and hence a short eclipse phase. Further, in vitro experiments have observed ZIKV viral production 6 h postinfection (46).

For each combination of *V*(0) and *k*, the distribution of other parameters across the population of monkeys was estimated using a nonlinear mixed effects approach (61). Under this approach, each parameter is assumed to be log-normally distributed in the population and can be expressed as *θ* is the median value in the population and *η*_{i} is the individual random effect, which is normally distributed as *N*(0, ω^{2}), accounting for variability between individuals in the population. Because the variance of this random-effect term could not be estimated with the small number of animals in this study, we set it to a fixed value for fitted parameters *R*_{0}, δ, and *p* (SD ω = 0.4 for each parameter), corresponding to a coefficient of variation just above 40% for all fitted parameters. We also tested the effect of varying the SD of the random effect (details are given in *SI Appendix*) and found that the models are equally supported (within 2 BIC points) for ω = 0.2, 0.3, and 0.4 (*SI Appendix*, Fig. S1*B*) and chose ω = 0.4 to allow the most parameter variation. Estimates for *R*_{0}, δ, and *p* (given a fixed value of *T*_{0}) were then estimated using the stochastic approximation of the expectation maximization algorithm implemented in Monolix (monolix.lixoft.com/), with data below the limit of detection accounted for as described in ref. 62. The combination of *V*(0) and *k* that provides the lowest BIC was carried forward for further analyses, and the individual parameters for each monkey were obtained using empirical Bayes estimates.

Using these individual parameters, the individual VL profiles were reconstructed, allowing us to calculate the area under the VL curve (AUC) for each animal. In this study, the plasma AUC is calculated as the area between the modeled log_{10} VL curve and the log_{10} detection limit of the experimental assay, calculated using the trapz function of NumPy in Python (63). The AUCs for other compartments were calculated in a similar fashion, using observed rather than modeled VLs.

To account for the fact that a different value of *R*_{0} could have been obtained with different pairs of values for *V*(*0*) and *k* we also assessed the distribution of *R*_{0} in the population by using a “model averaging” approach (64). We assume that *R*_{0} is distributed as a weighted sum of Gaussian distributions, each corresponding to one combination of *V*(*0*) and *k* considered, with weight given by *e*^{-BIC/2}. The mean of the Gaussian is given by the estimated population *R*_{0} with that *V*(*0*) and *k*, and the SD is the SE of that distribution. We then sampled 10^{6} values of *R*_{0} from this weighted sum of Gaussians and reported the median and 5th–95th percentiles.

Additionally, a simple individual fitting approach was used to confirm that conclusions similar to those in the mixed effects model with fixed variances were reached. For each monkey, the free model parameters were fitted to the observed VL via a nonlinear least squares approach using a differential evolution algorithm implemented in the DEoptim package in R (65). The parameters estimated by this approach across all nine monkeys were found to be in the same range as the distributions estimated from the mixed effects approach.

### ZIKV VL in Anatomic Compartments.

The correlations between ZIKV dynamics in the blood and the other tissue compartments were explored by comparing, for each NHP, the modeled AUC in the blood with the observed AUC in each other compartment, and the modeled peak VL in the blood with the observed peak VL in each other compartment. The correlations were assessed via a Spearman correlation coefficient (adjusted for multiple testing using the Bonferroni correction) and significance was assessed at the *P* = 0.05 level.

### Immune Markers.

Concentrations of a number of lymphocyte subsets and cytokines were measured in each of the monkeys during infection, and we refer to these as immune markers (see ref. 21 for details). For each monkey, we calculated Spearman correlations between characteristics of the immune markers and characteristics of the VL dynamics, as described in *SI Appendix*. We adjusted for multiple testing using the Bonferroni correction. To identify immune markers to be incorporated into our viral dynamic models, we considered those immune markers that demonstrated the strongest correlations with VL and seemed biologically relevant. We used NK CD16^{−} CD56^{+} cells and the cytokine MCP-1 in models G and H (*SI Appendix*). Additionally, we considered correlations between the concentration of IFN-α and the VL in each monkey, because IFN-α may contribute to immune control of ZIKV. IFN-α dynamics were incorporated in models D–F (*SI Appendix*).

### Effect of an Antiviral Treatment on Plasma VL.

Using repeated simulations with parameters selected at random from the fitted estimates of the parameter distributions, we explored the effects of an antiviral drug on the viral dynamics for different ratios of an assumed constant drug concentration (*C*) to the drug EC_{50,} where EC_{50} is the concentration providing a 50% viral growth inhibition in vivo. We assumed the drug reduced the viral production rate, *p*, by a factor (1− *ε*) where *ε = C*/(EC_{50} + *C*) is the drug antiviral effectiveness. Then we calculated the impact on the basic reproductive number *R*_{0}, peak VL, time to undetectable VL, and VL AUC for different (constant) values of *C*/EC_{50}.

### Effect of Favipiravir on ZIKV Plasma VL.

Our group identified favipiravir (T-705), a broad-spectrum RNA polymerase inhibitor (66), as a potential ZIKV inhibitor. The EC_{50} values of favipiravir against ZIKV strains was assessed as detailed in *SI Appendix* and provided values between 2.7 µg/mL and 6.6 µg/mL. Because the PK of favipiravir in cynomolgus macaques have been described in detail (40), we aimed to evaluate the effects of the drug on the VL dynamics by assuming the PK in rhesus macaques was similar to that described in Mauritian cynomolgus macaques (*SI Appendix*) and performing simulations as follows.

We simulated a dosing regimen including two loading doses of 250 mg/kg at hours 0 and 12, followed by 12 hourly maintenance doses of 150 mg/kg for 14 d. This dosing regimen was chosen because it has been shown to give no serious abnormalities in NHP, while maintaining a trough concentration at day 7 greater than the EC_{50} against ZIKV (40). Different treatment initiation times were tested, with dosing starting either at the time of infection or 2 dpi (median time of viral peak). For each scenario, 5,000 PK and viral kinetic parameters were drawn from their estimated distributions, giving simulated trajectories of response. Because the PK in male and female animals differed, this process was done for males and females separately. To summarize the simulations, the median of the VLs at each time point and the distribution of the AUC of log_{10} VL between 0 and 14 dpi over the 5,000 trajectories were reported for males and females separately.

## Acknowledgments

We thank France Mentré for advice on the model averaging approach. Portions of this work were done under the auspices of the US Department of Energy under Contract DE-AC52-06NA25396. This work was also supported by NIH Grants R01-AI078881, R01-OD0110955, and R01-AI028433 (to A.S.P.) and the European Union’s Horizon 2020 research and innovation program under Grants 666092 (to J.G.) and 734548 (to X.d.L.).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: asp{at}lanl.gov.

Author contributions: K.B., J.G., J.B.W., and A.S.P. designed research; K.B., V.M., X.d.L., S.-Y.L., C.E.O., J.B.W., and A.S.P. performed research; K.B., J.G., V.M., and A.S.P. analyzed data; and K.B., J.G., and A.S.P. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1704011114/-/DCSupplemental.

## References

- ↵.
- Lei J, et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Grischott F,
- Puhan M,
- Hatz C,
- Schlagenhauf P

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Abbink P, et al.

- ↵.
- Ferguson NM, et al.

- ↵.
- Suy A, et al.

- ↵.
- Mansuy JM, et al.

- ↵
- ↵
- ↵
- ↵.
- Osuna CE, et al.

- ↵
- ↵
- ↵
- ↵.
- Ribeiro RM, et al.

- ↵
- ↵.
- Baccam P,
- Beauchemin C,
- Macken CA,
- Hayden FG,
- Perelson AS

- ↵.
- Neumann AU, et al.

- ↵.
- Banerjee S,
- Guedj J,
- Ribeiro RM,
- Moses M,
- Perelson AS

- ↵.
- Clapham HE,
- Tricou V,
- Van Vinh Chau N,
- Simmons CP,
- Ferguson NM

- ↵.
- Ben-Shachar R,
- Koelle K

- ↵.
- Clapham HE, et al.

- ↵
- ↵
- ↵
- ↵
- ↵.
- Burnham KP,
- Anderson DR

- ↵
- ↵
- ↵.
- Madelain V, et al.

- ↵.
- Hadjichrysanthou C, et al.

- ↵.
- Kumar A, et al.

- ↵
- ↵
- ↵
- ↵.
- Hamel R, et al.

- ↵.
- Retallack H, et al.

- ↵
- ↵
- ↵
- ↵.
- Oliveira DBL, et al.

- ↵.
- Lim S-Y, et al.

- ↵
- ↵
- ↵
- ↵.
- Barr KL,
- Anderson BD,
- Prakoso D,
- Long MT

- ↵
- ↵
- ↵
- ↵
- ↵.
- Pinheiro J,
- Bates DM

- ↵
- ↵.
- Ascher D,
- Dubois PF,
- Hinsen K,
- Hugunin J,
- Oliphant T

- ↵
- ↵
- ↵.
- Sleeman K, et al.

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Biological Sciences
- Microbiology

- Physical Sciences
- Biophysics and Computational Biology