Modeling the Antibody Mediated Prevention trials

The antibody mediated prevention (AMP) trials are simultaneous global phase 2b randomized trials HVTN 703/HPTN 081 and HVTN 704/HPTN 085 (ClinicalTrials.gov #NCT02568215 and ClinicalTrials.gov #NCT02716675, respectively). To simulate the AMP trials, we adopted a slightly simplified study design as shown in Fig 1A. Here, each trial contains 3 equally populated arms: placebo, 10 mg/kg, and 30 mg/kg VRC01. Infusion of VRC01 occurs every 8 weeks for 10 total doses and HIV testing occurs every 4 weeks. Our modeling approach involved developing viral dynamics models, incorporating VRC01 PK/PD, simulating trials, and development of a clinical trial regression model (CTRM) capable of distinguishing reduced potency and insufficient breadth as causes of breakthrough infection (Fig 1B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. The model implementation of the Antibody Mediated Prevention (AMP) trials.

A) Two parallel trials in African women and North American MSM and TG individuals. Each trial contains 3 equally populated arms, placebo (gray), 10 mg/kg infusions (teal), and 30 mg/kg infusions (tan). Infusions occur every 8 weeks, for a total of 80 weeks and HIV testing occurs every 4 weeks. We refer to an 8 week interval as a ‘dosing interval’. B) The model workflow outlines the results.

https://doi.org/10.1371/journal.pcbi.1007626.g001

A deterministic mathematical model of primary HIV infection recapitulates natural primary infection viral loads and semi-quantitative host-cell dynamics

We designed and validated a deterministic mathematical model of natural HIV primary infection (Fig 2) [22]. Model equations are shown in Eq 5. In this model, susceptible cells S are naturally maintained at homeostasis with birth rate α_S and death rate δ_S. They become infected with rate β upon exposure to free virus V_f. Infected cells I_s,p have a state s ∈ [A, L] representing active or latent, and a phenotype p ∈ [U, P] representing unproductively or productively infectious provirus. The death rate of infected cells depends on their state, with viral cytopathic effects killing active cells rapidly and latent cells having a long (44 month) half-life. Only active, productively infectious cells produce virus (with rate π). An adaptive immune compartment E is recruited at rate ω depending on the total number of infected cells, but this recruitment is limited by a saturation constant E₅₀. Adaptive immune cells kill actively infected cells with rate κ and immunity wanes with rate δ_E. Virus is naturally cleared with rate γ. When VRC01 is present, free virus is reduced based on the concentration and neutralization potency of VRC01 against the virus. We focus on natural infection initially and include the model for VRC01 neutralization in the following sections.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. A mechanistic mathematical model of HIV primary infection including VRC01.

A) Model schematic: susceptible cells S become infected I_s,p. If the cell is activated s = A and the provirus is productive p = P, that cell produces virus with rate π. An adaptive immune compartment E is recruited to remove infected cells with rate ω, but saturates when E > E₅₀. The presence of VRC01 reduces free virus V_f. B) Best-fit models of natural HIV dynamics (without VRC01) using data from the RV 217 trial [23] demonstrates satisfactory fit to heterogeneous viral-load kinetics. C) Density plots of distributions of 6 estimated parameters (lines-medians, c-cells, d-days). 8 parameter values were set constant based on past studies; see Table 2 for constant values and initial guesses for estimated parameters. D) Calculated values of the basic reproductive number , the average number of new infected cells generated due to a single infected cell in a fully susceptible population. E) A representative model simulation using mean values from density plots in C recapitulates features of host-cell dynamics during HIV primary infection. Susceptible cells decrease by approximately 200 cells/μL and reconstitute by viral set-point. Latently infected cells are generated within the first week of infection and contain plausible levels of total HIV DNA (unproductively infectious latent provirus I_{s = L, p = U}) and replication competent provirus (productively infectious latent provirus I_s=L,p=P). Adaptive immune cells E become significant around the time of peak viral load (approximately 10 days after first positive) and adjust viral load set-point.

https://doi.org/10.1371/journal.pcbi.1007626.g002

We parameterized the natural viral dynamics model using human primary infection viral loads from the RV 217 study [23]. The model achieved satisfactory fit (matching peaks, times-to-peak, set-points) to heterogeneous viral load trajectories across all study participants (see Fig 2B). The complete model is not identifiable, so certain parameters were fixed at values based on experimental measurements and past modeling. In particular, latent reservoir parameters were estimated from separate data sources. Table 1 contains all fixed values and initial guesses for estimated parameters. Importantly, this approach does not aim to define the absolute values of viral dynamic parameters, but rather to identify plausible sets of parameters that recapitulate the natural viral dynamics in humans so that these dynamical systems can be studied under the perturbative force of passive immunization with VRC01.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Sensitivity analysis variables to assess deviations from the initial model predictions.

https://doi.org/10.1371/journal.pcbi.1007626.t001

Estimated parameters included the death rate of susceptible cells δ_S, the probability of productively infectious provirus τ given infection, and several adaptive immune parameters: the killing rate of the adaptive immune compartment κ, the removal rate of the adaptive immune compartment δ_E (inclusive of death and exhaustion), the half-maximal saturation constant of the adaptive response E₅₀, and the recruitment rate of the adaptive response ω. Population distributions for each of the 6 estimated parameters are organized into density plots showing fits are achieved with relatively narrow parameter ranges across participants (Fig 2C).

We combined estimated parameters with fixed parameters to calculate an approximate basic reproductive number (see Methods). This value provides the average number of infected cells generated by a single infected cell at the start of infection [31]. If this value is below 1, viral production cannot overcome clearance, and infected cells eventually disappear. Our estimates (Fig 2D) are different from previous estimates [32] as a consequence of using a differently structured model.

Beyond the quantitative fitting procedure, the average parameter values generated several semi-quantitative predictions that were consistent with existing data (Fig 2E). Initial susceptible cell concentrations were approximately 1/4 of typical CD4+ T cell counts in humans, agreeing with observations that not all CD4+ T cells are permissive for HIV replication [33, 34]. These cells were depleted by roughly 200 cells/μL, roughly agreeing with losses in total CD4+ T cell counts [23]. Soon after detection, numbers of latently infected cells agreed with previous measurements [35, 36] both in magnitude (a total of approximately 100 infected cells per million CD4+ T cells) and in relative proportions of replication competent provirus (approximately 1 in 100 latently infected cells). The relative timing and magnitude of the adaptive immune response coarsely matched experimental observations: our modeled adaptive immune compartment rose several orders of magnitude within 2 weeks following first positive viral load; further, by viral set-point the adaptive immunity has saturated and the infected cell death due to adaptive immunity was comparable to the infected cell death due to cytopathic effects [37–39]—as calculated using the ratio of κE/(κE+ δ_I).

The initial model predicts excellent VRC01 protection against HIV acquisition throughout the AMP dosing interval

To predict VRC01 impact on viral dynamics in the dosing arms of the AMP studies, we incorporated a common PK/PD approach into our viral dynamics model.

As documented by two phase 1 trials studying the pharmacokinetic (PK) properties of VRC01 [24, 40], the clearance of VRC01 from serum is bi-phasic with an initial rapid distribution phase and a prolonged secondary elimination phase having a terminal half life of roughly 2 weeks. For simple integration with the viral dynamics model, we modeled the VRC01 PK as (1) where is the infused VRC01 concentration (μg/mL of VRC01), and k₁ and k₂ are the distribution and elimination rates (per day). Our simplified PK model provided excellent fit to the data and best-fit parameters from all individual fits were consistent (S1 Fig).

We also retrieved in vitro pharmacodynamic (PD) data from the 47 studies containing HIV/VRC01 neutralization data hosted in the LANL CATNAP database [25] (S2 Fig). A logistic, or ‘Hill’, function was used to model the fraction of virus neutralized (bound) over continuous VRC01 concentrations (see Methods for derivation). We define the neutralized (or bound) fraction at any time t after a VRC01 infusion as (2) so that when the VRC01 concentration is much above the 50% inhibitory concentration (IC50) all virus is neutralized , ν → 1, and when , ν → 0 and no virus is neutralized. Note, taking a logarithm of this value (logν(t)) admits the instantaneous inhibitory potential (IIP), previously used to calculate the efficacy of antiretroviral therapy [41]. Using Eq 2 and IC50 and IC80 data from S2A Fig we computed each Hill slope h. The values of h centered around h ∼ 1 (S2A Fig), in agreement with previous results [42, 43]).

We simulated PK curves and PD curves separately (Fig 3A & 3B). The lowest concentrations found in simulated AMP interval PK models were roughly 1 μg/mL (Fig 3A). In comparison, while there are some individual curves with higher IC50s, the median PD curve suggested most PD parameter sets would have an IC50 above 1 μg/mL (thick black line in Fig 3B). By combining these models (cartoon representation in Fig 3C) we simulated 1000 examples (Fig 3D) to demonstrate the range of expected neutralization given exposures to different CATNAP strains at different times following a VRC01 infusion (25 random examples illustrated to avoid over-plotting).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. The serum PK and in vitro PD provide an initial estimate for AMP protection.

A) Modeled PK kinetics using estimates from S1C Fig and Eq 1. B) Modeled PD dose-response curves using data from S2A Fig and Eq 2. C) Illustration of the modeling framework: a single PK concentration kinetic curve () is combined with a single static PD neutralization (IC50, h) to model individual neutralization at some time post infusion (Eq 3). D) 1000 simulations of neutralization during an AMP dosing interval were performed to calculate average quantities (plotted as dashed lines). 25 randomly drawn examples are plotted as lighter lines. E) Modeling the reproductive number modulated by VRC01 neutralization given an exposure at some time. Each PK/PD curve is rescaled by a single reproductive number (Fig 2, Eq 4). If exposure occurs when (threshold highlighted by thin black line), protection from infection is likely. Some parameter ranges allow for breakthrough infections, typically late in the dosing interval and in the lower dose arm (10 mg/kg, teal). F) The fraction of simulations having over time demonstrates that in vitro potency measurements and serum concentrations predict excellent prevention efficacy. In all simulations thick dashed lines indicate median values and colors indicate study arm.

https://doi.org/10.1371/journal.pcbi.1007626.g003

Next, we modified the viral dynamics to include a reduction in free virus (3) so that as neutralized fraction increases, virus decreases. The approximate analytical expression for the basic reproductive number (4) now can be used to coarsely model an exposure at some time t following an infusion with VRC01 (Fig 3E). By combining parameter values from the natural infection model fitting (Fig 2B & 2C) inserted into Eq 17 with the variety of PK/PD values inserted into Eqs 1 and 2, we calculated 300 possible combinations for . Many parameter sets achieved complete protection (), but certain parameter sets permitted infection within the first weeks after infusion (Fig 3E). If a functional exposure occurs at time t when , that individual is likely to be protected. We use the fraction of simulations where (thick dashed lines) to estimate protection probability (Fig 3F. These values remained above 75% in both arms even at the end of the dosing interval. We emphasize that these initial model results use serum PK and in vitro PD parameters, very likely overestimating protection in the real trials.

Simulated clinical trials using the initial model predict high prevention efficacy

To study further characteristics of protection including breakthrough viral loads, we adopted a stochastic formulation to study the initial model in a full AMP simulated trial. This model simulates a single infected cell introduced into a discrete number of susceptible cells. The probability of each new cellular infection, viral production, etc. is now considered a stochastic event that occurs with an average rate determined from the deterministic model (See Methods). Importantly, this allows for extinction (viral clearance) after exposure in the context of bNAbs. As part of model validation, we verified that the first positive after exposure in the placebo arm fell into empirical windows of 7-10 days [44, 45], suggesting the parameterization and stochastic simulator are reasonably representative of early infection. To model bNAb prevention, we required a precise definition for the initial conditions of the simulation. In non-human primate studies, systemic virus is found before clearance with bNAbs, suggesting that bNAbs do not completely prevent infection of all cells [13, 15]. Therefore, we defined ‘functional exposures’ as those for which virus passes the mucosal barrier and infects at least one susceptible CD4+ T cell. Defined this way, functional exposures have a >90% chance of causing infection in the natural stochastic model. We used the incidence rates and population sizes of the AMP cohorts to estimate the expected number of functional exposures in each study arm. For example, a 3% incidence rate per year (as estimated for for the Americas and Europe [17]) with a 500 individual trial for 80 weeks admits approximately 40 functional exposures (Fig 4A). We assumed incidence was identical in each dosing arm and simulated exposures as uniform over time, drawing the parameters for each exposure from a set of viral dynamic (Fig 2C), PK (S1C Fig), and PD (S2A Fig) parameters.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. A representative AMP simulation using the initial model (serum PK and in vitro PD).

For each dosing arm in the trial (placebo, 10 mg/kg, and 30 mg/kg, respectively), we simulated 40 functional exposures (based on a 3% incidence rate per year) and plotted viral load data: colored solid lines are simulated viral loads and open circles represent viral load measurements given the sampling frequency of the AMP trials. Dashed black line indicates a limit of detection at 30 copies/mL. Each exposure is initialized with a single infected cell carrying a randomly drawn viral strain (IC50) and a VRC01 concentration which follows randomly drawn participant PK kinetics. To normalize correctly, we account for the fact that in the placebo model (A), some functionally exposed participants naturally clear the infection due to variability in host immune parameters and stochastic simulation. Many functional exposures are blocked by VRC01. The trial admits excellent prevention efficacy: PE = 1 − 6/36, or 83%, in B) the low dose arm and PE = 1 − 3/36, or 92%, in C) the high dose arm.

https://doi.org/10.1371/journal.pcbi.1007626.g004

In a representative simulated trial example with the initial model (Fig 4), VRC01 had a strong effect on total prevention efficacy and a dose-dependent prevention efficacy was observed: PE = 1 − 6/36, or 83%, in the low (10 mg/kg) dose arm and PE = 1 − 3/36, or 92%, in the high (30 mg/kg) dose arm, respectively (see Eq 20). The approach of viral dynamics also highlights a possible challenge for AMP. In theory, some individuals will receive another dose of VRC01 after infection but before first positive, thereby suppressing viral load and presenting a complicated scenario for confirmatory measurements. In some simulations we observed undetectable viral loads and latently infected cells, representing the possibility for occult infection.

A framework to detect the mechanistic causes of deviations from the initial model predictions

There are several reasons why true prevention efficacy in the AMP studies might be lower than the estimated PE from initial model simulations. Therefore, we studied four mechanistic hypotheses for reduced prevention efficacy to determine whether deviations from the initial model could be inferred from characteristics of breakthrough infections (summary in Table 1).

First, concentration of VRC01 at sites of HIV exposure may be lower than concentrations measured in serum. Additionally, in vivo neutralization in mucosa might theoretically be reduced due to anti-idiotype antibodies or competitive exclusion [46]. Both mechanisms represent a global reduction in VRC01 neutralizing impact based on deviations from serum PK and/or in vitro IC50 [27]. To study these mechanisms, we introduce the ‘potency reduction factor’ ϕ. Mathematically, this parameter decreases the ratio of VRC01 concentration to IC50, multiplying the term . This simple scaling is plausible because a significant correlation has been identified between in vitro IC50 and prevention efficacy in non-human primate studies [27, 28] and two studies demonstrated that VRC01 in human serum retained its neutralization ability [40, 47]. Based on evidence from non-human primate studies, in which antibody concentrations in the vaginal fluid were 1-2 orders of magnitude lower than serum concentrations [26], we varied the potency reduction from 1 (no change from initial model) to 1000.

Second, circulating strains in the AMP trials may not be well-described by the laboratory isolates collected in the CATNAP database. In other words, VRC01 may have less breadth against isolates in the AMP study than those in the database. To study this impact, we parameterized a distribution containing sensitive and resistant strains. While recognizing resistance is a continuum, we chose sensitive and resistant strains as those with IC50 <1 and >50 μg/mL, respectively. This definition alludes to the definition of breadth as the fraction of viruses with an IC50 above the detection threshold [48]. We then chose a bi-modal distribution to parameterize exposing sequences. A fraction f of IC50s were drawn from a lognormal distribution centered on IC50 = 10³ μg/mL with variance 1 log and the remaining 1 − f were drawn from a lognormal distribution centered on IC50 = 10⁻² μg/mL with variance 1 log. We varied the ‘resistant fraction’ f from 0 (no resistant strains) to 1 (all resistant). The values for distribution modes allow for edge cases where at the low end f = 0, most sequences are neutralized easily, and at the high end most sequences can infect over the entire dosing window (given initial infusion concentrations usually less than 10³ μg/mL.

Third, viral inoculum size [49] and/or transmission-site anatomy [50] appears to affect the number of founder viruses. To address the possibility that functional exposures with more than one infected cell could overwhelm antibody protection and reduce stochastic burn-outs, we varied the number of active, productively infected ‘founder cells’ I_A,P(0), shorthanded as A₀ from 1-20.

Lastly, available in vitro data have a detection threshold and do not quantify IC50 values above 50 μg/mL. Thus, it is unclear exactly how much VRC01 would be needed to neutralize these strains. In addition, ‘incomplete neutralization’ has been documented [51, 52]. To model incomplete neutralization, we introduced x₅₀ to represent the maximum IC50, which we varied from the detection threshold of 50 μg/mL up to 10,000 μg/mL.

Global sensitivity analyses of simulated clinical trials were employed to test the impact of the sensitivity analysis variables on clinical outcomes. 100 parameter combinations were created using Latin hypercube sampling (LHS) from the ranges of the four sensitivity analysis parameters (Table 1). For each dosing arm, 100 functional exposures were simulated and five outcomes were recorded. Variable 1 is prevention efficacy (PE). Variables 2 and 3 are the mean and standard deviation of the first positive viral load in all detected breakthrough infections. The last variables (4 and 5) are the mean and standard deviation IC50s of breakthrough viruses as would be measured by an in vitro assay. These outcomes are referred to as measurable because they will be available in a post-hoc fashion from the AMP trials [17].

By correlating the values of the four sensitivity analysis variables against the five measurable trial outcomes Fig 5A & 5B, we observed several important relationships. The variable correlated most strongly to prevention efficacy was the resistant fraction f. The standard deviation of breakthrough IC50s was most correlated with the potency reduction factor ϕ. Both resistant fraction and potency reduction factor impacted the mean IC50 of breakthrough strains comparably. The maximum IC50 x₅₀ and the initial number of infected cells A₀ did not correlate strongly with any clinical outcomes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Global sensitivity analysis links unknown biological variables to measurable trial outcomes.

A) Results from 100 simulations of each dosing arm correlating sensitivity analysis variables (see Table 1 for definitions) against measurable trial outcomes. For each simulation, a value of each sensitivity analysis variable was chosen from a LHS sample. 100 functional exposures were simulated in each trial. Correlations agree in both dosing arms—teal 10 mg/kg, tan 30 mg/kg. B) Absolute values of Pearson correlation coefficients (marker shape indicates dosing arm) show that different variables correlate more strongly with different outcomes. Potency reduction factor and resistant fraction are the most strongly correlated overall, but each correlates more strongly with a different sensitivity analysis variable: potency reduction factor most strongly correlates with the standard deviation of breakthrough virus IC50, whereas resistant fraction correlates most strongly with prevention efficacy.

https://doi.org/10.1371/journal.pcbi.1007626.g005

To help interpret the impact of potency reduction factor and resistant fraction, we examined all breakthrough IC50s in a given trial (S3 Fig). Decreased potency allowed a wider distribution of strains (inclusive of lower sensitivity) to infect: thus the variability of breakthrough IC50 increased. Raising the resistant fraction increased the mean IC50 of breakthrough strains and decreased the variability.

A clinical-trial regression model (CTRM) accurately distinguishes partial efficacy due to reduced global potency vs. insufficient breadth against circulating strains

Based on the finding that outcomes were distinctly correlated with different biological variables (Fig 5), we developed a framework to assess the reverse correlation. That is, we sought to use the measurable trial outcomes to distinguish which unobservable biological variable was most-responsible for partial trial efficacy.

Using the global sensitivity analysis results, we performed linear multivariate multiple regression [53] using the measurable outcomes as the predictor variables and the sensitivity analysis variables as the dependent variables. Using only the five measurable outcomes of each trial, we assessed our prediction of potency reduction and resistant fraction using a cross-validation train/predict scheme. We randomly chose 40 simulations to train, and 10 to predict. This was repeated 5 times and errors were averaged. Training and testing was repeated for each dose arm. Estimates were reliably accurate (see true vs. predicted in Fig 6A). Absolute error (true−predicted, Fig 6B) was mostly unbiased, with median values for potency reduction slightly below 0. Relative error (|true−predicted|/true × 100) indicates estimation error is typically at or below 50% with some outliers (Fig 6C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. A clinical trial regression model (CTRM) uses trial outcomes to infer potency reduction and resistant fraction.

We simulated 100 trials in each dose arm, and performed 5-fold cross validation. We used sets of 40 randomly selected trials to train the CTRM, and then tested the CTRM to predict potency reduction and resistant fraction only using measurable trial outcomes from 10 different randomly selected trials. Errors were then averaged across 5 replicates. A) True vs predicted values for both variables. B) Total error (true−predicted) indicates reasonably unbiased estimates. C) Relative error (|true−predicted|/true × 100) indicates estimation error is typically at or below 50%. D) Heat maps indicate no single variable distinguishes reduced potency from increased resistance. E) Illustration of sequential inference using trial outcomes using an example simulation with a certain potency reduction and resistant fraction (ϕ = 100, f = 0.4, red dot). The inference proceeds by applying trial outcomes c_i sequentially to filter out unlikely simulated parameter combinations. Darker boxes indicate higher likelihood of a certain parameter combination. First positive viral dynamic outcomes (mean c₁ and std dev c₂) do not discriminate strongly. However, by filtering by prevention efficacy c₃ and IC50 outcomes (mean c₄ and std dev c₅), the true value is identified correctly.

https://doi.org/10.1371/journal.pcbi.1007626.g006

To demonstrate the approach intuitively, we plotted heat-maps relating the potency reduction factor and resistant fraction to each of the five measurable outcomes (Fig 6D). No single outcome is sufficient to identify the value of both sensitivity analysis variables. However, in practice, by using the value of each outcome found from an example trial, we can sequentially constrain the possibilities for the two values, until a reasonable estimate is developed. In Fig 6E, we perform this process by beginning with a simulation with a certain potency reduction and resistant fraction (ϕ = 100, f = 0.4, red dot). This value of potency reduction was chosen because it is plausible based on nonhuman primate bNAb concentration studies [26]. Then, using each heat map from Fig 6D, we calculated the percent error error where more precise estimates are represented by darker shades. Clinical outcomes c_i are applied sequentially. Viral dynamic outcomes (first positive mean c₁ and std dev c₂) rule out some parameter combinations, but are not very discriminatory. However, prevention efficacy c₃ and IC50 outcomes (mean c₄ and std dev c₅) strongly constrain the possible parameter combinations. Gray shades indicate the likelihood of parameter combinations Fig 6B & 6C. After all 5 outcomes have been used to filter out unlikely parameter combinations, the darkest square correctly identifies the true value.

Ethics statement

This paper uses de-identified data obtained previously and no new observations requiring patient consent or institutional review board approval have been performed.

Data and computational code availability

Code and data to replicate all analyses are hosted at https://github.com/dbrvs/AMP.

Mechanistic mathematical model for primary HIV infection in the context of bNAb intervention

We designed a mathematical model to describe primary HIV infection in the context of protection with broadly neutralizing antibodies (bNAbs). Our model (Fig 2A) has evolved from the canonical models of viral dynamics [66–68] and was developed previously to model SHIV kinetics [22].

The model contains susceptible cells (S; likely CCR5+/CD4+ T cells), which are produced constantly with rate α_S and are removed with rate δ_S. Free virus V_f infects these cells and generates infected cells I_s,p. These cells have a cell state and phenotype (subscript s and p). We use a vector notation to express the four possible cell states, that is I_s,p = [I_A,P, I_A,U, I_L,P, I_L,U], where active/latent is denoted by A, L and productive/unproductive is denoted by P, U. The infectivity and subsequent designation of cell state and phenotype is given by β_s,p = β × [τ(1 − λ), (1 − τ)(1 − λ), τλ, (1 − τ)λ]. Only productively infected cells produce virus, so π_s,p = [π, 0, 0, 0]. Both productive and unproductive infected cells die rapidly based on viral cytopathic effects and/or bystander killing [69, 70]. Latently infected cells persist with extremely long half-lives [71], but for this model we ignore reactivation from latency because it is rare and thus negligible during primary infection. Therefore, .

The adaptive immune system is modeled as an additional effector cell compartment E, which is inclusive of CD8+ T cells and other cells [72, 73]. Effector cells help control viral replication by killing all non-latent cells with rate κ_s,p = [κ, κ, 0, 0]. The recruitment rate of new adaptive immunity (ω_s,p = [ω, ω, 0, 0]) governs this recruitment and a saturation constant E₅₀ limits recruitment. Effector cells are naturally produced with the rate α_E to mimic the extremely low precursor frequency [74] and are also naturally removed (through death or exhaustion) with rate δ_E.

The model is expressed as ordinary differential equations, (5) where time derivatives are denoted by the overdot, β = ∑_s,p β_s,p, and a ⋅ b denotes the inner product. We incorporate the effect of VRC01 into the model above by adjusting the fraction of free virus V_f = V[1 − ν(t)] where ν(t) is the fraction neutralized VRC01 (described below).

Natural primary infection viral load data

Human primary infection viral load data were collected during the RV217 study; complete details are found in Ref. [23]. The study followed 2276 high-risk volunteers, and participants having detectable HIV-1 were followed such that plasma HIV-1 RNA could be quantitated twice weekly. For modeling purposes, we used data from 30 individuals from Thailand and Uganda whose viral load was observed at least 5 times, who had a single viral load measurement above 10⁵ copies/mL and had viral load measurements for viremia for at least 10 weeks following first positive viral load.

Estimation of model parameters for natural HIV primary infection

We fit the model to the RV217 viral load data using SciPy’s optimize package, which uses a least-squares Levenberg-Marquardt approach equivalent to maximum likelihood with a normally distributed variance for the data. 6 out of 14 parameters are fit, with the remaining values fixed (see Table 2 for fixed values and initial guess of fit parameters).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary of model parameters found in the model, Eq 5.

Parameters in () indicate initial values for estimation; others were held fixed.

https://doi.org/10.1371/journal.pcbi.1007626.t002

A more complete discussion of the identification of parameter values is provided in a previous work using a similar model [22]. Briefly, viral clearance rate γ was estimated from a human apheresis experiment. Infected cell death rate δ_I was estimated from human viral load decay data after initiation of ART. Viral production was estimated from a single cycle SIV experiment in macaques. Productive virus fraction was estimated in ex vivo cultures of human tonsil tissue. Susceptible cell rates (α_S and δ_S) were estimated from viral load data in a human treatment interruption trial. Importantly, that model has the same structure as our own. Adaptive immune rates are difficult to estimate because the model does not specify a specific cell phenotype. Human experiments have estimated anti-HIV T cell precursor frequencies, which we use to estimate α_E. in vivo microscopy was used to estimate CD8+ T cell killing rates, which we use to estimate κ. Maximal fractions of virus-specific CD8+ T cells in mouse LCMV experiments are 20–70%, which we combine with typical human CD8+ T cell concentrations to estimate E₅₀. The recruitment and death rates of the adaptive response are estimated from a T cell proliferation study in HIV-infected humans.

The rates governing susceptible and effector cells are the least well-understood, as well as the source (in our model) of the large variation in viral set-point equilibrium. Thus when fit, these values vary the most across participants. For instance, see the ranges of values of κ and δ_E shown in Fig 4C. There are some correlations among parameters and for the purposes of this work we do not claim to have identified absolute values. Rather, we have developed a variety of parameter sets for viral dynamics to be adjusted in the presence of VRC01.

For the deterministic model fitting, the initial viral load V(t = 0) is specified at a 30 copies/mL detection limit. Susceptible cells and adaptive immunity are initialized to equilibrium values, S(0) = α_S/δ_S, which assumes no real T cell depletion before first positive, and E(0) = α_E/δ_E, which assumes no HIV-specific immune response before first positive. The infected cells initial conditions are calculated using a quasi-static equilibrium so that the total number of infected cells is I(0) = V(0)γ/π. Then, we use the proportions to calculate I_s,p = I(0)β_s,p/β.

Modeling VRC01 pharmacokinetics (PK)

A detailed pharmacokinetic (PK) model for VRC01 was published by Huang et al. [24]. Here we adapt a simpler analytical formulation so that the concentration of VRC01 at any time t after a dose can be expressed as (6)

The terminal half-life of VRC01 can be calculated using the second phase decay constant t_1/2 = ln(2)/k₂.

Using SciPy’s optimize package, we estimated all four parameters based on data from HVTN 104 in which HIV-uninfected adults received multiple-dose intravenous VRC01 at 10 and 30 mg/kg every 4 or 8 weeks. Each individual (n = 12) received 3 infusions (S1A Fig). We used data for each infusion to estimate 32 and 31 parameter sets for each dose, respectively. S1B Fig shows all fits, and S1C Fig the estimated parameters. These parameter sets could be used to simulate continuous VRC01 trajectories (as in Fig 3A).

Curation of pharmacodynamic data

We use tabulated pharmacodynamic (PD) data including measurements of IC50 and IC80 from 782 viral isolates as tabulated in the LANL CATNAP database [25]. Some of these isolates are repeated, but measurements are not identical, so we included all data. Using IC50 and IC80 we calculated Hill coefficients for each isolate as (7)

In cases where both IC50 and IC80 are identical (h → ∞). In those cases, we set the Hill coefficient to a randomly selected hill coefficient taken from the distribution of all other isolates. Further, for these cases, the IC50 value is ‘saturated’ by experimental constraints. It may be that the virus is much more resistant, or that the real value of IC50 is much greater. Thus, in S2A Fig, we illustrate what these distributions might look like when saturated IC50s are drawn from a uniform distribution ranging from the saturated IC50 value up to 10³ μg/mL. Further, in Fig 3, we adjust the dose response curves to incorporate these theoretically possible resistant strains.

In S2B Fig we illustrate all PD data across virus clades.

Modeling VRC01 pharmacodynamics

To incorporate VRC01 neutralization into our model, we employ a stoichiometric neutralization model [30]. To derive this model, we assume VRC01 creates immune complexes C, proportional to the concentration of free virus V_f and the concentration of VRC01 , given the rate constant r as (8) where h represents a stoichiometric factor indicating how many VRC01 molecules must be present to neutralize a single virion. The total amount of virus V equates to the sum of the number of immune complexes and the number of free virions, that is V = C + V_f. Writing the fraction of bound virus as ν we have the relation (9) which can be solved for the bound fraction as (10)

We can then solve for a useful quantity, the concentration of VRC01 for which ν = 0.5, or the 50% inhibitory concentration IC50 = r^h. Thus, the neutralized fraction (ranging from 0 to 1) at some time t after infusion is (11)

Plots of this quantity after drawing random values for PK an PD are shown in Fig 3D. Then, during the stochastic viral dynamics simulation we interpret the reduction in free virus due to neutralization (12) to model the protective impact of VRC01.

Derivation of early infection basic reproductive number

The basic reproductive number defines the average number of infected cells generated by an infected cell in a completely susceptible population. We derive an approximate basic reproductive number in the context of bNAb therapy from Eq 5 by employing several approximations. First, we assume the viral dynamics are much faster than cellular dynamics so that at all times and thus (13)

Similarly, the binding of antibodies is assumed to be rapid compared to infection dynamics such that the free virus is (14)

We insert this expression into the differential equation for the active productively infected cells so that (15)

Next, we assume that early in infection the depletion of target cells is minimal, and the recruitment of adaptive cells has not occurred significantly. That is, we assume these cells remain at equilibrium S = S* = α_S/δ_S and E = E* = α_E/δ_E. Inserting these fractions, the concentration of infected cells can be factored out to leave (16)

Dividing through by δ_I + κα_E/δ_E allows us to identify the value (17)

As in the ‘survival function’ approach, we can think of this quantity as the product of the average lifespan of an active, productively infectious cell [δ_I+ κα_E/δ_E]⁻¹, the rate at which virus is produced from this cell π, the average lifespan of these virions 1/γ, d) the rate at which each virion infects susceptible cells to produce new infections βα_S/δ_S, and the fraction of these infections that are active and productive (1 − λ)τ [31].

We now define the time varying reproductive number that depends on VRC01 concentration and exposing strain as (18) such that the exponential solution to Eq 16 has the form (19)

Therefore, if an exposure occurs when , the number of active, productively infectious cells decreases exponentially. Likewise, if the number of infected cells grows exponentially. This important threshold can be used to estimate the circumstances when exposure will result in infection.

Modeling the infection probability during VRC01 infusion windows using the approximate early infection basic reproductive number

The basic reproductive number defines the stochastic extinction probability and therefore can be used to estimate the probability of infection, [58]. However, in our model grows monotonically as VRC01 declines. Therefore, while in practice at exposure typically results in extinction, extinction is not guaranteed.

Using the estimated parameters, we can calculate a range of possible values for (see Fig 2C & 2D). Putting in values for ν(t) based on PK/PD estimates, the value of is plotted in Fig 3E. Then, to calculate the probability of infection, we calculated the fraction of simulated trajectories at every time for which . This is plotted in Fig 3F.

Stochastic simulations of clinical trials

Using the ordinary differential equation system Eq 5, we developed a stochastic branching process simulation inspired by previous modeling of HIV [58, 59, 84]. Our implementation in Python, which employs the τ-leap approach [85], is publicly available. In each time interval Δt = 0.01 days, a Poisson number of each transition type occurs. For instance, the number of susceptible cells grows S = S + ΔS where ΔS is a Poisson random variable . Or, a certain number of active productively infected cells are created, , which meanwhile removes the same number of susceptible cells and a free virus. When an active productively infectious cell dies, a poisson number of virions with mean π are created. Such stochastic bursting has been shown to be similar to a model with continuous viral production [84].

The stochastic implementation allows us to simulate infection beginning with a single infected cell. Unlike the deterministic model, this approach permits burnout. Each of the initial state variable concentrations and rates involving concentration are converted to discrete values by multiplying by a volume. We choose this volume to be 10⁸ μL based on the observation that there is approximately 1-10 L of blood in an adult human and that there are approximately 10-100 times more T cells in lymph tissue than blood [86]. Our model for transmission is based on several elegant non-human primate studies [13, 15, 61, 87].

We assume that most sexual exposures to HIV are non-functional. In those exposures, sexual contact never results in infected cells in the recipient. This probabilistic process is not modeled, though it likely is responsible for the extremely rare per-coital transmission rate in humans [88]. Instead, we model functional exposures, defined as those in which a founder virus has crossed the mucosa and infected at least one local target cell in the recipient. These cells then are trafficked rapidly to lymph tissue where viral dynamics occurs.

To model the AMP trials, we assume repeated infusions of VRC01 every 8 weeks, and HIV viral load measurements (with detection threshold 30 copies/mL) every 4 weeks. Based on incidence data (and ignoring the influence of daily oral PrEP), we estimate approximately 40 functional exposures per trial arm [17]. This model does not distinguish between an individual who is doubly exposed and two individuals singly exposed. Then, for each ‘trial’ we simulated the three dosing arms, and computed prevention efficacy (PE) compared with placebo. We use a simplified calculation of PE using the ratio of infections that occur in a VRC01 treatment arm n^T to the number of infections in the placebo arm n^P. That is, (20)

A more sophisticated estimator using a Nelson-Aalen estimator is described in [17]. Another formulation like the ‘averted infections ratio’ could be used alternatively [89]. Our simple measure can be interpreted as the total fraction of participants at the end of the trial that would likely have been infected but were instead protected due to VRC01. For example, if 7 individuals were infected in the 10 mg/kg dosing arm and 39 were infected in the placebo arm the prevention efficacy would be 1 − 7/39, quoted as PE = 82%.

Global sensitivity analysis

Using the stochastic simulator of the clinical trials, we performed global sensitivity analyses on several key biological parameters (see Fig 5 and Table 1). For each simulated trial, we randomly chose values of all parameters simultaneously from ranges: potency reduction factor ϕ ∈ [1, 10³], resistant fraction f ∈ [0, 1], maximum IC50 x₅₀ ∈ [50, 10⁴], initial number of infected cells A₀ ∈ [0, 20]. Variables covering several orders of magnitude were drawn uniformly on a logarithmic scale. Clinical outcomes including prevention efficacy, mean and variance of first positive viral loads, and mean and variance of breakthrough IC50s were recorded.

The resistant fraction f parameterizes a simulated distribution of IC50s. The distribution is bi-modal with a fraction f of IC50s drawn from a lognormal distribution centered on IC50 = 10³ μg/mL with variance 1 log and the remaining 1 − f drawn from a lognormal distribution centered on IC50 = 10⁻² μg/mL with 1 log variance. We vary f from 0 to 1 because of we are uncertain how minor resistant variants may manifest in a realistic exposure.

Clinical trial regression model (CTRM)

Based on the correlation between biological variables and measurable trial outcomes (Fig 5B), we developed a multi-variate regression model that we term a clinical trial regression model (CTRM) to estimate the sensitivity analysis variables (v_j) from the five clinical outcomes (c_i). That is, for each biological variable, we identified the regression coefficients w_ij from the expression (21) where each z_j is a constant and is a normally distributed error term with a certain variance . The regression was accomplished using sklearn in Python.

The CTRM approach indicates that the full stochastic trial simulations can be reasonably well-approximated by a set of 5 regression coefficients for each sensitivity analysis variable.