article: To achieve the 90-90-90 goals set by UNAIDS, the number of new HIV infections needs to decrease to approximately 500,000 by 2020. One of the ‘five pillars’ to achieve this goal is pre-exposure prophylaxis (PrEP). Truvada (emtricitabine-tenofovir) is currently the only medication approved for PrEP. Despite its advantages, Truvada is costly and requires individuals to adhere to the once-daily regimen. To improve PrEP, many next-generation regimen, including long-acting formulations, are currently investigated. However, pre-clinical testing may not guide candidate selection, since it often fails to translate into clinical efficacy. On the other hand, quantifying prophylactic efficacy in the clinic is ethically problematic and requires to conduct long (years) and large (N>1000 individuals) trials, precluding systematic evaluation of candidates and deployment strategies. To prioritize- and help design PrEP regimen, tools are urgently needed that integrate pharmacological-, viral- and host factors determining prophylactic efficacy. Integrating the aforementioned factors, we developed an efficient and exact stochastic simulation approach to predict prophylactic efficacy, as an example for dolutegravir (DTG). Combining the population pharmacokinetics of DTG with the stochastic framework, we predicted that plasma concentrations of 145. 18 and 722. 23nM prevent 50- and 90% sexual transmissions respectively. We then predicted the reduction in HIV infection when DTG was used in PrEP, PrEP ‘on demand’ and post-exposure prophylaxis (PEP) before/after virus exposure. Once daily PrEP with 50mg oral DTG prevented 99–100% infections, and 85% of infections when 50% of dosing events were missed. PrEP ‘on demand’ prevented 79–84% infections and PEP >80% when initiated within 6 hours after virus exposure and continued for as long as possible. While the simulation framework can easily be adapted to other PrEP candidates, our simulations indicated that oral 50mg DTG is non-inferior to Truvada. Moreover, the predicted 90% preventive concentrations can guide release kinetics of currently developed DTG nano-formulations.
  HIV-1 continues to be one of the greatest public health challenges. While it is possible to suppress virus replication with antiretroviral combination treatment, the virus can persist in cellular and anatomical reservoirs for decades, precluding a cure [1–3]. Because of the inability to cure HIV, preventing its transmission is of utmost importance. The 90-90-90 target formulated by UNAIDS aim to end AIDS by 2030. An intermittent goal will be a drastic reduction of new HIV infections: While approximately 2. 1 million individuals became infected with HIV in 2015 [4], the intention is to reduce this number to 500,000 cases by 2020 and to fewer than 200,000 by 2030. Currently, pre-exposure prophylaxis (PrEP) for high-risk individuals is one of the five ‘pillars’ set by UNAIDS to achieve a drastic reduction in HIV infections. Of the available agents, tenofovir and emtricitabine (Truvada) have been extensively studied and were approved by the FDA and EMEA in 2012 and 2016 respectively. Most studies agree that Truvada can potently prevent HIV infection, if individuals adhere to the once-daily regimen [5,6]. However, major shortcomings of Truvada-based PrEP are its costs [7], the fact that it imperfectly protects from infection, and the necessity for daily drug intake, which often leads to inadequate adherence. These deficits of Truvada-based PrEP may be overcome by next-generation PrEP regimen, including more cost-efficient drugs, and drug formulations that require antiviral injections only every few month (see [8] for an overview of the PrEP pipeline). Quantifying prophylactic efficacy in the clinic is ethically problematic and extremely expensive, since it requires to conduct large (N > 1000 individuals) trials over very long time spans (years) to obtain statistically evaluable results. On the other hand, pre-clinical PrEP experiments only allow to study certain aspects in isolation. While PrEP efficacy is the result of a multivariate interplay of viral- and host factors there is a lack of translational tools to integrate available knowledge and to rationalize which agents and dosing regimen to test clinically. Having predictive models at hand that allow to rule out-, or prioritize certain regimen can avoid putting individuals at harm and greatly reduce clinical failure rates and associated costs. Our intention was to develop a method that integrates pharmacokinetic and pharmacodynamic (PK/PD), as well as viral characteristics to a priori assess the per contact prophylactic efficacy of arbitrary PrEP strategies against HIV. Recently, modelling approaches have been developed to predict the per-contact PrEP efficacy [9,10] by integrating various host and viral factors. Despite their advantages, these approaches conventionally neglect the pharmacokinetic–pharmacodynamic characteristics of HIV drugs and are therefore unable to simulate drug dosing, and dose frequency in order to ascertain how PrEP can most effectively be deployed. We recently developed a novel approach which fully integrates the pharmacology of nucleotide reverse transcriptase inhibitors (NRTIs) [11]. This approach, however, approximates virus extinction by its elimination probability during the first replication cycle. While this assumption is reasonable for moderately potent prophylactic compounds (like all investigated NRTIs), it underestimates the prophylactic efficacy of highly potent drugs and fails to predict efficacy in post-exposure prophylaxis (PEP). We overcome the aforementioned limitations by building on recent developments for the simulation of stochastic processes [12], implementing a numerically exact simulation approach to assess PrEP/PEP efficacy for time-varying drug concentrations (pharmacokinetics). The benefit of such integrative framework is to elucidate the relative importance of the distinct pharmacological and viral factors, such as the on- and offset of prophylactic protection, sensitivity to missed dosing events, -virus inoculum size and -timing of virus exposure. In the current work, we combine population pharmacokinetics with the novel stochastic simulation method to analyze different PrEP/PEP schedules with the second-generation integrase inhibitor dolutegravir.
  It has been shown that intrinsic stochastic fluctuations may be negligible in the so-called large copy number regimen (when Xt >> 1) and consequently, an ODE approach that models concentrations of molecules Xt/Ω, suffices [22]. In the current work, we re-formulate the above stated CME into a hybrid stochastic-deterministic (discrete-continuous) system Xt = (Yt; Zt), where Yt denotes the discrete-stochastic and Zt denotes the continuous-deterministic part. In our application, we view the antiviral pharmacokinetics (concentration time profile) as the external dynamical environment Zt, which we will model in terms of a low-dimensional set of ODEs (exemplified below). This is common practice in the pharmacometrics/systems pharmacology field and is supported by the fact that typically large quantities of drug molecules reach the target site. The target site concentrations Dt ⊆ Zt affect reaction propensities of our internal stochastic system, which models the stochastic events after exposure (viral replication and clearance reactions). The internal system Yt represents the state of the viral compartments represented by the number of free viruses V, early stage infected T cells T1 and late infected T cells T2, i. e Yt = [V, T1, T2, ]T, exemplified below. While this reduces the dimensionality, propensity functions of the stochastic sub-system are subsequently time-dependent, i. e. the continuous-deterministic drug concentration-time profile may affect reaction rates of the discrete-stochastic subsystem Yt on an infinitesimally small time scale. We will tackle this problem algorithmically, enabling the numerically exact estimation of prophylactic efficacy. Our goal is to estimate the prophylactic efficacy φ of particular medication regimen SD. The prophylactic efficacy denotes the reduction in infection probability per contact, φ (Y 0, S D) = 1 - 1 - P E (Y 0 | S D) 1 - P E (Y 0 | ⌀) (prophylactic efficacy), (1) where PE (Y0|SD) and PE (Y0|⌀) denote the virus extinction probabilities for a particular prophylactic regimen SD and in the absence of prophylactic drugs ⌀ respectively. The extinction probability is defined as P E (Y 0) ≔ P (Y t = [ 0 0 0 ] | Y 0 = [ V T 1 T 2 ]) (2) for t → ∞ and where Y0 denotes the initial state of the stochastic viral dynamics subsystem. In other words, the extinction probability is the probability that a stochastic trajectory eventually reaches the absorbing state [0,0, 0]T of the viral subsystem Y. Naturally, the infection probability is the complement of the extinction probability, PI (Y0) = 1 − PE (Y0). The terms PE (Y0|SD), PE (Y0|⌀) will be computed using a mathematical model of the viral dynamics that mechanistically considers the direct effects of antivirals on their respective target processes (outlined in a related article [23]), as well as individual drug pharmacokinetics following particular prophylaxis regimen. We adopted the viral dynamics model described in [24,25]. Long-lived and latently infected cells are only implicitly considered (outlined at the end of the section), motivated by the observation that transmitted viruses are not macrophage-tropic [26,27] and in line with related modelling approaches [9,10,28–30]. Although this model is a coarse representation of the molecular events happening during virus replication, it allows to accurately and mechanistically describe the effect of all existing antiretroviral drug classes on viral replication, as demonstrated in e. g. [31], and can be parameterized by in vitro and clinical data, Table 1. The modelled viral replication cycle consists of free infectious viruses V, uninfected T-cells, early infected T-cells (T1) and productively infected T-cells (T2). Early infected T-cells (T1) and productively infected T-cells (T2) denote T-cells prior- and after proviral integration respectively, where the latter produces virus progeny. During the onset of infection the number of viruses is relatively low and the number of uninfected T-cells Tu is fairly unaffected by virus dynamics [28,32,33]. We thus consider Tu = λT/δT to be constant over the course of simulations. The dynamics of the stochastic viral replication model after virus exposure are then defined by six reactions (the model is depicted in S1 Fig): a 1 = (CL + CL T · T u) · V t (clearance of free virus; V → *) (3) a 2 = (δ PIC + δ T 1) · T 1, t (clearance of early infected cell; T 1 → *) (4) a 3 = δ T 2 · T 2, t (clearance of late infected cell; T 2 → *) (5) a 4 = β · T u · V t (infection of a susceptible cell; V → T 1) (6) a 5 (D t) = (1 - η D (t) ) · k · T 1, t (proviral integration; T 1 → T 2) (7) a 6 = N T · T 2, t (production of virus; T 2 → V + T 2), (8) with C L T = (1 ρ r e v − 1) ⋅ β in eq (3), as outlined in [24] where ρrev = 0. 5 denotes the probability to successfully complete reverse transcription in the absence of inhibitors [34,35]. Free viruses are cleared by the immune system with a rate constant CL. Further, free viruses can be also cleared during unsuccessful T-cell infection CLT through the destruction of essential viral components of the reverse transcription-, or pre-integration complex [34,35]. The term β represents the lumped rate of infection of T-cells, including the processes of virus attachment to the cell, fusion and reverse transcription, leading to an early infected cell T1, before proviral integration. The term k denotes the rate by which early infected T1 cells are transformed into productively infected T2 cells, involving proviral integration and cellular reprogramming. The term NT denotes the rate of production of infectious virus progeny by productively infected T2 cells. The terms δ T 1 < δ T 2 denote the rates of clearance of T1 and T2 cells respectively and δPIC denotes the rate of intracellular destruction of the pre-integration complex. Parameters for the viral model are summarized in Table 1 and a mechanistic derivation of the dynamics from first principles is given in [24] (Supplementary Text therein). In this article, we study distinct prophylactic schemes with the second-generation integrase inhibitor dolutegravir (DTG). Integrase inhibitors act intracellularly by preventing proviral integration. In our virus dynamics model (eqs (3) – (8) ), this translates into a decrease in propensity function a5 by a factor (1 − ηD). Notably, DTG is active in its adminstered form (does not require biotransformation) and has physicochemical attributes that allow the unbound drug to rapidly cross the cellular membrane. We modelled the direct effect of dolutegravir using the Emax-equation [36] η D (t) = D t m IC 50 m + D t m, (9) where Dt is the target site concentration of the drug at time t and the term IC50 and m denote the drug concentration at which the targeted process is inhibited by 50% and a hill coefficient [37] respectively. Note that the equation above couples the stochastic viral dynamics subsystem Yt to the deterministic subsystem Zt, where the latter propagates the drug concentrations Dt ⊆ Zt. We used non-linear mixed effects modelling techniques [43] to derive a descriptive pharmacokinetic model that accurately captures the observed pharmacokinetic variability within- and across different patients. In this framework, both a minimal structural model f (θi, ⋅) is fitted to clinical data, alongside with statistical models describing the distribution of pharmacokinetic parameters θ within the population, as well as the measurement- or unexplained noise. Let D i, t be the measured plasma concentration of a drug in the ith individual at time point t. The likelihood of that measurement is defined through D i, t = f (θ i, t) · (1 + ϵ i, t) + ϵ ˜ i, t (10) where f denotes the solution of the structural model (a low dimensional set of ordinary differential equations; see below) that corresponds to the measurement. The vector θi contains the pharmacokinetic parameters for the ith individual. The variables ϵi, t and ϵ ˜ i, t denote proportional and additive error terms (measurement- or unexplained noise), which are typically assumed to be normal distributed, i. e. ϵ i, t ∼ N (0, σ 2) and ϵ ˜ i, t ∼ N (0, σ ˜ 2). The prior probability is typically assumed to be a multivariate log-normal distribution with log θ i = log θ + r i, (11) where θ denotes the vector of mean population parameters (fixed effects) and ri is normal distributed, i. e. r i ∼ N (0, Ψ). As mentioned before, direct computation of the extinction probabilities in eq (1) may not be possible as Yt still contains a prohibitively large number of states. In the following, we will utilize the results from a related article [23] in combination with the novel EXTRANDE algorithm [12] to compute the extinction probabilities for time-varying drug effects, i. e. taking drug pharmacokinetics into account. We consider a0 (Yt, Dt) = ∑k ak (Yt, Dt) to be the sum of the K reaction propensities changing the internal stochastic system. Note that in the exact SSA [47], no external input exists and therefore the propensities stay constant in between two reaction firings. In this case, the time to the next reaction event is exponentially distributed with parameter a0. In our case a0 (t) changes between two stochastic reaction firings because of pharmacokinetic inputs. Solving this problem requires to compute a0 (t) by numerical integration each time after a stochastic reaction has fired, which can be computationally expensive. Instead, in EXTRANDE, an upper bound B for a0 (t) is estimated and thinning techniques (rejection steps) are employed. For a look-ahead time horizon L, the upper bound Bt+L is chosen, such that B t + L ≥ a 0 (Y t + u, D t + u) (16) holds for all u ≤ L and assuming no stochastic reaction fires within the time interval t + L. When L is fixed, it is possible to solve for Dt+L, since it is assumed that stochastic reactions do not affect D. E. g. if the dynamics of D ⊆ Z are determined by a set of ordinary differential equations, numerical integration from t to t + L enables to predict Dt+L, which in turn allows to compute a0 (Yt+u, Dt+u). The internal stochastic system is then augmented with an extra reaction (a K+1th reaction), whose firing does not change the state of the stochastic subsystem Y. The probability of firing this extra reaction at time t + L is proportional to the ratio B t + L − a 0 (Y t + u, D t + u) B t + L. Obviously, it has to be guaranteed that eq (16) is true for the entire look-ahead time horizon L. On the other hand when Bt+L is chosen to be too large, many extra reactions will be fired (rejection/thinning step) and the algorithm becomes inefficient. A key to efficient simulation with EXTRANDE is therefore a good choice of Bt+L, which in turn depends on the look-ahead time horizon L. Codes were written in MatLab R2016b (MathWorks, Natick, MA; v. 9. 1, including optimization, parallel computing and statistics toolboxes). Individual pharmacokinetic model parameters for healthy individuals were drawn from the distributions defined by the parameter estimates from the final dolutegravir population pharmacokinetic model (Table 2; NONMEM $SIMULATION, n = 1000 individuals; eqs (12) – (14) ). We then simulated individual pharmacokinetic profiles for the prophylactic schedule SD under consideration using ode15s in MatLab. To simulate different adherence levels, a sequence of uniformly distributed random numbers with ri~U (0,1) was drawn and the ith dose was missed if ri > adherence level. The number of viruses to be inoculated were drawn from a previously parameterized distribution that accurately resembles the relation between transmitter virus loads and recipient infection probabilities [11]. In brief, we used a two-stage process: First, we sampled the viral load VL in a potential transmitter. Earlier work [11] (Supplementary Note 4 therein) showed that the virus load distribution in potential transmitter populations [51] follows a log-normal distribution, i. e. VL~logN (μ, σ) with μ = 4. 51 and σ = 0. 98. Secondly, we used the virus load in the transmitter to determine the number of viruses V0 entering a replication-competent compartment in the virus-exposed individual using a binomial model: V0~B (f (VL), rhomo) with rhomo = 3. 71 × 10−3 for homosexual exposure [11], f (VL) = ||VLc||, where ||⋅|| is the next integer function, and c = 0. 389 [52,53]. For PrEP simulations with different adherence levels, a time of virus exposure was randomly drawn within a 3 month interval. The corresponding drug concentrations at this time and the number of transmitted viruses reaching a target cell compartment were used as the initial states for EXTRANDE and simulated until stopping criteria were satisfied, illustrated in Fig 2 for virus infection (panels A & B) and -elimination (panels C & D). For ‘PrEP on demand’ simulations, the time of virus exposure was fixed as indicated in the corresponding graphics. In case of PEP, virus was inoculated as stated above and the stochastic viral dynamics were simulated in the absence of drugs until the time of PEP initiation to determine the initial condition of the system and henceforth simulated until a stopping criterium was reached. In total, for each prophylactic strategy, 5000 simulations were run.
  A total of 354 plasma concentration measurements from 56 individuals were used to build the population pharmacokinetic (PK) model for dolutegravir (DTG). Healthy volunteers (N = 17) contributed rich PK profiles with a total of 270 samples taken between 0 hours (pre-dose) and 216 hours after a final DTG dose. In addition, eighty-four measurements, randomly drawn between 1-25. 75 hours post-dose were available from HIV patients week 1,2, 3 and 4 weeks post-efavirenz switch. A two-compartment pharmacokinetic model best described the data (Fig 1A) and was fitted in a Bayesian context to fully capture inter-individual pharmacokinetic variability. Following multivariate analysis, allometric scaling (centered on 70kg) of weight was considered as a fixed effect in the model. Different values of apparent oral clearance (CL/Fbio) were estimated for DTG alone in healthy volunteers and in patients following 1,2, 3 and 4 weeks post-treatment switch. Residual variability was described by a combined proportional-additive model for healthy volunteers and a proportional error model for HIV-infected patients. All parameter estimates for healthy volunteers are displayed in Table 2. The model was used to generate PK parameters of virtual patients populations, whose PK-profiles are summarized in Fig 3B–3D alongside observed DTG concentrations, Fig 3C and 3D. As can be seen in Fig 3B, dolutegravir is rapidly absorbed after oral administration and maximal concentrations are achieved after tmax = 1. 58 hours (population 5–95% range 1. 53–1. 63). Pharmacokinetics reach a steady state after about 4 doses. During steady state, minimum- (pre-dose) and maximum concentrations were Cmin = 2918nM (1916–4336) and Cmax = 8471nM (6353–11331) for 50mg oral DTG and the half life of the drug was 14. 5h (5–95% range 13. 5–15. 9). Fig 4A shows the relation between the plasma concentration of DTG and its prophylactic efficacy after homosexual virus exposure. For these simulations, the number of viruses reaching a replication-competent compartment after homosexual contact were sampled from a previously parameterized virus exposure model [11] (Fig. 2 and Supplementary Note 4 therein). The estimated concentration ranges achieved at steady state for 2-, 10- and 50mg oral DTG once daily were 117–339,583–1694 and 2918–8471nM respectively as indicated on the y-axis of Fig 4A. Within this concentration range, the median prophylactic efficacy for once daily 2mg ranged from 43. 6 to 75. 7%. For 10mg, efficacies ranged from 87. 1 to 97. 5%, and for 50mg almost complete (99. 5 to 100%) protection was achieved. The estimated concentrations to prevent 50- and 90% infections, E C 50 (V ^) and E C 90 (V ^), were 145. 18 and 722. 23nM respectively. During pre-exposure prophylaxis, medication adherence may be incomplete. Fig 4B displays the prophylactic efficacy of once daily 50mg, 10mg and 2mg oral dolutegravir, considering varying levels of adherence (25-, 50-, 75-, 95- and 100% of doses taken). Viral challenges were simulated to randomly take place during a 3 month interval with inoculum sizes drawn from a previously parameterized distribution [11], see Fig 2 for two examples. The mean predicted prophylactic efficacies for 50mg with 25-, 50-, 75-, 95- and 100% adherence were 60% (95% confidence bounds: 55. 15–64. 84), 85. 54% (82. 58–88. 50), 96. 63% (95. 19–98. 07), 98. 88% (98. 04–99. 71) and 99. 36% (98. 73–99. 99), respectively. Notably, the prophylactic efficacy of 50mg oral DTG becomes saturated, and exceeds 95%, if at least 75% of the pills were taken. Conversely, 2mg and 10mg oral dolutegravir allow for considerable residual infection events and 2mg oral dolutegravir efficacy increases almost linear with increasing adherence levels. ‘PrEP on demand’ denotes a short-term pre-exposure prophylaxis, initiated only a few hours before a potential viral exposure. This strategy has recently been evaluated in the IPERGAY study using Truvada [54]. We simulated a dosing scheme similar to the IPERGAY protocol [54]: An individual at risk initiates PrEP only a few hours before a potential viral exposure and takes two consecutive doses 24 and 48hours after the first dose. Fig 4C depicts the predicted prophylactic efficacy of DTG when taken ‘on demand’. Population pharmacokinetic profiles for ‘PrEP on demand’ using 50mg are depicted in Fig 3A. The mean prophylactic efficacies for 50mg varied, depending on the timing of the first dose with respect to viral exposure, between 78. 63–83. 93% (95% confidence bounds: 75. 05–87. 05), for 10mg between 64. 49–73. 01% (59. 92–77. 02) and for 2mg it was 36. 86–46. 34% (30. 87–51. 90). The prophylactic efficacy decreased with decreasing dose and also with increasing time difference between the initiation of ‘PrEP on demand’ and viral exposure, clearly visible for 50 and 10mg. This trend is opposite to the trend for ‘PrEP on demand’ with Truvada [11]. A reason for this is the fast uptake of systemic DTG (compare Fig 4B), whereas the Truvada’s active components tenofovir diphosphate (TFV-DP) and emtricitabine triphosphate (FTC-TP) require intracellular phosphorylation after cellular uptake of the parent compound, which delays the time until maximal concentrations are achieved at the target site [31,55]. Since DTG does not need to undergo any chemical modification, and according to ‘free drug hypothesis’ the unbound intracellular concentration largely reflects the unbound plasma concentration, implying that the drug rapidly reaches the target site. Lastly, we wanted to assess the efficacy of 50mg oral DTG in preventing infection when taken as post-exposure prophylaxis (PEP). We assessed the prophylactic efficacy with regard to different durations of PEP and with regard to the timing of initiation after virus exposure in Fig 4D. Fig 4D indicates that 50mg oral DTG can effectively prevent infection (> 80%) when initiated shortly (within 6 hours) after exposure and when continued for as long as possible. The graphics indicate that the efficacy starts to drop when PEP is initiated later than 6 hours and when it is taken shorter than 7 days. Also, our simulations suggest that initiating the prophylaxis earlier has a more pronounced effect than prolonging PEP, arguing for the immediate start of PEP in case of known- or suspected HIV exposure. We previously estimated that once daily PrEP with Truvada provides ≈ 96% protection in fully adherent individuals [11]. While it is difficult to quantify PrEP adherence clinically [56], a surrogate measure is often calculated based on the percentage of individuals with detectable drug. Corresponding clinical efficacy estimates in apparently highly adherent individuals were 86-100% in the IPERGAY study [57], 58-96% in the PROUD study [6] and 96% in the Partners PrEP OLE study. In comparison, we predicted almost complete (99-100%) protection when 50mg DTG was taken once daily for prophylaxis. The VOICE [58] and FEM-PrEP [59] studies indicated that Truvada may not prevent infection in poorly adherent individuals, i. e. if ≈ 30% of individuals had detectable drug. In comparison, with 25% adherence to once daily 50mg DTG, we estimated about 60% protection and over 85% protection if at least half the drugs were taken. For ‘PrEP on demand’ with Truvada, we recently estimated [11] that about 74–92% infections can be averted, depending on the time of viral exposure relative to the initiation of Truvada dosing. The corresponding efficacy estimate in the IPERGAY trial was 86% [54], in line with our previous work. Herein, we predicted that ‘on demand’ PrEP with 50mg DTG is non-inferior to Truvada, providing 78. 63–83. 93% protection. Lastly, while PEP with Truvada is not recommended due to the slow intracellular accumulation of pharmacologically active NRTI-triphosphates, PEP with 50mg DTG can prevent about 80% infections when initiated no later than 6hours post exposure. In summary, our simulations indicate that prophylaxis with 50mg DTG is non-inferior to Truvada and that it may outperform Truvada when individuals’ medication adherence substantially deviates from a once daily PrEP protocol, as in the case of poor adherence and post-exposure prophylaxis.
  While pre-exposure prophylaxis with Truvada can prevent sexual HIV-1 transmission, it has severe limitations that are to be overcome by next-generation PrEP regimen [8]. However, drug candidate selection for next-generation PrEP and selection of administration schemes in clinical trials are prone to high failure rates. One reason is the poor translatability of animal- and ex vivo/in vitro experiments, whereas the statistical requirement of large sample sizes (N > 1000 individuals) and long trial durations to detect prophylactic efficacy clinically leads to exorbitant costs. Thus, there is an unmet need for tools to a priori assess the utility-, select and prioritise next generation PrEP candidates for human trials. Our intent was to develop a method that integrates pharmacokinetic and pharmacodynamic (PK/PD), as well as viral characteristics (inoculum size, timing of exposure) to a priori assess prophylactic strategies against HIV. Such integrative framework allows for the intelligent design of once-daily, episodic PrEP, as well as simulating the on- and offset of event-driven PrEP [54] or long-acting injectable PrEP formulations [60–62]. Moreover, it allows to assess how stable a prophylaxis is when individuals poorly adhere to the planned prophylaxis scheme, i. e. when individuals miss dosing events or start to take the drugs shortly before- or only after exposure to HIV. The use of population pharmacokinetic models allows for an accurate description of the observed pharmacokinetic variability within- and across patients, which can be used to assess the effectiveness of PrEP coverage in metabolically diverse population (Figs 3 and 4). While the attributes which make any compound favourable in the clinic extend beyond PK/PD, our approach is particularly useful to rule out- or prioritize PrEP candidates and/or strategies for further clinical investigation. To enable this ambitious goal we had to develop the theoretical basis that would enable us to accurately predict infection/extinction events for any arbitrary prophylaxis dosing regimen. The theoretical context is to compute solutions of the chemical master equation (CME), which is intractable due to the curse of dimension. Most naturally, instead of directly solving the CME, it is possible to sample trajectories using Monte Carlo techniques such as the stochastic simulation algorithm (SSA) [47] and to empirically reconstruct ℙ (X t = x i) from the trajectories. However, it is not clear when to classify a trajectory as an infection event. Moreover, stochastic simulations become inefficient when reaction propensities are high, which is typically encountered in a large copy number regimen, e. g. when Xt >> 1. To solve the latter issue (large copy number regimen), a number of hybrid methods have been proposed that partition a system X = (Y, Z) into species Y that are represented by a discrete-stochastic (CME) model and species Z that are approximated by their concentrations [63,64]. Since Z then evolves on an infinitesimally small time scale, a natural consequence of this partitioning is that stochastic propensity functions a (Yt, Zt) evolve between two stochastic reaction firings (t, t + τ). Thus, numerically exact computation requires to solve an integration problem to compute a (Yt, Zt+u). In sampling-based methods, i. e. the integral-based methods employed in [25,55], frequent initialization of numerical integrators can become a major computational burden (i. e. τ small). Direct hybrid methods [65–67] overcome this problem, however these methods still require to directly solve the CME part, which can be prohibitively large (involving thousands of states), limiting their applicability. In S2 Fig we show a comparison of simulations results obtained by integral-based (e. g. as used in [25,55]) sampling methods vs. the EXTRANDE method, as well as their respective run times. As expected, the results are identical, while the simulation time is lower for EXTRANDE. The method we are using is the recently developed EXTRANDE (extra reaction algorithm for networks in dynamic environment) method [12]. The EXTRANDE approach is based on a rejection, or thinning and allows for the numerically exact simulation of intrinsically stochastic kinetics embedded in a dynamical environment, i. e. the stochastic sub-system Y may constitute time-dependent propensities that are affected by the dynamic environment Z. We adapted EXTRANDE to allow the numerically exact prediction of infection/virus clearance in the context of time-varying drug inhibition. The key adaptation was to define proper stopping criteria that classify trajectories as infection events. To this end, utilizing results from a related article [23], we were able to compute and dynamically update an extinction simplex (eq (18) ), i. e. a part of the state space in which extinction can occur, see Fig 1. Whenever a trajectory leaves the extinction simplex, it is safe to classify it as an infection event and to stop the simulation. Classifying trajectories as infection events is not a straightforward choice: Using arbitrary thresholds Ytr = [Vt = c0, T1, t = c1, T2, t = c2]T can significantly distort predictions (e. g. ‘PrEP efficacy’) or increase computational run-time: While too small thresholds overestimate the number of infection events, large thresholds increase the computational run-time unnecessarily, since stochastic simulations become very inefficient when Y >> 1. Furthermore, there is no control- or knowledge of the numerical error made (the probability to falsely classify a trajectory as an infection event). Using the extinction simplex method proposed herein, it is guaranteed that the probability to falsely classify a trajectory stays below a user-defined threshold ε << 1. Moreover, the algorithmic run-time is optimal for providing this user-defined precision. Particularly during ‘PrEP on demand’ or ‘PEP’ simulations, dynamic adaptation of the extinction simplex can be algorithmically harnessed: Note that there is a positive relation between the size of the extinction simplex and the drug concentration D (compare Fig 1). After the last dosing events in ‘PrEP on demand’ or ‘PEP’ simulations the size of the extinction simplex shrinks as the drug concentration tapers, making it more likely that stopping criteria are met, which minimizes runtime. We used our framework to assess the utility of dolutegravir (DTG), which may be suitable for prophylaxis, since it has a good safety profile [68], a high resistance barrier [69] and a long half life in the blood plasma. Foremost, utilizing in vitro and in vivo parameters in our mathematical framework, we estimated that concentrations of E C 50 (V ^) = 145. 18 and E C 90 (V ^) = 722. 23 nM prevent 50- and 90% infections respectively. These concentrations can guide dosing and release kinetics of nanoformulated long-acting dolutegravir, which is currently in preclinical development [62]. Moreover, as soon as human pharmacokinetic data is available, our framework can easily be adapted to predict the PrEP utility of the long-acting formulation (by updating the pharmacokinetic model). As an example, we focussed on oral DTG herein, for which we had sufficient pharmacokinetic data to build a population PK model. Combining this model with the EXTRANDE framework, allowed to assess different prophylactic strategies: Overall, our simulations suggested that oral 50mg OD DTG may have a potential for PrEP with an estimated efficacy of 99 to 100% (perfect adherence). For comparison, we previously estimated that Truvada may prevent 96% of infections when taken once daily (perfect adherence) [11]. Our model suggests that DTG’s protective efficacy remains high (> 80%) even at adherence levels as low as 50%. This apparent forgiveness to poor adherence is due to DTG’s prophylactic potency rather than its halflife: I. e. if most DTG doses are taken, concentrations ranges are achieved where the concentration-prophylaxis profile is saturating (compare Fig 4A). Consequently, concentration changes do not proportionally translate into changes in prophylactic efficacy. In event-driven PrEP, we predicted that prophylactic efficacy reflects the drug profile in the blood plasma (compare Fig 3), which is characterized by a rapid absorption (tmax ≈ 1. 58 [h]) and slow elimination. This is in contrast to tenofovir-emtricitabine (Truvada) whose activity is not reflected by their plasma levels, since these drugs require conversion to intracellular diphosphates [31,55] to exert their antiviral activity, delaying the overall onset of activity. The main advantage of DTG over Truvada is in the context of post-exposure prophylaxis (PEP), where we predicted that it can potently prevent infection if initiated no later than 6hours post-exposure and taken for at least 5 days. We do not suggest to use single drug PEP with DTG; the sole purpose was to determine whether DTG would be effective, if individuals fail to take it before virus exposure. This kind of assessment allows to quantify the risks of prophylaxis in situations where the regimen is not taken as intended and moreover provides a scientific basis to include it in PEP multi-drug combinations. Our model has several limitations, but also a number of important advantages. Our simulations do not take into account drug concentrations at the site of mucosal exposure (e. g. cervix, rectum) [60,70]. These concentrations have, however, not been validated as targets for successful prevention or treatment, whereas data exist (albeit limited) for plasma drug concentrations. Instead, we modelled based on unbound concentrations, in line with the widely accepted ‘free drug hypothesis’, stating that unbound concentrations at the target site are responsible for pharmacological action. For drugs highly bound to plasma protein (> 90%), naturally since plasma protein concentrations are lower in tissues other than plasma, their total concentrations at sites other than the plasma are magnitudes lower [70]. Strikingly, however, the unbound plasma concentrations coincide with the unbound tissue concentrations [71], strongly arguing for the validity of the ‘free drug hypothesis’ [40,41]. Therefore, throughout the work, we assumed, according to the ‘free drug hypothesis’ [39] that the unbound concentrations in plasma and at the target site coincide. Note that dolutegravir is highly lipophilic (logP ≈ 2), enabling the unbound drug to rapidly cross cellular membranes, generating an equilibrium between the unbound drug on either side of the cellular membrane [42]. We estimated the probability of virus clearance (and the prophylactic efficacy φ) as a function of the number of viruses ultimately reaching a target cell environment after sexual exposure, and not as a function of mucosal exposure. The utilized virus exposure model [11] is calibrated to reflect the per-contact infection risks for typical transmitter virus loads and different routes of sexual exposure exposure. However, it should be noted that in an accompanying article [23], we also observed that increasing the inoculum size decreases the prophylactic efficacy, i. e. estimates of prophylactic efficacy depend on the route of transmission: For example, if exposure to HIV occurs via blood transfusion (large inoculum size), most prophylactic drugs may fail to offer protection. Our framework can be adapted or developed in a number of ways. The separate impact of treatment as prevention [72] (reduction of donor virus load to decrease contagiousness) versus prophylactic efficacy in the recipient individual can be simulated by calibrating the virus inoculum distribution [11]. The effect of PrEP on resistance transmission from a donor to a recipient can be incorporated in the framework by increasing IC50 in eq (9) (fold resistance) and possibly reducing certain reaction constants (fitness deficits). Likewise, the effects of PrEP on resistance emergence can be considered. However, during the early events after virus exposure (when infection can still be averted), the population size may be too small for resistance to appear de novo in the exposed individual. For example, a particular point mutation appears with probability 1 − (1 − μ) n, where μ ≈ 2. 2 ⋅ 10−5 [73] is the per base mutation rate of HIV per cell infection (= reverse transcription event) and n the number of cell infection events. According to these numbers, it requires ≈ 30000 cells to be infected for resistance to arise with 50% probability. Considering these numbers, a likely scenario for de novo resistance to appear is when PrEP had not been taken at the time of exposure, such that the infection expanded exponentially and a resistant mutant may have been generated at random. When PrEP is (re-) initiated at some later time it could provide the necessary pressure to select out the resistant type from the quasispecies population. Modelling these events is out of the scope of the recent work, as it requires distinct (and more coarse) simulation approaches, e. g. [25,74]. It is well known that the establishment of a latent reservoir is the major barrier to viral extinction during treatment [49] and this reservoir may be established as early as 3 days post infection [48,50]. We considered infection of long lived cells during our simulations, as outlined in the Methods section. I. e. , whenever long-lived cells became infected during simulations, viral extinction was considered infeasible. Notably, there are two sources: one with a half life of ≈ 14 days is responsible for the second phase of viral decay that can be observed in infected patients on combination treatment [75,76]. This reservoir has been attributed to infected macrophages [75,77]. Based on the mentioned half life of this reservoir, it would require T (x) = − ln (1 − x) ln (2) / t 1 / 2 days to eliminate a single cell with probability x under complete virus inhibition, e. g. ≈ 45 days to eliminate a single infected macrophage with 90% probability (x = 0. 9) and over 90 days to eliminate this reservoir with 99% probability. Latent infected T cells decline even slower, with a half life of ≈ 6-44 month [78–80]. This reservoir is partly responsible for the third phase of ‘decay’ and is assumed to prevent HIV cure during effective combination treatment. In summary, we have developed an innovative modelling approach to a priori assess prophylactic roll-out strategies by fully integrating individual PK/PD profiles and viral dynamics into a hybrid stochastic-deterministic framework. We used this framework to assess the prophylactic efficacy of the second-generation integrase inhibitor dolutegravir with respect to poor adherence, in event-driven prophylaxis and post-exposure prophylaxis. Overall our simulations showed that oral prophylaxis with 50mg DTG is non-inferior to Truvada and has profound advantages with respect to post-exposure prophylaxis. Moreover, we predicted that concentrations above EC90 = 722. 23nM can prevent > 90% infections after sexual exposure. These target concentrations can guide loading doses for novel long-acting nanoformulations of DTG [62]. By adapting the pharmacokinetics model, the framework can easily be used to predict the prophylactic utility of other candidate drugs currently under development, such as oral maraviroc (MVC) [81], and raltegravir (RAL) long-acting injectable rilpivirine [60] or cabotegravir [61], or it may be adapted to predict vaccine efficacy.
summary: In 2012, pre-exposure prophylaxis (PrEP) with Truvada was approved. It is considered one of the ‘five pillars’ by UNAIDS to drastically reduce HIV transmission. However, Truvada provides imperfect protection, is costly and individuals often fail to adhere to the once-daily regimen. Next-generation PrEP compounds, including long-acting formulations, are currently developed to improve PrEP. However, clinical trials using next-PrEP often fail. Since they involve many (>1000) individuals and long durations this incurs unacceptable costs, apart from individual tragedies. While animal- and ex vivo/in vitro experiments poorly translate into human efficacy, predictive tools are urgently needed that allow for PrEP candidate prioritisation. We developed an efficient simulation tool to predict the prophylactic utility of arbitrary dosing regimen. After developing population pharmacokinetic models for dolutegravir (DTG), we set out to predict its prophylactic utility in PrEP, PrEP ‘on demand’ and post-exposure prophylaxis. We found that 50mg DTG is non-inferior to Truvada in all aforementioned prophylaxis schemes. Moreover, we determined concentrations-prophylaxis profiles, which can guide release kinetics of currently developed DTG nano-formulations.
keywords: blood cells
  medicine and health sciences
  hiv prevention
  hiv prevention pills
  immune cells
  pathology and laboratory medicine
  pathogens
  post-exposure prophylaxis
  immunology
  microbiology
  retroviruses
  viruses
  preventive medicine
  pre-exposure prophylaxis
  immunodeficiency viruses
  rna viruses
  pharmacology
  public and occupational health
  white blood cells
  animal cells
  medical microbiology
  hiv
  t cells
  microbial pathogens
  viral replication
  pharmacokinetics
  prophylaxis
  cell biology
  virology
  viral pathogens
  biology and life sciences
  cellular types
  lentivirus
  organisms
title: Hybrid stochastic framework predicts efficacy of prophylaxis against HIV: An example with different dolutegravir prophylaxis schemes