Numerical approaches for the rapid analysis of prophylactic efficacy against HIV with arbitrary drug-dosing schemes

Pre-exposure prophylaxis (PrEP) is an important pillar to prevent HIV transmission. Because of experimental and clinical shortcomings, mathematical models that integrate pharmacological, viral- and host factors are frequently used to quantify clinical efficacy of PrEP. Stochastic simulations of these models provides sample statistics from which the clinical efficacy is approximated. However, many stochastic simulations are needed to reduce the associated sampling error. To remedy the shortcomings of stochastic simulation, we developed a numerical method that allows predicting the efficacy of arbitrary prophylactic regimen directly from a viral dynamics model, without sampling. We apply the method to various hypothetical dolutegravir (DTG) prophylaxis scenarios. The approach is verified against state-of-the-art stochastic simulation. While the method is more accurate than stochastic simulation, it is superior in terms of computational performance. For example, a continuous 6-month prophylactic profile is computed within a few seconds on a laptop computer. The method’s computational performance, therefore, substantially expands the horizon of feasible analysis in the context of PrEP, and possibly other applications.


Implementation and pseudo-code of Probability Generating System (PGS)
The PGS computes the extinction probabilities backwards in time, but uses an ODE solver to do so. Modern ODE-solvers automatically adapt the time step to achieve a default-, or user-defined accuracy. Hence, the time-step does not need to be defined a priori, but depends on the tolerance of the ODE-solver. The pseudo-code of the PGS is shown in Algorithm 1.

• Precomputation:
-Define an extra time τ that has to be added to the time interval of interest [T s , T e ]. In our predictions we use τ = 100hours ≈ 7 · t 1/2 , where t 1/2 ≈ 14.5hours denotes the half life of DTG. This choice of τ guarantees that the concentrations of the drug at time T e + τ are < 1% of the drug's trough levels.
-Pre-compute the pharmacokinetic profile D t for time interval [T s , T e + τ ] by solving eq. (11)-(13) (main article). The pharmacokinetic solution can be packed into a function so that the drug concentration at any time point within the interval [T s , T e + τ ] can be called.
-The reaction propensities for the single educt molecules are precomputed. Propensity a 5 should also be time-continuous like D t , and could be called "on the fly" using the function for D t .
-Pre-calculate the extinction probabilities P E (Y t =V , ∅), P E (Y t = T 1 , ∅) and P E (Y t =T 2 , ∅), for the three unit vectors in eq. (4) (main manuscript), in the absence of drugs, e.g. given in [1].
• Initialization: -The initial values of the ODE (eq. (23), main manuscript) are set to [P E (Y t =V , ∅), P E (Y t =T 1 , ∅), P E (Y t =T 2 , ∅)] and the time point to begin the integration to T e + τ .
The set of ODEs eq. (23) (main manuscript is solved backwards: in this work we used the ODE solver solve ivp in SciPy. The function to compute the pharmacokinetic profile D t , was called within the function that was passed to the ODE-solver, so that the value of a 5 can be updated during the integration process.
1 Input: start time point T s , end time point T e + τ 2 Result: extinction probabilities profile P E Y t =V , S , 3 Precomputation: 4 compute target-site pharmacokinetic profile D t for t ∈ (T s , T e + τ ); 5 compute values of a 5 based on D t : a 5 (t) for Y t =T 1 ; 6 compute the reaction probabilities a j , j ∈ {1, 2, 3, 4, 6} for single educt molecules respectively in the absence of drugs; 7 compute the extinction probabilities in the absence of drugs: at time point T e + τ 10 Solve the ODEs backwards: