The impact of infected T lymphocyte burst rate and viral shedding rate on optimal treatment scheduling in a human immunodeﬁciency virus infection

—We consider a mathematical model of human immunodeﬁciency virus (HIV) infection dynamics of T lymphocyte (T cell), infected T cell, and viral populations under reverse transcriptase inhibitor (RTI) and protease inhibitor (PI) treatment. Existence, uniqueness, and characterization of optimal treatment proﬁles which minimize total amount of drug used, viral, and infected T cell populations, while maximizing levels of T cells are determined analytically. Numerical optimal control experiments are also performed to illustrate how burst rate of infected T cells and shedding rate of virions impact optimal treatment proﬁles. Finally, a sensitivity analysis is performed to detect how model input parameters contribute to output variance.


I. INTRODUCTION
HIV, also known as the human immunodeficiency virus, is a virus which impairs the immune system by entering vital cells (e.g. dendritic cells, microglial cells, CD4+ T cells, and macrophages) through CD4 and coreceptors on the cell membrane and destroying them [Cunningham et al. 2010;Cenker et al. 2017]. Once the CD4+ levels drop to a critical level, cell-mediated immunity is lost and viral load increases, which leads to the development of acquired immunodeficiency syndrome (AIDS), the final stage of HIV infection. When this occurs, due to low immunity, the patient's body becomes a breeding ground for opportunistic infection and AIDS defining cancers such as Kaposi sarcoma, non-Hodgkin lymphoma and cervical cancers. Without treatment, average survival time after HIV infection is 9-11 years [UNAIDS 2007].
Pioneering research done by scientists in the 1990s determined that direct cell-to-cell contact is an important, and perhaps the predominant, contributor to the propagation of HIV within a host. Specifically, HIV spreads by entering immune cells and subverting intercellular communication to enable their own viral replication through structures such as tunneling nanotubes and filopodia. Immunological synapses and phagocytosis of infected cells are also hijacked by HIV and used as gateways to infect target cells. Furthermore, HIV is able to elicit fusion between infected donor and target cells, forming horrific infected cell clusters with a very high capacity for viral reproduction and survival. There are many mechanisms by which cell-cell transfer of HIV occurs, some of which are summarized in Figure 1 from [Bracq et al. 2018]. For more details on how these mechanisms work, refer to [Bracq et al. 2018].
The main treatments for HIV are a class of drugs called antiretrovirals (ART). Often, these drugs are given in combination to suppress viral replication and reduce plasma HIV viral load in a treatment called highly active antiretroviral therapy (HAART) [Autran et al. 1997;Komanduri et al. 1998;Lederman et al. 1998]. Two common such ARTs often used in HAART are reverse transcriptase inhibitors and protease inhibitors [Arts et al. 2012].
Reverse transcriptase inhibitors work to terminate transcription, thus inhibiting viral replication. They do this in one of two ways: nucleoside/nucleotide reverse transcriptase inhibitors incorporate into the nascent DNA strand during reverse transcription, whereas non-nucleoside reverse transcriptase inhibitors bind to non-catalytic enzyme sites and is usually mediated through steric hindrance that impedes structural changes in HIV reverse transcriptase [Gulnik et al. 1995;Gotte et al. 2000]. Protease inhibitors are usually substrate-based inhibitors which are designed specifically against the viral protease based on its crystal structure. They act on this viral protease, inhibiting maturation of new viral particles. By doing this, they attack the already formed HIV before the next cycle of infection begins [Michaud et al. 2012]. Figure 2, reproduced from [Michaud et al. 2012 Page 2 of 12 A. Bukkuri, The impact of infected T lymphocyte burst rate and viral shedding rate on optimal ... Currently, physicians primarily use three markers to assess the stage of the disease and create treatment strategies: CD4 T cell count, viral load, and drug resistance [PENTA 2004]. Through optimal control and sensitivity analyses, we propose the inclusion of two more criteria into this assessment: viral shedding rate and infected T cell burst rate. Medically, it is possible to measure the rate of viral shedding using polymerase chain reaction (PCR) from collected samples as done in [Tronstein et al. 2011]. Burst rate of infected T cells can be measured using standard burst rate analysis techniques, such as through viral inhibitors, washouts of infected cells [Dimitrov et al. 1993;Eckstein et al. 2001], quantitative image analyses with in situ hybridization, quantitative competitive RT-PCR of bulk tissue or single cells at limiting dilution [Chun et al. 1997;Haase et al. 1996;Hockett et al. 1999;Chen et al. 2007;Hyman et al. 2009].
In this paper, we first summarize a model of HIV infection with RTI and PI treatment developed by [Mobisa et al. 2018]. Then, in section 3, a detailed optimal control analysis will be performed, including existence, uniqueness, and analytical characterization results. In section 4, parameter estimates are provided and several numerical simulations were performed, specifically analyzing the impact of changes to the burst rate of infected cells and shedding rate of virions on optimal treatment profiles. In section 5, a sensitivity analysis is performed to quantify the amount of impact each parameter value has on healthy T cell, infected T cell, and viral population dynamics.

II. MODEL DESCRIPTION
The model for HIV infection dynamics among T cells (T), infected T cells (U), and free virus particles (V) under RTI (u) and PI (v) treatment is taken from a model created by [Mobisa et al. 2018] and is reproduced below: In this model, the T cells are produced at rate r and die naturally at rate µ; the growth of the T cells is bounded with the enforcement of a carrying capacity. The healthy T cells become infected by the free virus and actively infected T cells via the mass infection terms βV T and αT U , respectively. The infected T cells die naturally at a rate of k. Furthermore, these cells produce free viruses, as captured by the first term in the viral equation, which are then removed from circulation at a rate c per virus. Note that the viral decline depends on the amount and efficacy of the RTI and PI treatments.

III. OPTIMAL CONTROL ANALYSIS
To determine best treatments strategies for HIV infection, we formulate our problem as an optimal control one. We aim to maximize the number of healthy T cells, while minimizing the number of infected T cells, free viruses, and drugs u and v used over a fixed therapy horizon [0,T].
Our control class are piecewise continuous functions defined for all t such that 0 ≤ u(t) ≤ 1 where u(t) = 1 represents maximal RTI treatment and u(t) = 0 represents no RTI treatment. Similarly, v(t) = 0 represents no PI treatment and v(t) = 1 represents maximal PI treatment. We can then describe the class of admissible controls as Now, we define our objective functional and optimal control problem. For a fixed therapy interval [0,T], maximize the objective functional subject to the above ODE dynamics and initial conditions of T(0) = 1000, U(0) = 1, V(0) = 0.1. Note that the objective functional is quadratic, instead of linear, in u and v to ensure that the resulting optimal controls can take on a continuum of values from 0 to 1, rather than being binarily confined to either 0 or 1.

A. Existence of Optimal Control
Using the theory developed by [Fleming et al. 1975], we can determine the existence of an optimal control for our state system. Specifically, boundedness of solutions of the system for finite time is required for existence and uniqueness of an optimal control. Using the technique of supersolutions and keeping in mind that T ≤ T max , upper bounds on the solutions of the state system are found. Consider the following system: Note that the supersolutionsT ,Ū , andV of this system are bounded on a finite time interval. We can also re-write this system in matrix form as follows: Since this is a linear system in finite time with bounded coefficients, the supersolutions T ,Ū , andV are uniformly bounded. We can now prove existence of an optimal control as done in [Bukkuri 2019].
Theorem 1: There exists optimal controls u and v that maximizes the objective functional J(u,v) if the following conditions are met: 1) The class of all initial conditions with controls u and v such that u and v are Lebesgue integrable functions on [0,T] with values in the admissible control set along with each state equation being satisfied is not empty 2) The admissible control set is closed and convex 3) The right hand side of the state system is continuous, is bounded above by a sum of the bounded control and the state, and can be written as a linear function of u and v with coefficients depending on time and the state variables 4) The integrand of the functional is concave on the admissible control set and is bounded above by c 3 −c 2 |u| θ −c 1 |v| ψ , where c 2 , c 1 > 0, and θ, ψ > 1.
Proof: First, from an existence result in [Lukes et al. 1982], since our state system has bounded coefficients and any solutions are bounded on the finite time interval, we obtain the existence of the solution of the state system. Second, the admissible control set is closed and convex, by definition. For the third condition, the right hand side of the state system is continuous since each term with a denominator is nonzero. Moreover, the system is bilinear in the controls and can be rewritten as where X = (T, U, V ) and γ is a vector-valued function of X. Since the solutions are bounded, we have where C 1 depends on the coefficients of the system. Also, note that the integrand of J(u,v) is concave on the admissible control set. The existence of optimal control follows from the fact

B. Characterization of Optimal Control
Now, we will characterize the optimal control pair (u,v) using a version of the Pontryagin Maximum Principle (PMP), as done in [Bukkuri 2019;Bukkuri 2020]. First, let's define the Lagrangian associated with J(u,v): Here, the ν and δ terms have been added in as penalties for non-optimal dosing patterns, i.e.
Theorem 2: Given optimal controls u * and v * and solutions of the corresponding state system, there exist adjoint variables λ i for i = 1, 2, 3 such that where λ i (T ) = 0 for i = 1, 2, 3 by the PMP transversality condition. Moreover, u * is given by: while v * is similarly given by:  Proof: To maximize the Lagrangian (with respect to the optimal control pair variables), we differentiate L with respect to u and v. From this, we get: Thus, the representation of u Since 0 ≤ u ≤ 1 and 0 ≤ v ≤ 1, the explicit control profiles given in equations 5 and 6 are determined. Moreover, since both state and adjoint solutions are L ∞ -bounded, the right side of the adjoint and state equations are Lipschitz for those solutions. This ensures that the solutions to the optimality system are unique (and thus the optimal controls are unique), assuming the final time is not very large. A rigorous proof of such an argument can be found in [Burden et al. 2004] and [Fister et al. 1998].

IV. NUMERICAL SIMULATIONS
We now aim to numerically implement our optimal control solutions to visualize optimal drug treatment strategies. The optimality system can be thought of as a two-point boundary value problem, which was solved using a fourth-order iterative Runge-Kutta scheme, as done in [Jung et al. 2002]. In this scheme, a forward sweep of the state equations with initial guesses for u and v was performed. Then, a backward calculation using the adjoint equation and an update of the controls was made as done by [Duda 1997] and [Deininger et al. 2003]. This was iteratively performed until convergence was obtained.
The parameter values in Table I represent standard HIV infection values, were obtained from [Mobisa et al. 2018], and were used in the following optimal control simulations.
First, we run a control simulation. In this simulation, we let w = 1000 and c = 2.4. As explained in [Mobisa et al. 2018], w varies greatly and we shall later explore how optimal treatment protocols change when the burst rate is much lower or higher. The value of 2.4 for c is based on viral infectivity assays which found that HIV-1 strains IIIB and RFII which lost half of their infectivity in 4-6 h at 37 • C [Layne et al. 1989]. However, it's clear that different strains of HIV in different individuals and at different temperatures may lead to different viral shedding rates. We shall soon see the impact of lower viral shedding rates on optimal protocols. Figure 3 is a depiction of optimal RTI and PI treatment for our control case. Note that, in the simulations below, 100 time steps is equivalent to one day; thus, the following simulations were run for over 600 days. From this, it's clear to see that the optimal treatment calls for full intensity of both RTI and PI treatment for about 290 days. Now, let's take a look at what happens when we change the burst rate. In Figure 4, the upper image represents w = 1 and the lower image is for w = 100, 000.
In the case of the low burst rate, we notice that the optimal profiles does not change much-the only change is that it takes slightly longer to reach the maximal dosage for both RTI and PI treatment. However, it is clear that the high burst rate has a significantly different profile. Specifically, we see the emergence of bang-bang controls for both RTI and PI treatment; moreover, we notice that the RTI treatment is given for slightly longer (≈ 20 days) in the beginning than the PI treatment. Both treatments then have a drug holiday until ≈ day 210, after which both drugs are given at full potency for ≈ 80 days, before all treatment is stopped.
Let's now consider the impact of different shedding rates, for the following values for c: 0.1, 0.2, 0.3, and 5. In Figure 5, the upper three images represent c values of 0.1, 0.2, and 0.3, respectively, while the bottom image represents a c value of 5. In these simulations, we again see a departure from our control optimal treatment recommendations. First, consider the high shedding rate value. In this case, optimal treatment recommends full dosage treatment of RTI and PI for ≈ 130 days instead of ≈ 280 days, before ceasing all treatment. Now, consider the lower shedding rates. When c is 0.1, we notice that the RTI treatment is given at full dosage for ≈ 190 days, before switching to rapid bang-bang controls (between full and no treatment) until ≈ 300 days. A drug vacation is then given for ≈ 90 days, before giving the RTI Fig. 5. Impact of Shedding Rate of Virions (c) treatment at full potential for ≈ 20 more days, then ceasing treatment. For the PI treatment, we notice that the treatment is prescribed at full potency for the first ≈ 80 days, before it gradually decays to no treatment over the next ≈ 20 days. The treatment again picks up around day 210, and switches to rapid bang-bang controls like the RTI until about day 300. Then we note the presence of bang-bang controls (mostly remaining at no treatment, though) during the RTI drug holiday, before continuing full treatment at ≈ day 380 for around 30 days, then ceasing all treatment.
When the c value is 0.2, the optimal protocols change drastically again. The RTI treatment is prescribed at full potency for the first ≈ 145 days, before declining to 0 treatment over the next 5 days or so. A drug holiday is then taken for around the next 250 days, after which the treatment is gradually increased (over ≈ 90 days) until the end of the considered therapy horizon. The PI treatment is given at full potential for the first ≈ 90 days; it then rapidly declines and stays at no treatment for the rest of the 600 day interval. Finally, let's consider what happens when the shedding rate is 0.3. In this case, the optimal profile includes an RTI treatment which stays at full potency for the entire interval. The PI treatment is given at its maximum for ≈ 190 days, before rapidly declining to no treatment, and then gradually increasing back to full potency (over an ≈ 80 day period) after about 10 days, and staying at this maximum level for the rest of the therapy horizon.
Thus, from the above numerical experiments, we can clearly see that the standard optimal treatment protocol, which calls for full treatment of both RTI and PI for ≈ 290 days, can drastically change depending on each patient's specific infected T cell burst rate and viral shedding rate.
V. SENSITIVITY ANALYSIS To further assess the impact of parameter values on HIV dynamics we perform a sensitivity analysis on our ODE system using the Sobol-Martinez method. Other than u and v, which were given ranges of 0 to 1, all other parameters were given Page 9 of 12 a range of 10% lower to 10% higher than the values indicated in the table in Section 5. The Sobol-Martinez method uses variance decomposition techniques to measure the contributions of input parameters to output variance. The algorithm outlined by [Zhang et al. 2015] was implemented here.
The left panel of each image is the first order Sobol index: the contribution to the output variance by the single input alone. The right panel is the total Sobol index: the contribution to the output variance caused by a model input, including the first-order effects as well as all higher-order interactions.
Note that, in the Figure 6, uther is equivalent to u and vther is equivalent to v. Also, all parameters with sensitivity orders less than 0.01 were omitted from the figures. First, consider the healthy T cells. In this case, we see three main parameters which impact the T cell population: the burst rate of the infected T cells, the RTI treatment, and the viral shedding rate. The other parameters which had minor effects on T cell population dynamics are the natural death rate of healthy T cells and the viral infection rate by free virions. Thus, the most effective treatments, solely in terms of the CD4+ cell population, are the RTI treatment and those that target the burst rate of infected T cells and the shedding rate of virions. Considering the infected T cells, we note that almost all the parameters are effective in changing its population dynamics-thus the sensitivity analysis cannot be very enlightening for treatment planning. However, when we consider the viral graph, we see that the the shedding rate is most significant, followed by the burst rate of infected cells.
Thus, for the overall desired dynamics, we deem that further medical research into the development of drugs which specifically target the shedding rate of virions and the burst rate of infected cells is advised, as these are the most effective ways of minimizing the virus and infected T cell population, while maximizing the healthy T cell count. These results are in accordance with our nu-merical optimal control experiments, which show great changes in optimal treatment protocols when shedding and burst rate parameters are changed.

VI. CONCLUSION
In this paper, we performed optimal control and sensitivity analyses of an HIV model with reverse transcriptase inhibitor and protease inhibitor antiretroviral treatments. Proofs of existence, uniqueness, and analytical characterization of the optimal control profile were given. Numerical optimal control simulations were performed and it was found that burst rate of infected T cells and viral shedding rate have great impacts on optimal control profiles, often suggesting intricate bang-bang controls for the treatment. We hope that the results of this analysis will help physicians more effectively assess and treat HIV patients. Specifically, we hope physicians will measure and include viral shedding rates and infected T cell burst rates into treatment consideration. Moreover, we hope that increased pharmaceutical efforts go into developing drugs which specifically target burst and shedding rates. Finally, we hope that this inspires future work into creation of optimal treatment profiles for other similar viral infections.