Introduction

HIV-1 infection remains one of the major global health challenges with an estimated 33 million infected and a continuing spread [1]. Currently, an efficient vaccine remains to be developed, while at the same time the complete elimination of replication-competent virus within the host can not be achieved due to the persistence of the virus in inducible, latent cellular reservoirs [2, 3], as well as insufficient pharmacological suppression of actively replicating virus in some anatomical/cellular reservoirs [4, 5]. However, the current situation urges for methods that could bring the epidemic to a halt, or possibly end it. Currently, the most promising strategies are based on the use of antiviral drugs:

1. Pre-exposure prophylaxis (PrEP) [6–9] aims to protect uninfected individuals ‘at risk’ by decreasing the probability of infection upon virus exposure, e.g. [10]. PrEP may however be too costly to be broadly implemented in resource-poor countries [11].
2. Currently, the decision to initiate treatment against HIV is largely guided by CD4+ cell levels [12, 13]. However, the viral load, which is the primary determinant of infectiousness [14, 15], may be very high within the time-window between HIV infection and initiation of treatment. ‘Treatment for prevention’ [16] aims to put infected individuals on therapy as early as possible. This can reduce the infectiousness of a patient by decreasing within-host virus levels, which reduces the amount of transmitted viruses per contact and the probability of infection upon exposure. Analysis of the only completed clinical study to date, HPTN052 [16], estimated that ‘treatment for prevention’ may reduce the number of linked HIV-1 transmissions by 96% and the number of total HIV-1 transmissions by 89% relative to delayed treatment initiation and subsequently it was nominated as the “breakthrough of the year 2011” by the Science magazine [17].

In the aftermath of the HPTN052 trial, the cost-efficacy of ‘treatment for prevention’ was analyzed by many mathematical modeling approaches (reviewed in [18]). One problem is that most of these approaches focused solely on the epidemic level and did not model drug resistance development within the hosts, which indirectly assumes that the efficacy of ‘treatment for prevention’ is constant over time. However, because viral transmission is strongly correlated with viral levels in the transmitting individual [14, 15, 19–21], it is reasonable to assume that also the efficacy of ‘treatment for prevention’ is intimately connected with viral suppression. One major challenge during HIV treatment lies in the virus’ tendency to develop drug resistance [22], which in turn can lead to virus rebound and promote HIV transmission for the time it remains unrecognized. An earlier treatment initiation may thus demand an improved therapeutic strategy, that allows long-term control of virus replication (beyond the typical duration of a clinical trial). While sophisticated patient monitoring and timely treatment changes may allow to minimize windows of unrecognized viral breakthrough, they require significant monetary funds, good infrastructure, diagnostic facilities and the availability of alternative treatment options. Only few of these may be available in resource-constrained countries, where the requirement of resources may strongly dominate the possibility to implement a reasonable ‘treatment for prevention’ strategy. Obviously, scaling ‘treatment for prevention’ requires careful examination of various aspects and a policy maker should strike a proper balance between societal and individual perspectives [23].

This work addresses the scaling of ‘treatment for prevention’ by suggesting optimal treatment strategies for the long-term control of HIV within its host (as recommended by [24]). Optimality will be defined from a national economic perspective, taking into account that a diseased individual implies an economic loss. By considering the national economic perspective, we do not evaluate what should be done, but rather what is already worthwhile. However, we also evaluate the derived optimal strategies from an individual perspective and in terms of their utility in prevention, i.e. whether a strategy prolongs the life of an infected person and whether the risk of HIV onward transmission is reduced.

We hereby focus on two distinct approaches to handle treatment decisions: The first assumes that treatment decisions (i.e. when to change therapy) are made on an individual basis, guided by infrequent diagnostics (referred to as diagnostic-guided strategy). This represents a medical scenario in which a treating physician decides based on the diagnosed status of the patient that he encounters. The second approach suggests pro-active treatment decisions (referred to as pro-active strategy), i.e. does not require diagnostic ascertainment of the patients’ disease status. The two approaches are modeled and solved by two distinct mathematical frameworks. The former is addressed using the recently developed framework of ‘Markov Decision Processes with Rare State Observations’ [25]: For each disease state, it computes the optimal treatment and the next time of medical diagnostics, minimizing viral burden as well as treatment- and diagnostic costs. The latter approach (the pro-active strategy) is modeled as an open-loop switched system, where the decision to change the treatment depends on the initial disease state of the patient and the anticipated, (treatment-)induced stochastic dynamics up to some time t. The later strategy allows to switch treatment before drug resistance is detectable in the individual (pro-active) and may be easier to implement in resource-constrained settings, where poor infrastructure and the costs of diagnostics limit their applicability. By assessing these two distinct frameworks side-by-side, we can rigorously evaluate the different treatment paradigms in terms of their optimality. Algorithms to solve these problems were developed and are stated in the supplementary materials.

Several other groups have suggested optimal [26–28] or sub-optimal [29, 30] treatment strategies to mitigate drug resistance in HIV-1. All authors treated the underlying system deterministically, which fails to capture the intrinsic stochastic nature of HIV drug resistance development [31] and the time-scales on which drug resistance develops. None of the previous work focused on HIV prevention, and neither work questioned the analyzed treatment philosophy, either focusing on pro-active treatment switching strategies [26–28, 30], or diagnostic-driven strategies [29]. In contrast, we used a stochastic model of HIV long-term dynamics after drug application [25] to more realistically capture the underlying dynamics. Also, we evaluate different assumptions for the controllability of the disease dynamics, by evaluating the two different optimal control frameworks, which allows for an objective assessment of alternative treatment philosophies.

The manuscript will be organized as follows: We will extend- and parameterize the model introduced in [25] for our needs. After recapitulating essential theory for the diagnostic-guided strategy, we introduce the mathematical concepts behind the pro-active strategy, solve both optimal control problems and evaluate them with respect to monetary costs, patient survival and reduction of onward transmission. All algorithms that we developed to solve the optimal control problems will be provided in the S1 and S2 Text for the interested reader.

Materials and Methods

Within this work, we investigate optimal treatment strategies in silico by formulating- and solving two optimal control problems referred to as the optimal diagnostic-guided strategy and the optimal pro-active strategy. In general, an optimal control problem requires a mathematical model of the controlled process and a performance- or cost criterion. Likewise, our problem will be broken down into these ingredients.

Model of Controlled HIV Dynamics

The two addressed optimal control approaches share an identical model (Fig 1) that reflects the short-term dynamics of viral decay- and rebound (Fig 2), as well as the stochastic HIV long-term dynamics after drug application, see Fig 3. Within this work, we put a focus on viral kinetics and will only indirectly relate to the patient’s health. This is because we are interested in ‘treatment for prevention’ and particularly its efficacy in decreasing onward transmission, which is correlated with the viral load [19–21] and not necessarily with the immune status of the HIV infected patient.

Download:

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

Fig 1. Simplified HIV Model.

A: Transitions between copy number states n_C. B: Transitions in between viral strains M.

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

Download:

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

Fig 2. (Short-term) viral dynamics.

Left panels (A-C): Viral decay. Right panels (D-F): Viral growth. A: Data used for estimating viral decay parameters δ_h, δ_m. Blue circles indicate viral decay profiles from [41], green upward pointing triangles denote data from [42], black squares denote data from [43] and red downward pointing triangles denote data from [4]. Horizontal dashed lines and background shading indicates the assignment of the depicted data to the sets ℋ (> 4000 viral RNA/mL)), ℳ and ℒ (≤ 50 viral RNA/mL) of our model. B: When assuming 100% effective treatment (η = 1), the model shown in panel B is derived. This model is used to identify decay parameters δ_h and δ_m (circled parameters in panel B). C: Data-derived (error bars, dashed lines) and predicted (solid lines) probabilities of states ℋ and ℒ using the model in panel B with estimated parameters δ_h and δ_m. D: Data from treatment interruption trials [43–45] used for estimating viral growth parameters k_ℓ,∅, k_m,∅. Magenta diamonds indicate viral rebound profiles from [44], cyan pentagrams indicate data from [45] and black left-pointing triangles indicate data from [43]. Horizontal dashed lines and background shading indicates the assignment of the depicted viral growth data to the sets ℋ, ℳ and ℒ. E: We assumed the absence of treatment (η = 0), such that the model shown in panel E is sufficient to decribe the data and allows identifying growth parameters k_ℓ,∅ and k_m,∅ (circled parameters in panel E). F: Data-derived (error bars, dashed lines) and predicted (solid lines) probabilities of states ℋ and ℒ using the model in panel E with δ_h and δ_m and estimated parameters k_ℓ,∅ and k_m,∅. The parameter estimation procedure is exemplified in the Material and Methods section.

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

Download:

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

Fig 3. Long-term viral suppression.

Long-term data was used to estimate clinical drug efficacy η(a₁, {WT, R2}), η(a₂, {WT, R1}) and rates of drug resistance emergence μ_R1,∅, μ_R2,∅, using the model depicted in Fig 1, after parameters for viral growth and decay were estimated from data in Fig 2. A: Predicted (solid black line) and clinically observed probability of viral suppression (states ℒ; ≤ 50 viral RNA/mL) after treatment with efavirenz (EFV) based HAART (first line therapy). Clinical data was derived from [46] (red dots), [47] (orange squares), [48] (green diamonds), [49] (magenta upward pointing triangles) and [50] (blue downward pointing triangles). In all studies, the NRTI backbone consisted of 3TC + AZT. B: Goodness-of-fit plot for first line therapy. C: Predicted (solid black line) and clinically observed probability of viral suppression (states ℒ; ≤ 50 viral RNA/mL) after treatment with ritonavir boosted lopinavir (LPV/r) based HAART (second line therapy). Clinical data was derived from [51] (brown right-pointing triangles), [52] (cyan pentagrams) and [52] (grey left-pointing triangles). In all studies, the NRTI backbone consisted of a deoxycytidine analog + abacavir or tenofovir or stavudine, reflecting clinical practice (the exact choice of the backbone may depend on prior exposure [13]). D: Goodness-of-fit plot for second line therapy.

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

State space.

HIV can be successfully suppressed if drug resistance does not develop. Thus, any model that aims to represent the long-term HIV dynamics upon treatment should include drug resistance development. The process of drug resistance development denotes an intrinsically stochastic process, which is determined by random mutation events (point mutations, recombinations). Long term HIV-dynamics in the context of drug treatment may therefore be dominated by these intrinsically stochastic events [31], necessitating stochastic modeling approaches [32–34]. The fundamental evolution equation for intrinsically stochastic kinetics is the chemical master equation (CME). Each state described by the CME comprises a combination of discrete numbers of individuals of the respective species (e.g. viral strains), resulting in state space dimensions ℕ₀×ℕ₀×…×ℕ₀, i.e. [35, 36]. A major mathematical drawback is the fact that the CME cannot be solved directly due to this complexity. Therefore, a modeler can either approximate the solution of the CME by Monte-Carlo schemes [37], aim at hybrid approaches [38–40], which can yield particular characteristics of the CME, or perform a state space reduction (lumping). In this manuscript, we adapt a model [25] that relates to the latter approach. For this model we can solve the coarse-grained CME directly when computing optimal control strategies.

In brief, the HIV model contains four lumped viral copy number states for each of the four virus strains. The set of states 𝒮 thus has dimension 4⁴ = 256 states + 1 [patient death]: If the respective virus type is absent, we denote the respective state by 0, if it is present in low copy numbers, i.e., for < 50 virus copies/mL blood (detection limit of assays used in the clinic), the respective state is denoted by ℓ, for medium copy numbers between 50 and 4000 virus copies/mL blood we denote the lumped states by m and for high copy numbers with more than 4000 virus copies/mL blood, it is h. This coarse graining is in line with the levels of virus produced in the distinct cellular reservoirs of HIV, see e.g. [34]. The following four viral strains M are considered: a strain WT (wild type) that is susceptible to all treatment lines, a strain R1 which is susceptible to a second treatment line, but unaffected by (resistant to) the first treatment line, a strain R2 that is susceptible to the first treatment line, but unaffected by the second, and a highly resistant strain HR, which is resistant to all treatments. In order to describe a virologic state x we choose a compact vector notation of the form where n_C ∈ {0, ℓ ,m, h} denotes the viral copy number of each viral strain WT, R1, R2 or HR. For example, the state $x=[\ell ,m,0,\ell ]$ describes the situation of a ℓow number of wild type strains, a medium number of R1-mutants, the absence of R2-mutants, and a ℓow number of highly resistant viruses. Mutations from one strain to another can give rise to novel viral populations, as shown in Fig 1.

Control actions.

The actions describe ‘what the controller can do to influence the system’. In terms of HIV therapy, a physician can e.g. choose what treatment(-line) to apply and when to change it. In resource-constrained settings, only few treatment lines are available. In the case of South Africa these may include a first- and a second-line therapy [13]. Taking these considerations in account, we consider two distinct treatment lines (actions) a₁,a₂ ∈ 𝒜. Each action a ∈ 𝒜 induces unique disease dynamics, related to a unique Markov Jump Process that is entirely determined by its infinitesimal generator L_a. The entry L_a[x,y] ≥ 0 represents the rate of transition from state y ∈ 𝒮 to state x ∈ 𝒮, y ≠ x, given an action a and it holds that L_a[y,y] = −∑_{x ≠ y} L_a[x,y]. We define a probability space Ω and let p ∈ Ω denote a probability distribution vector on the state space 𝒮 with the entry p[x](t) referring to the probability of being in the state x ∈ 𝒮 at time t, i.e. (1) where ℙ is the probability measure. Obviously, the number of components of a probability vector p is equal to |𝒮|. For a given action a ∈ 𝒜, the dynamics of the probability vector are given by (2)

The above equation is known as the Master Equation. We introduce the transpose of the transition matrix on 𝒮 for some time lag τ and action a (3) where e denotes the matrix exponential. The component T_a,τ[x,y] refers to the transition probability from state y to state x for a time lapse τ under the application of action a and will be used later in the cost functionals of the closed-loop optimal control problem (diagnostic-guided strategy) and the open-loop optimal control problem (pro-active strategy).

Generator entries.

The distinct treatments a ∈ 𝒜 are related to distinct generators L_a of our HIV-model. The basic transitions between copy number states for each viral strain M, n_C(M), are shown in Fig 1 and exemplified for the highly resistant strain HR below. (4) (5) where * indicates an arbitrary number of the respective virus strain (WT, R1 and R2 in the example above). The parameters k_ℓ,a and k_m,a denote the reaction propensities of going from copy number ℓ to copy number m and from copy number m to copy number h respectively (viral growth), which are decreased depending on the treatment a ∈ {a₁, a₂} because treatment essentially suppresses viral growth. The parameters δ_m and δ_h denote the reaction propensities for going from copy number m to copy number ℓ and from copy number h to copy number m respectively (virus elimination). The parameters d_h > d_m > d_ℓ denote the propensity for the death of the patient. We assume that high viral burden (states h and m respectively) increases the risk of death, whereas d_ℓ equals the propensity for “natural death”. The propensity for death was computed according to d = 1/(residual life expectancy), and is exemplified in [25].

The considered transitions between viral strains M are depicted in Fig 1. Specifically, transitions between viral strains generate a ℓow number of viral particles from either a medium or high number of viruses belonging to a distinct strain. Note, that transitions between viral strains may involve several distinct point mutations (indicated by blue and red bars in Fig 1B). Exemplified for the wild type strain WT those are: (6) (7) (8) (9) where the first two lines indicate drug resistance arising from the wild type strain and the remaining two lines indicate transitions from resistant strains yielding the wild type strain. The parameters μ_R1,a and μ_R2,a denote the propensity for the emergence- and disappearance of drug resistance to treatment 1 or 2 (a₁,a₂), respectively, emanating from copy number state h or m. Note, that we consider only the following transitions: WT ↔ R1, WT ↔ R2, R1 ↔ HR and R2 ↔ HR, which is motivated by the fact that a direct transition from WT ↔ HR is very unlikely, because the genetic distance between the two viral strains is too large to be overcome at once.

The effect of treatments a₁ and a₂ on the viral growth & transition rates is considered in the following way: (10) (11) (12) where M ∈ {WT, R1, R2, HR} denotes the strain of the reactant virus. $\tilde{M}\in \{\mathrm{\text{WT}},\mathrm{\text{R1}},\mathrm{\text{R2}},\mathrm{\text{HR}}\}$ denotes the event related to a particular drug resistance emergence/disappearance, see Fig 1B. The parameter η(a,M) denotes the efficacy of treatment a on the reactant viral strain M; i.e. if strain M is susceptible to treatment a ∈ {a₁, a₂}, then 0 < η(a,M) ≤ 1, and if the viral strain M is insusceptible to treatment a ∈ {a₁, a₂} then η(a,M) = 0. In the absence of medical intervention a = a_∅, η(a,M) = 0. Therefore, the parameters k_ℓ,∅, k_m,∅ and ${\mu }_{\tilde{M},\varnothing }$ denote the growth rates and respective transition rates in copy number states m and h in the absence of intervention, as shown in Table 1.

Download:

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

Table 1. Parameters of the HIV-model.

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

Parameter estimation.

In order to estimate model parameters, we proceeded in a step-wise approach: We first estimated parameters related to viral decay (δ_h, δ_m) and then used these estimates in order to estimate parameters related to viral growth in the absence of treatment (k_ℓ,∅, k_m,∅), using data from [4, 41–45]. Finally, we used the estimated decay- and growth parameters along with data on the long-term (> 2 years) suppression of HIV-1 in order to estimate parameters related to the drug efficacy (η(a₁, {WT, R2}), η(a₂, {WT, R1})) and to drug resistance development (μ_R1,∅, μ_R2,∅) [46–52].

Parameters were estimated in MATLAB using lsqcurvefit by minimizing the following weighted least squares criteria, with θ denoting the set of estimable parameters. (13) where π[x](t_i) denotes the data-derived probability distribution on the model-defined state-space (computed using the ecdf function in MATLAB), p[x](t_i,θ) defines the solution of Eq (2) for time t_i with parameter set θ and ω_i denotes the weight parameter. Parameter estimation was performed 50 times respectively with random start parameters to verify the convergence to globally optimal parameter estimates θ^*.

Viral decay

A total of 311 data points from 31 patients and 4 independent clinical studies were available from [4, 41–43], which accurately assess the dynamics of viral decay after initiation of treatment (see Fig 2A). For the data analyzed, we assumed 100% effective treatment (η = 1), as proposed by others who estimated viral decay parameters [41, 53]. The lumped viral model (see Fig 1) then further reduces to the model shown in Fig 2B, which allows to identify decay parameters δ_h and δ_m. The data-derived probabilities π[x](t_i) were computed as 1− the cumulative probability to leave set ℋ (> 4000 viral RNA/mL) and the cumulative probability to enter set 𝓛 (≤ 50 viral RNA/mL). Error bars were computed using Green’s formula. In line with the data, we assumed that the initial HIV virologic status is represented by high copy numbers of susceptible virus.

Viral growth

A total of 89 data points from 17 patients and 3 treatment interruption trials [43–45], was used to estimate viral growth parameters k_ℓ,∅ and k_m,∅. In line with the data, we assumed the absence of treatment (η = 0), such that the model shown in Fig 2E is sufficient to describe the data. Data-derived probabilities were computed as 1− the cumulative probability to leave set 𝓛 and the cumulative probability to enter set ℋ, respectively, and error bars were computed using Green’s formula.

Drug efficacy and -resistance

Using the full model (Fig 1), we estimated parameters relating to the clinical drug efficacy of both treatment lines η(a₁, {WT, R2}) & η(a₂, {WT, R1}) and rates of drug resistance emergence μ_R1,∅ and μ_R2,∅.

In analogy with the South African treatment guidelines, we assumed that the first-line therapy consists of efavirenz (EFV) + zidovudine (AZT) + lamivudine (3TC). Long-term studies usually evaluate the probability of viral suppression, which is defined in terms of undetectable virus loads (≤ 50 viral RNA/mL). Translated to our model, this refers to the condition in which all viral mutants are in state ℓ or absent; i.e. $[\le \ell ,\le \ell ,\le \ell ,\le \ell ]$, which we denote by the set of states by 𝓛. Probabilities of viral suppression from 5 clinical studies [46–50] were used for parameter estimation. As a second-line treatment we assumed a ritonavir-boosted lopinavir (LPV/r) based HAART, see [13]. Since the exact choice of the NRTI backbone may depend on the prior exposure of the individual patient, we used data evaluating the long-term efficacy of LPV/r + an NRTI backbone consisting of a deoxycytidine analog + stavudine [51] or abacavir [52] or tenofovir [52].

All model parameters are shown in the Table 1. The original data and model predicted dynamics of viral decay and -rebound are shown in Fig 2 (A: raw viral decay data; B: model to evaluate viral decay; C: model-predicted vs. clinical decay profiles; D: raw viral growth data; E: model to evaluate viral growth; F: model-predicted vs. clinical growth profiles). Data for the long-term control of HIV-1, predicted dynamics and goodness-of-fit are shown in Fig 3A–3D for the two treatment lines (a₁ and a₂). As can be seen in Figs 2 and 3, the model appropriately captures both the short-term viral dynamics, as well as long-term dynamics of viral suppression.

Cost assignment.

Public health initiatives are often constrained by available funds. The countries with the highest HIV burdened are also among the poorest and financial commitments from donors have stagnated or decreased [54] in recent years. Thus, the requirement of resources may strongly dominate the policy making process in a resource-constrained context. Because of these conditions, we designed the performance criterion from a national economic perspective.

The performance criterion valuates the induced system dynamics and controls, i.e. the viral status of the patient and the costs of treatment. We will consider both the direct costs due to the applied treatments c_𝒜 and indirect costs due to the virologic/health status of a patient c_𝒮. Our analysis will be conducted from a country’s public health-care/monetary perspective. Therefore, the costs related to the different states c_𝒮 will be computed based on the average productivity loss pL(n_C) times the average daily monetary contribution of one individual (assessed in terms of the daily per capita GDP), i.e. c_𝒮(x) = pL(x)·GDP, with $\text{pL}(x)=\underset{n_C}{\max }\text{pL}(n_C)$, which implies that the total virus load reflects the cost of the individual infection status at any point in time. Death is interpreted in terms of a complete loss in productivity. Furthermore, we take diagnostic costs into account, which applies only in the diagnostic-guided strategy, the standard of care and the HPTN052 protocol (the latter two are modeled for comparison). The cost of diagnostics will be set to a fixed value and closely reflect the cost of a drug resistance test for the diagnostic-guided strategy and the cost of a virus load determination in the case of the standard of care and the HPTN052 protocol.

The integration of momentary/running costs yields the objective function (performance criterion) for the optimal control problem. While performance criteria generally depend on the particular application at hand, we decided to consider expected discounted costs on an infinite time horizon. We chose an infinite time horizon, because HIV treatment does not have a previously known endpoint (i.e. time of death). At the same time, a differentiated weighting of immediate and later costs is reasonable due to an upper limitation of life expectancy and aspects of inflation. Costs arising at time t > 0 are thus weighted by a discount factor 0 < e^−λt < 1. In this regard, the concrete choice of a discount factor λ will depend on the presumed annual inflation in the considered setting. For all calculations, we consider the inflation rate in South Africa as a representative of a resource-constrained country with a large HIV burden, see Table 2. The discount factor also guarantees convergence of the cost functional and therefore allows the numerical solution of the optimal control problem.

Download:

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

Table 2. Cost parameters for South Africa.

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

The costs per unit time comprise both the direct costs due to the applied treatments and indirect costs due to the virologic/health status of a patient. Thus, we can write (14) where c_𝒜:𝒜 ↦ [0, ∞) is the direct cost of action per unit time and c_𝒮:𝒮 ↦ [0, ∞) is the indirect cost produced by the state per unit time with parameters given in Table 2.

We define a cost function (15) which denotes the expected discounted costs for the time interval (0,τ] when starting in state x and choosing an action a ∈ 𝒜 for propagation of the stochastic process for the entire interval τ. Further, we define the cost vector 𝒦_a ∈ ℝ^|𝒮|, where its x^th component denotes the direct and indirect cost per unit time for the state x ∈ 𝒮 as shown below (16) so that it holds that (17) where the vector φ_x denotes a point-distribution, i.e. a single realization X_t of the Markov Jump Process. If the initial state is described by an arbitrary distribution p on the state space 𝒮, we get (18) where p[x] denotes the probability of the x^th state.

Results

Optimal Treatment Strategy

The optimal diagnostic-guided strategy is given in S1 Table. In brief, for the considered parameters (Tables 1 and 2), it is suggested to use the first-line treatment a₁ in all states, except those where the virus is resistant against treatment a₁, but susceptible to a₂. In the later case treatment line a₂ is suggested. In line with this treatment recommendation, patient monitoring is only suggested as long as the patient is infected with drug-susceptible (“wild-type”) virus. If the patient has a high or medium virus load, the next diagnostic test should be within 25 days, if the patient has a ℓow/non-detectable virus load, after 152 days.

These results may indicate that the cost for diagnostics is too high in relation to the economic benefit resulting from more close monitoring and informed treatment adaptation (this will be discussed later in the Discussion). An exemplary trajectory that highlights the treatment strategy is shown in Fig 4A. The blue line indicates a patient-specific trajectory. The filled black marks indicate the times when a diagnostic test is performed and the background shading indicates the applied treatment (white: a₁, gray shading: a₂). In the example, the patient initially has a high copy number h of wild type (WT) virus, while none of the drug resistant viruses are present. This state is represented by the vector notation $X_{t_0}=[h,0,0,0]$. For this state, the optimal treatment policy (see S1 Table) suggests to use treatment a₁ and to perform the next diagnostic test in 25 days (the second black marking in panel Fig 4A). At the next diagnostic test, the patient is in state $[m,0,0,0]$ for which continuation of treatment a₁ is recommended and the next diagnostic test is scheduled after 25 days (the 3^rd–9^th black marking in panel Fig 4A). In the following, the virus remains suppressed, with a small detected ‘blip’ after about 500 days. After about 600 days of treatment, during the time lapse between diagnostic tests, the a₁ resistant strain R1 emerges. Notice transitions from the state $[m,0,0,0]\to [m,\ell ,0,0]\to [m,m,0,0]$, then $[\ell ,m,0,0]$ and finally $[\ell ,h,0,0]$ in the Fig 4A, where the copy number of a₁ resistant strain R1 increases from a ℓow copy number to a high copy number (virus rebound after resistance development). At the time point of the next diagnostic (at around 700 days), the emergence of resistance is identified $[\ell ,h,0,0]$ and a switch to treatment a₂ is suggested (marked by gray region in Fig 4A). After the initiation of treatment a₂, a transition to state $[\ell ,\ell ,0,0]$ can be observed in the trajectory, which implies a decrease in the a₁ resistant strain (viral suppression).

Download:

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

Fig 4. Disease progression for the diagnostic-guided strategy (individual trajectory, panel A) and pro-active strategy (probabilistic measure, panel B).

The white region denotes application of treatment a₁ and the gray region denotes the application of treatment a₂. We assumed that the initial HIV virologic status is represented by a treatment-naive patient with high copy number of wild type virus [n_C(WT) = h,n_C(R1) = 0,n_C(R2) = 0,n_C(HR) = 0]. In panel A, the blue line represents a stochastic realization of HIV dynamics in a single individual treated with the diagnostic-guided strategy and black dots indicate diagnostic assessments. In the y-axis, all states belonging to the set of viral states ℋ, ℳ and are indicated. denotes an undetectable total viral load, i.e. this is the set of states for which condition n_C(M) ≤ ℓ for all possible virus mutants M holds ($[\le \ell ,\le \ell ,\le \ell ,\le \ell ]$). Likewise, ℋ denotes a high total viral load, i.e. refers to all states for which for at least one viral strain M, n_C(M) > m. The remaining viral states belong to ℳ. Only the initial part of the trajectory is presented (day 0–800 after treatment initiation) and details of transitions to each state are labeled for clarity. In panel B, the black, red, magenta and blue lines represent the probabilities of states ✠ (patient death), ℋ, ℳ and after application of the pro-active strategy. Note, that for the pro-active strategy, the x-axis is logarithmically scaled.

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

The optimal pro-active strategy depends on the initial probability state of the patient p₀. We assumed that the patient is treatment naive and has high virus copy numbers, i.e. $p[h,0,0,0](t_0)=1$ and $p\left[x\right](t_0)=0$ for $x\in 𝒮\backslash [h,0,0,0]$. For this scenario, it is suggested to start with treatment line a₁ and to switch to a₂ after 14 days, which is then maintained. The trajectories of the patient probability states are depicted in Fig 4B. For the ease of interpretation, we illustrate only the sets of viral states 𝓛, ℳ, ℋ and patient death ✠. 𝓛 denotes an undetectable total viral load. Translated to our model, this is the set of states for which condition n_C(M) ≤ ℓ for all possible virus mutants M holds, i.e. the current state has to fulfill $[\le \ell ,\le \ell ,\le \ell ,\le \ell ]$ to belong to this set. Likewise ℋ denotes a high total viral load, i.e. refers to all states for which for at least one viral strain M, n_C(M) > m is fulfilled. The remaining viral states belong to ℳ. One can nicely see that after approximately 260 days, maximum viral suppression may be achieved in the sense that the probability to have undetectable virus load (𝓛) is maximal (64.19%), while the patient may have intermediate viral loads ℳ with 15.57% probability and high viral loads ℋ with only 14.40% probability (the probability of death is 5.84%). After this time, it becomes more likely to fail treatment, as indicated by an increase in states ℳ and ℋ relative to 𝓛. We also assessed the sensitivity of the optimal pro-active strategy to variations in parameter values and found it to be fairly insensitive to parameter perturbations, see S3 Text. For comparison, we also show the dynamics for the case when no treatment switches were conducted in S4 Text in relation to the optimal pro-active strategy.

Cost of Strategy

In our model, the cost incurred by a treatment strategy can be divided into two types: The direct costs, which include treatment- and diagnostic costs, and indirect costs incurred by the virologic/health status of a patient (state costs). The pro-active strategy does not comprise diagnostic tests, whereas the protocol for the current standard of care (S.O.C.), as well as the protocol used in the HPTN052 [16], which we simulate for comparison, require viral load measurements. Currently, the expensive resistance tests are not part of the protocol for the standard of care, nor were they used for treatment decisions in HPTN052. The protocol for S.O.C. recommends changing treatment, if viral load (which is measured at month 6 and then every 12 months) is detectable and confirmed in a follow up testing after 2 months. The protocol for the HPTN052 trial recommends changing treatment, if two consecutive viral load measurements were greater than 1000 copies/mL, 16 weeks after treatment initiation. Viral load was measured at week 2, at month 1, 2, 3 after treatment initiation and then every 3 month. The cost of virologic testing is roughly 30 US$ per test [57, 58]. In contrast to S.O.C. and HPTN052, the diagnostic-guided strategy requires drug resistance testing. We set the cost of the diagnostics for the diagnostic-guided strategy to 200 US$ per test, in line with the recent literature [57, 59].

Download:

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

Table 3. Expected discounted costs for an infinite time horizon.

https://doi.org/10.1371/journal.pcbi.1004200.t003

Table 3 displays the expected discounted costs for an infinite time horizon for different strategies and highlights the direct- and indirect costs of each strategy, respectively. This comparison shows that the pro-active strategy performs best (83,819 US$), followed closely by the diagnostic-guided strategy (83,858 US$), the HPTN052 protocol (84,600 US$) and then by the standard of care (85,641 US$). The total expected discounted costs for the pro-active- and the diagnostic-guided strategy are 2% less than that of the standard of care. The state costs (indirect cost related to patient-well being) are the major determinant of the total cost, making up roughly 98%, 97%, 97% and 93% of total cost for the S.O.C., the HPTN052 protocol, the pro-active—and the diagnostic-guided strategy respectively. In terms of state costs, the diagnostic-guided strategy performs best.

The direct costs (treatment and diagnostic costs) are highest for the diagnostic-guided strategy (5,539 US$) followed by the pro-active strategy (2,772 US$), the HPTN052 protocol (2,390 US$) and the standard of care (1,871 US$). The direct costs make up only 2%, 3%, 3% and 7% of the total costs for S.O.C., the HPTN052 protocol, the pro-active and the diagnostic-guided strategy respectively. The direct costs of the pro-active and the diagnostic-guided strategy are roughly 48% and 196% more than that of S.O.C.

Patient Survival

Clearly, the primary goal of any treatment strategy is to improve and prolong the life expectancy of the treated individual. We therefore compare the distinct treatment strategies in terms of patient survival. For that purpose, we define the following term: which denotes the probability of death ✠ at time s given that the patient was treated according to treatment strategy stg. Given two distinct strategies; stg and a reference treatment strategy stg_ref, the term $T_{0\to t}^{+}(\mathrm{\text{stg}},{\mathrm{\text{stg}}}_{\mathrm{\text{ref}}})$ refers to the expected years of life gained (life prolongation) when the treatment strategy stg is used, relative to the reference treatment stg_ref at time t after initiation of treatment: (32)

In other words, given a patient is treated with stg and another patient is treated with stg_ref for time t, the terms $T_{0\to t}^{+}(\mathrm{\text{stg}},{\mathrm{\text{stg}}}_{\mathrm{\text{ref}}})$ refers to the expected time that a patient treated with stg will live longer than the patient treated with stg_ref.

We compared all strategies with the following reference strategies stg_ref: i) no medical intervention, ii) the standard of care treatment, iii) treatment according to the HPTN052 protocol and iv) the diagnostic-guided strategy. Fig 5A and 5D show the trajectories of expected life prolongation by different strategies in relation to i)-iv). Table 4 displays the expected life-years gained after 1 -, 2 -, 5 -, 8 -, 12—and 13.7 years of treatment respectively, where we additionally show the expected life prolongation in relation to the uninfected state.

Download:

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

Fig 5. Relative expected life prolongation [years] for different treatment strategies.

The purple solid lines, green dashed lines, blue dash-dotted lines and black dots represent the pro-active strategy-, the diagnostic-guided strategy, the current standard of care and the HPTN052 protocol respectively. The thin black line denotes the line of unity (no improvement/worsening). Panels A-D show the expected life-time prolongation for the distinct treatment strategies in relation to the no treatment, standard of care, the HPTN052 protocol and the optimal diagnostic-guided strategy respectively.

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

Download:

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

Table 4. Expected relative life-time gained using different strategies.

https://doi.org/10.1371/journal.pcbi.1004200.t004

The first five rows of Table 4 show the expected loss-of-life-time of an HIV infected person treated with distinct strategies in relation to an HIV uninfected person. After 13.7 years, an HIV patient receiving no treatment lives on average 6.2 years less than a healthy person. An HIV patient receiving treatment according to S.O.C., the pro-active strategy, the diagnostic-guided strategy or according to the HPTN052 protocol lives on average 3, 2.66, 2.3 and 2.82 years less than a healthy person. Fig 5A shows that all treatment strategies are better than receiving no treatment at all and prolong the life of an HIV patient by at least 3.2 years in relation. Fig 5B shows that the diagnostic-guided, pro-active strategy and the HPTN052 protocol are better at increasing patient survival than the standard of care. Further, Fig 5C shows that the optimal strategies are slightly better than the HPTN052 protocol and Fig 5D shows that the pro-active strategy and the HPTN052 protocol are slightly worse than the diagnostic-guided strategy. Table 4 shows that during the initial 2–3 years of treatment, there is almost no difference between the diagnostic-guided and the pro-active strategy with regard to patient survival. After 13.7 years of treatment, the difference between the two optimal strategies is less than 5 month (0.358 years).

Expected Reduction in Secondary Cases

Besides the primary goal of improving the life of the HIV patient, ‘treatment for prevention’ has gained interest in recent years. ‘Treatment for prevention’ strategies reduce onward transmission of the virus by reducing the infectiousness of HIV positive individuals. In order to measure the efficacy of the treatment strategies in preventing HIV-1 transmission, we estimated the incidence rate per 100 person-years associated with each HIV lumped state (ℓ, m, h) from a meta-analysis by Attia et al [14] (see S5 Text). The meta-analysis summarizes the outcome of 11 clinical studies on HIV-1 transmission in heterosexual sero-discordant couples, primarily from Africa.

For a strategy stg applied for a time t, the following equation gives a measure of the expected number of secondary cases/transmissions per survivor (33) where 𝕀ℝ(x) is the incidence rate per 100 person-years for a state x in our virus dynamics model, as explained in S5 Text and given in Table 1. Given two strategies, stg₁ and stg_ref, the percentage of potential infections prevented by strategy stg₁ in comparison to the reference strategy stg_ref is given by the quotient: (34)

We computed the expected reduction of secondary cases for different strategies taking either no treatment or the current standard of care as the reference strategy. In comparison to no treatment, the maximal reduction of secondary cases for the pro-active -, the diagnostic-guided strategy, the HPTN052 protocol and S.O.C. are achieved roughly 1.5–3 years after treatment initiation with values of 86%, 87%, 82% and 79% respectively, see Fig 6A. The relative reduction of secondary cases per survivor for the diagnostic-guided and the pro-active strategy are very similar, with an increase for the first 2 years, followed by a slow decline (see Fig 6A and Table 5). The relative reduction of secondary cases per survivor for the HPTN052 protocol is slightly better than that of S.O.C, with a tendency to decline over time, see Table 5. Note, that the computed relative reduction of secondary cases with the HPTN052 protocol was 82% (Table 5), which is slightly lower than the reported relative reduction of transmission events in the actual HPTN052 study [16] (reduction of 96% of linked and 89% of total transmission events). We have discussed reasons for this apparent under-prediction later in the manuscript. The difference between the optimal strategies (diagnostic-guided and the pro-active strategy) and S.O.C. becomes evident, when looking at the relative risk reduction by the optimal treatment strategies in relation to S.O.C. in Fig 6B. The reduction in secondary cases per survivor by the optimal strategies in comparison to S.O.C. is highest at the beginning and then slowly decreases over time.

Download:

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

Fig 6. Comparison of the relative reduction of secondary cases per survivor.

The purple solid, green dashed, blue dash-dotted lines and black dots represent the expected relative reduction of secondary cases per survivor by the pro-active-, the diagnostic-guided strategy, S.O.C. and in the HPTN052 protocol. In panel A, the reference strategy is no treatment and in panel B it is S.O.C.

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

Download:

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

Table 5. Expected relative reduction of secondary cases per survivor using different treatment strategies after different treatment durations.

https://doi.org/10.1371/journal.pcbi.1004200.t005

Discussion

The main aim of this work was to develop a rigorous mathematical framework that allows to compare different treatment paradigms in terms of monetary costs, treatment benefit and efficacy for ‘treatment for prevention’. It was previously stated [60], that the durability of ‘treatment for prevention’ should be assessed. Our simulations over a long time horizon (up to 5000 days/13.7 years) indicate that the effect of ‘treatment for prevention’ is significant and remains relatively stable beyond the time horizon typically assessed in clinical studies, see Fig 6A and Table 5, and that it may even be improved. We estimated that a standard of care therapy in e.g. South Africa can achieve a 66–79% reduction of HIV-1 onward transmission, in comparison to delivering no treatment. We also implemented the HPTN052 protocol, as stated in [16] and predicted that it would achieve up to 82% reduction of HIV-1 transmission, being more effective than the current standard of care, as shown in Fig 6B.

Statistical assessment of the actual HPTN052 trial [16] yielded estimates for the relative reduction of transmission of 96% for linked transmission and 89% for any transmission. Our simulated HPTN052 protocol yielded a 82% reduction of onwards transmission, which is within the confidence range of the reported estimates (CI: 73–99% for linked transmission and CI: 68–96% for any transmission) [16]. Note, that only one linked transmission event (1/1585 person-years) was observed in the early therapy arm of HPTN052 [16], giving rise to the statistical uncertainty in the reported estimate. Nevertheless, our simulations may under-predict the efficacy of HPTN052 due to several factors:

1. The reported treatment efficacy in the HPTN052 study was higher than predicted by our model: Virologic failure was only observed in 5% of participants in the early-therapy group of HPTN052, possibly explaining the difference between the outcome of the simulation vs. the clinical trial.
2. Despite only 5% failing to suppress the virus in the HPTN052 study, 66% initiated a second line therapy [16], meaning that a significant proportion of patients switched treatment before/without virologic failure. In our simulations of the HPTN052 protocol, patients only switched treatment when they showed signs of virologic failure. However, one may speculate that these treatment switches before/without virologic failure may have an impact on the long-term viral suppression that could be similar to a pro-active treatment switch.
3. The primary measurable endpoint of the HPTN052 study was the infection of the sero-discordant partner. Onward transmissions to other individuals could not be quantified for obvious reasons.

While a number of trials are now underway to confirm the results of HPTN052, see e.g. [61, 62], our in silico approach specifically addresses the need for an improved treatment strategy, particularly taking affordability into account, which suggests strategies that are suitable for scaling up.

Our work may indicate that if ‘treatment for prevention’ is scaled up and implemented using the current standard of care treatment strategy, its efficacy may not be as high as expected from HPTN052. Unlike in HPTN052, where monitoring of treatment success (viral suppression) and timely execution of treatment changes were realized, in resource-constrained countries close patient monitoring is currently not implemented in a routine setting and may be difficult to realize due to infrastructural and economic requirements.

Two alternative strategies for the immediate initiation of therapy were assessed in our work that take into account the mentioned limitations. Both suggested strategies (the diagnostic-guided strategy and the pro-active strategy) yielded better results in our simulations in terms of the reduction of onward transmission (see Table 5) at a lower price (Table 3). Both optimal strategies could yield a 72–87% reduction in HIV onward transmission in comparison to no treatment, see Fig 6A and Table 5. In comparison to the standard of care, we estimated that the diagnostic-guided strategy and the pro-active strategy could yield another 33–38% reduction of onward transmission after 2 years of treatment, but the advantages of the diagnostic-guided strategy and the pro-active strategy over the standard of care slowly declined over time, see Fig 6B. This indicates that both optimal strategies have a particular strength in reducing early transmissions (shortly after treatment initiation) in comparison to the current standard of care. This may be of particular utility, if transmission occurred primarily during early infection [63, 64]. In our work, we did not take behavioral factors into account, which would lead to a time-dependency of the infection rate. Rather, we assumed that the infection rate 𝕀ℝ(x) is constant over time, but dependent on the total virus load as reported earlier [14, 15, 19–21]. If transmission would primarily take place during an early infection, the advantages of the diagnostic-guided strategy and the pro-active strategy over the standard of care would be even more pronounced than indicated in Fig 6B.

The optimal diagnostic-guided strategy suggested patient-specific diagnostics, i.e. dependent on the patient’s virologic status (see S1 Table), unlike fixed intervals as in S.O.C, or the protocol stated in [16]. In summary, the optimal diagnostic-guided strategy suggests to take frequent diagnostics (≈ every month) if the patient is infected with a high or medium number of treatment-susceptible virus and less frequent (≈ every 5 month) diagnostics if the patient is infected with a ℓow/undetectable number of virus. No diagnostics are recommended for the remaining virologic states. Altogether, a very sparse diagnostic schedule for individual patients is suggested. Previous work [25] indicated that price reductions for the diagnostic tests would yield a better patient-outcome, which indicates that available drugs may not be utilized optimally in resource-poor settings, because diagnostics are currently too expensive. Of note is the fact that albeit treatment being available at very low expense (due to negotiations by the Clinton Foundation [65]), diagnostics may not be.

Furthermore, we suspected that allowing treatment change only after diagnostic confirmation of treatment failure (i.e. some time after drug resistance has occurred) may limit future treatment options [34]. Since the optimal diagnostic-guided strategy suggested rare diagnostics, and because it only allows to change treatment after resistance is detectable, we evaluated pro-active switching strategies. Note, that pro-active treatment switching strategies tested in the clinic increased virologic suppression and lowered rates of drug resistance emergence in HIV-1, when compared to conventional strategies [66, 67]. Similar strategies are also used against bacterial infections and cancers.

The computed pro-active strategy suggests a single treatment change without prior ascertainment of the viral status within a treated patient. Surprisingly, this strategy could yield comparable outcomes in terms of monetary costs, patient health and reduction of onward transmission, see Tables 3–5 and Figs 5 and 6. Our work thus indicates that pro-active strategies, may be as effective as diagnostically-driven ones, when diagnostics are expensive or inaccessible. Note, that unlike other optimal control approaches, i.e. [28] that suggest infinitely fast switching between regimens to mitigate drug resistance emergence, our predicted pro-active strategy actually only recommends a single treatment change, which is clinically more realistic. We also analyzed the sensitivity of the pro-active strategies with respect to the timing of this switch (see Fig 7). The graphic illustrates, that the switch is optimal after 14 days, however the difference in the performance measure is marginal, as long as the treatment switch is performed before ≈ 30 days (1 month) after treatment initiation.

Download:

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

Fig 7. Sensitivity of the pro-active strategy with respect to timing of the treatment switch.

The purple solid line represents the total expected discounted costs for the pro-active strategy with respect to the switching time (shown on the x-axis). The horizontal green dashed, blue dash-dotted lines and black dots represent the total expected discounted costs for the optimal diagnostic-guided strategy, S.O.C. and in the HPTN052 protocol. The inset shows a zoom into the first 30 days after treatment initiation.

https://doi.org/10.1371/journal.pcbi.1004200.g007

Obviously, pragmatic and clinical considerations need to be taken into account to translate our results into practice. Also, several assumptions have been made, which require careful evaluation. For example, we used a very coarse-grained model of HIV within-host dynamics, which was required to enable the numerical computation of optimal controls, particularly for the closed-loop system employed in the diagnostic-guided strategy. Most models of viral dynamics, e.g. [33, 68, 69], were developed to accurately predict short-term viral dynamics after drug application and are unable to predict virologic failure after long time intervals, in contrast to our coarse-grained model, which was developed and parameterized in order to predict short-term viral dynamics as well as virologic failure after very long time-intervals, see e.g. Figs 2 and 3. It is therefore more suitable than existing models in estimating the long-term response to antiviral treatment. However, in the future we will focus on developing more elaborated HIV-models and on algorithms to solve the control problem for the chemical master equation directly, without state-space lumping. Note, that other computationally efficient numerical approaches, such as model predictive control [30], could be used to approximate the optimal treatment strategies. However, there is no guarantee that the computed control using these approaches will be optimal.

In our approach, we modeled treatment change as a switched system, which neglects the pharmacokinetics of the distinct drugs [10, 70–72] and may only indirectly reflect drug adherence in an average population (drug efficacy η is a constant term in our model). Neglecting pharmacokinetics may, however, be a justifiable step in this modeling exercise, because of the considered time-scales (on the order of years), and also because optimizing e.g. drug adherence was not an objective of this study. However, if the main interest is for example in optimizing the switch between two treatment lines by optimal dosing in order to prevent time frames of insufficient viral suppression or drug over-exposure, or to include patient-specific or time-dependent drug adherence, we advise to consider a different control system, for example [73]. Within such a framework, monitoring (e.g. viral load assessment) may also be incorporated as a tool to assess individual drug adherence and to allocate resources to improve it.

We did not explicitly consider costs related to contraindications caused by the treatment. For example, the second treatment line a₂ may be less tolerable. Mathematically, this can be modeled in terms of increased treatment costs for a₂, in comparison to the first treatment line a₁. In order to test the sensitivity of the optimal pro-active strategy to this parameter, we conducted the necessary computations in S3 Text and found that the computed strategy was fairly insensitive to changes in treatment costs. This may indicate that the benefits of the treatment switch outweigh these potential shortcomings.

Also, we did not include screening costs or the costs of the initial virologic assessment, thus our calculations refer to the public health costs that accrue from the start of HIV treatment. These costs will, however, only enter as a constant to each of the tested strategies and will not change the results beyond the addition of this constant to the values stated in Table 3. Additional costs (personnel, infrastructure, transportation) may come along with diagnostic assessments. It is likely that hidden costs for diagnostics are substantial. With a higher cost of diagnosis, the pro-active strategy may outperform the diagnostic-guided strategy, which may suggest an even less frequent diagnostic schedule, supporting our claim that current diagnostics may be too expensive to be appropriately used.

We used the price of a drug resistance test (k_dia ≈ 200 US$ [57, 59]) to account for diagnostics in the diagnostic-guided strategy. This had the following reason: Current guidelines recommend to measure the total virus load [13] and to switch treatment, if, based on this partial information, virologic failure is anticipated. As reported by others [57], this may lead to unnecessary treatment switches. In contrast, a resistance test directly informs the physician about the necessity of treatment change. Mathematically, partial information, i.e. the total virus load, would lead to a distinct control framework, namely Partially Observable Markov Decision Processes (POMDP) [74], which are extremely challenging to solve, particularly for larger models like the one used herein (Fig 1). In POMDPs, partial information may be mapped into a ‘believed’ full virologic status, for example observing a high total virus load may be due to some resistance development, e.g. the viral state $[\ell ,h,0,0]$ with some probability. However, it is hard for us formalize the physicians intuition (i.e. the relation between observation, belief and truth) regarding this ‘mapping’ of partial measurements to viral states x.

As a primary outcome of our modeling exercise, we estimated the expected relative number of secondary infections prevented (Table 5 and Fig 6); -unlike many other approaches (summarized in [18]), which take the absolute number of secondary cases into account. Estimating absolute numbers of secondary cases would require to model complex behaviors, i.e. sexual relationships, etc. over time, for which we do not have data for validation, nor was it the main focus of the current work. This also prevents us from predicting the course of the epidemic or deriving its reproductive number R₀ in relation to distinct treatment strategies. However, the primary aim of this study was to compare the efficacy of different treatment strategies, which is nicely quantified in terms of the expected relative number of secondary infections prevented. Note, that this relative estimate requires no assumptions on the underlying transmission dynamics, except that it assumes that these dynamics are similar for a tested strategy versus its comparator.

In addition to insights in HIV ‘treatment for prevention’ strategies, the developed mathematical/control theoretic framework may already be applicable to many medical phenomena. Further developments may improve its applicability to even more complex processes, which can be accurately described by intrinsically stochastic dynamics. For example, the open-loop optimal control approach (used to determine the optimal pro-active strategy) may be turned into a closed-loop system, if diagnostics are taken from time-to-time to determine the viral state of a patient, i.e. p[x](t_j). Also, the closed-loop system that requires state determination (the diagnostic-guided strategy) may be combined with the open-loop system in order to allow for pro-active treatment changes in between diagnostic assessments.

Supporting Information

Acknowledgments

We thank Maureen Smith for proof-reading the manuscript.