Dynamics of Drug Resistance: Optimal Control of an Infectious Disease

Abstract

Antimicrobial resistance is a significant public health threat. In the U.S. alone, 2 million people are infected and 23,000 die each year from antibiotic resistant bacterial infections. In many cases, infections are resistant to all but a few remaining drugs. We examine the case where a single drug remains and solve for the optimal treatment policy for an SIS infectious disease model incorporating the effects of drug resistance. The problem is formulated as an optimal control problem with two continuous state variables, the disease prevalence and drug's "quality" (the fraction of infections that are drug-susceptible). The decision maker's objective is to minimize the discounted cost of the disease to society over an infinite horizon. We provide a new generalizable solution approach that allows us to thoroughly characterize the optimal treatment policy analytically. We prove that the optimal treatment policy is a bang-bang policy with a single switching time. The action/inaction regions can be described by a single boundary that is strictly increasing when viewed as a function of drug quality, indicating that when the disease transmission rate is constant, the policy of withholding treatment to preserve the drug for a potentially more serious future outbreak is not optimal. We show that the optimal value function and/or its derivatives are neither C ¹ nor Lipschitz continuous suggesting that numerical approaches to this family of dynamic infectious disease models may not be computationally stable. Furthermore, we demonstrate that relaxing the standard assumption of constant disease transmission rate can fundamentally change the shape of the action region, add a singular arc to the optimal control, and make preserving the drug for a serious outbreak optimal. In addition, we apply our framework to the case of antibiotic resistant gonorrhea.

Keywords: antimicrobial resistance; dynamic health policy; health care management; infectious disease models; optimal control.

Figure 1:

Action/Inaction boundary p_b along with the solution to the first-order optimality condition, i.e., $\overline{p}$, for (a) long-term focused planner (ρ₁ ≤ β − r) and (b) short-term focused planner (ρ₂ > β − r). Gray arrows indicate the vector field defined in the system of ODE’s (5)–(6) when W (t) = 1, t ≥ 0. In panel (a) where ρ₁ > β − r, observe that the infimum of the set of drug qualities for which it is economical to prescribe the drug for some p > 0 is zero (i.e., $\underset{\_}{q}=0$). In addition, since ${\overline{p}}^{\prime }(q)=-v_c^{\prime }(\overline{p}(q))/q{v}_c^{\prime \prime }(\overline{p}(q))$ converges to zero when $\overline{p}(q)$ is close to zero, $\overline{p}$ is tangent to the q axis. In contrast, in panel (b) where ρ₂ > β − r, it is $\underset{\_}{q}>0$. Similar to panel (a), since ρ₂ > 2β − 2r, the second derivative of v_c to diverges infinity when p approaches to zero. This causes ${\overline{p}}^{\prime }(\underset{\_}{q})=0$ leading to the tangency of $\overline{p}$ to the q-axis When ρ = ρ₃, however however, the second derivative of v_c is infinite at p = 0 and hence ${\overline{p}}^{\prime }(\underset{\_}{q})>0$. This case is shown in Fig. 9a in Appendix A.

Figure 2:

(a) The value function for an arbitrary switching time τ when the initial point (p, q) is (0.75, 0.8) and(b) the system trajectory when the switching time τ = 3. For any switching time τ < 3, the state of the system remains on the blue line for t ≤ τ after which moves vertically toward p_∞ = 0.75, i.e., Q(t) = Q_a(τ; 0.8) and P (t) = P_c(t − τ; P_a(τ; 0.75, 0.8)), t > τ. The two vertical dashed lines indicate the time when the state of the system enters $\mathcal{A}$ and exits $\mathcal{A}$. Panel (a) and (b) together illustrate that although initially v_s(0.75, 0.8; τ) is increasing in τ, it will eventually decrease as the system trajectory is “strongly entering”.

Figure 3:

Action/Inaction regions for (a) long-term focused planner (ρ₁ ≤ β − r) and (b) short-term focused planner (ρ₂ > β − r). In these two panels ${\mathcal{Q}}_n$ is indicated by gray the dark rectangle. The white area denotes the optimal inaction (continuation) region while the gray shading on the right-hand side of the solid line denotes the optimal action region. The optimal value function as a function of q and its first derivative when ρ₁ ≤ β − r and p = 0.68 (the horizontal dashed line in panel (a)) as well as when ρ₂ > β − r and p = 0.61 (the horizontal dashed line in panel (b)) are shown in panels (c) and (d), respectively. In panels (c) and (d), the left vertical axis represents v* and the right vertical axis represents its first derivative. In these panels a vertical dashed line has been added (at q = 0.64 and q = 0.72, respectively) drawing attention to the the non-differentiability of v* with respect to q at these points. Moreover, in panels (c) and (d) v*(0.68, q) and v*(0.61, q) are constant when q ≤ 0.64 and q ≤ 0.72 respectively as in the inaction (continuation) region ${\mathcal{C}}^{*}$, v*(p, q) = v_c(p).

Figure 4:

The optimal value function as a function of p and its first derivative for (a) long-term focused planner (ρ₁ ≤ β − r) when q = 0.64 (vertical line in Fig. 3a) and (b) short-term focused planner (ρ₂ > β − r) when q = 0.72 (vertical line in Fig. 3b). In each panel, the left vertical axis represents v* and the right vertical axis represents its first derivative. On each panel a vertical dashed line has been added (at p_b (0.64) = 0.68 and p_b (0.72) = 0.61, respectively) drawing attention to the the non-differentiability of v* with respect to p at (p_b(q); q), with $q\in {\mathcal{Q}}_n$. Note that when ρ₁ ≤ β − r, the derivative of v*(p, q) with respect to p is unbounded when p is equal to zero. In contrast, when ρ₂ > β − r, the derivative of v*(p, q) remains bounded when p = 0.

Figure 5:

(a) The action/inaction boundary for the gonorrhea case study. The shaded area indicates the region in which it is optimal to treat all patients and the unshaded area indicates the region in which it is optimal to treat no one (because the treatment is effectively obsolete due to high resistance levels). Consistent with our main analysis in section 3, it is never optimal to only treat a subset of infected individuals. The gray arrows depict the vector field for the gonorrhea case study when W (t) = 1, t ≥ 0, indicating the direction of the movement in the state space when treatment is offered to all infected individuals. The solid blue line represents the state trajectory for the initial point (p, q) = (0.27%, 99.6%) when treatment is provided to all patients (W (t) = 1, t ≥ 0), even after it is no longer optimal to provide treatment to anyone. The three dashed vertical lines correspond to q = 95%, q = 90%, and q = 6.4%. (b) Projected annualized gonorrhea total expenditures over an infinite time horizon for W (t) = 1_{t≤τ} as the function of τ. The three dashed vertical lines mark the times when Q_a(t; 99.6%) = 95%, Q_a(t; 99.6%) = 90%, and Q_a(t; 99.6%) = 6.4%.

Figure 6:

(a) Sensitivity of the optimal action/inaction region to the uncertainty in the model primitives for the gonorrhea case study. The shading represents the likelihood of a point being in the optimal action region with solid black representing 100% and solid white 0%. (b) The distribution of the stopping time corresponding to the optimal treatment policy.

Figure 7:

(a) The state space partition to action/inaction region when γ = −1.5 along with two optimal state trajectories. (b, c) The optimal state trajectory as a function of time for two initial starting points. (d) The optimal control as a function of time for two different starting points. To enable comparison, p-axis and q-axis in panels (b) and (c) are kept parallel to these axes in panel (a).

Figure 8:

Action/Inaction boundary p_b, denoted with a solid black line; the solution to the first-order optimality condition, i.e., $\overline{p}$, denoted by a dashed red line; optimal action region ${\mathcal{A}}^{*}$ marked by gray shading; optimal inaction (continuation) region ${\mathcal{C}}^{*}$ in white; set ${\mathcal{A}}^{*}\backslash \mathcal{A}$ marked by dark gray shading; and two hypothetical trajectories (P_a(t; p, q), Q_a(t; q)) staring from point a (p, q) in the interior of $\overline{{\mathcal{C}}^{*}}$

Figure 9:

(a) Action/Inaction boundary p_b along with the solution to the first-order optimality condition, i.e., $\overline{p}$, for 2(β − r) < ρ. (b) Illustration of step 2, part (a) and (b) of the proof of Thm. 2.

Figure 10:

The autonomous value function and its first derivative when (a) long-term focused planner (ρ₁ ≤ β − r), and (b) short-term focused planner (ρ₂ > β − r). In each panel, the left vertical axis represents v_c(p) and the right vertical axis represents its first derivative.

Figure 11:

The optimal state trajectory for three different starting points from four perspectives: (a) in the state space; (b) in the p-dimension as a function of time; (c) in the q-dimension as a function of time; and (d) of the optimal control as a function of time. Panel (a) also illustrates iso-curves of the optimal value function. To enable comparison, p-axis and q-axis in panels (b) and (c) are kept parallel to these axes in panel (a).