Optimal Controls of the Highly Active Antiretroviral Therapy

and Applied Analysis Volume 2020, Article ID 8107106, 23 pages https://doi.org/10.1155/2020/8107106 Abstract and Applied Analysis 2and Applied Analysis 2 is the oldest and the largest class of antiretroviral drugs that includes drugs such as Zidovudine, Didanosine, Stavudine, Lamivudine, Abacavir, Tenofovir, and combinations of these. e NRTI and NtRTI inhibit reverse transcription by being incorporated into the newly synthesized viral DNA strand as faulty nucleotides. (3) Nonnucleoside reverse transcriptase inhibitors, or NNRTI, (Nevirapine, Delavirdine, Efavirenz, Etravirine, Rilpivirine) also inhibit reverse transcriptase, but in a different way, by binding to an allosteric site of the enzyme. (4) Integrase inhibitors (Raltegravir, Elvitegravir, Dolutegravir) inhibit the enzyme integrase, which is responsible for integration of the viral DNA into the DNA of the infected cell. (5) Protease inhibitors (PI) inhibit the activity of protease, an enzyme that cut nascent proteins into specific pieces for the final assembly of new functioning virions. ere are ten approved protease inhibitors, namely: Saquinavir, Indinavir, Ritonavir, Nelfinavir, Amprenavir, Lopinavir/Ritonavir, Atazanavir, Fosamprenavir, Tipranavir, and Darunavir. (6) Maturation inhibitors inhibit the last step in gag processing where the viral capsid polyprotein is cleaved, blocking the conversion of the polyprotein into the mature capsid protein. As a result, the virions released have a defective core and are mainly noninfectious. e guidelines for use of the antiretroviral drugs are a source of controversy within the medical community. us, there is diversity of opinion regarding the stage of infection and the CD4+ T cells count (which are the primary target cells for HIV, and which depletion leads to the immune system collapse and development of AIDS) when the therapy should be started. It is generally accepted that the objectives of the treatment are to maximally suppress the viral load and to maintain CD4+ T cells at an acceptable level. Other factors that have to be also taken into consideration are: (i) a possibility of individual intolerance and adverse side-effects; (ii) a possibility of developing drug resistance for a patient, and (iii) a possibility of emergence of drug-resistant strains, which thereby can be passed onto others. Moreover, the cost of the therapy, which is high, should be also taken into consideration. To reduce the possibility of developing drug resistance and of the emergence of drug-resistant strains and to minimize the side effects, a combination of several antiretroviral drugs (typically, of three or four) is usually used. is approach is known as highly active antiretroviral therapy, or HAART. Also, it is generally accepted that, in order to reduce a possibility of emergence of drug resistant strains, antiretroviral therapy, once started, should never be stopped: as current US Guidelines states, “Patients initiating antiretroviral therapy should be willing and able to commit to lifelong treatment” [6]. e above-mentioned controversy in formulating the guidelines implies that some sort of optimization can be applied, and that the problem should be quantified and explored mathematically. Of course, at the present stage, it would be hardly possible to take into consideration details such as individual sensitivity of a specific patient to a specific drug. However, some common principles that can assist in developing better guidelines can be established. A natural approach for this sort of problems is to use the tools and methods of the classic optimal control theory. e problem of HIV control and the applicability of the optimal control theory to this problem attracted attention of mathematicians fairly early: seminal works by Butler et al. [7] and Kirschner et al. [8] initiated extensive literature in this direction. A mathematical techniques suggested in these publications proved to be successful and was subsequently applied to a variety of mathematical model of HIV (see, e.g., [9–22] and bibliography therein). ese results brought important insight into the problem. However, since this advance utilized essentially the same mathematical techniques, all these results also share some drawbacks, that we have to discuss before we proceed to our study. 1.2. Form of the Objective Functional. Starting from [7, 8], all the above-mentioned publication (and, in fact, the vast majority of the literature on the optimal control in mathematical biology and medicine) deal with optimal control problems of maximization of the cumulative (integral) concentration of healthy CD4+ T cells and simultaneous minimization of the cost of the treatment (represented by the weighted squares of the considered controls). Control problems with such functionals are mathematically convenient because they never deal with singular control phenomenon and can be solved numerically as a two-point boundary value problem. A setback of such objective functionals is that the optimal controls, obtained for such objective functionals, are very sensitive to a specific form of the functional and, in particular, to the value of the exponent. at is, the optimal controls obtained for such objective functionals are not robust, and a small variation of the exponent or even the weights leads to significant changes of the optimal controls. Furthermore, the above-mentioned approach does not guarantee the optimality of the obtained solutions. Instead, it ensures only necessary conditions of the optimality and a local extremum of the functional. At the same time, while the objective functionals that are independent of the controls are free of the mentioned drawbacks, their use involves dealing with considerably more difficult problems, as the control-independent functionals usually lead to bang-bang controls and can have singular arcs. ere are also other issues of mathematical nature that are related to a form of the objective functional. A comprehensive discussion of mathematical aspects of this problem can be found in recently published book by Schättler and Ledzewicz [23] (p. 48–50). Nevertheless, despite this complexity, there is a certain progress in application of the control-independent functionals in biomedical control problems. Results that affirmatively demonstrate that the use of such functionals is possible and potentially highly beneficial, include publications of Ledzewicz, Schättler and their coauthors on the optimal control in anticancer therapy (see, e.g., [23, 24] and bibliography therein), and of the authors on the optimal waste water biotreatment [25, 26–28] and on the optimal control of epidemics [29–34]. 3 Abstract and Applied Analysis e objective functionals that include squares of the controls naturally occur in problems originated in mechanics and engineering, where square of a control is usually proportional to the energy required for the control action. A typical motivation for the use of the quadratic objective in HIV treatment (and in biomedical applications in general) is minimizing the “costs” in terms of both financial costs and the adverse side effects. at is, this specific form of objective functional implies that the cost and side effects increase as squares of the drug concentrations in patient’s bloodstream. However, while there is little doubt that the severity of side effects grows faster than linear with the growth of drug concentration, and, hence, an exponential functional response may be a good representation for this, the actual value of the exponent is uncertain. Moreover, it depends on patient’s individual response and varies for different patient and different drugs. is uncertainty combined with the outcome sensitivity to the exponent implies that for clinical practice the value of results, obtained for this particular type of objective functional, is questionable. Furthermore, the quadratic “cost” is hardly able to adequately represent the financial cost of the therapy either. Taking into consideration the uncertainty of the form of the real-life “cost” function and the sensitivity of the results to this, we consider an optimal control problem for HIV treatment disregarding the “costs”. Instead, we assume that the “costs”, in terms of both the side effects and the financial cost of the therapy, is limited by the maximal acceptable levels of medication. An argument in support of such a formulation is that the actual objectives of the therapy are: (i) prevention of the development of AIDS by a patient, and (ii) control of the epidemic in a society. For an individual patient, these objectives imply reducing virus load to an acceptable safe level. Considering the severity of the possible consequences of the control failure, individual “costs” of the therapy can be disregarded. We also assume that the both, financial cost of the therapy and the side effects are known and accepted prior the start of the therapy. An initial advance in this direction was made in [35], where the authors considered an optimal control problem for a reasonably simple 3-dimensional antiretroviral therapy model. e model in this paper was based on the 2-dimensional Wodarz model of HIV dynamics [36], where an equation representing the concentration of a drug inhibiting the incidence rate was added. e drug intake was considered as the single control. e optimal control problem for the model with a control-independent objective functional was state and solved, and the principal properties of the optimal controls were established. 1.3. Interaction of Several Simultaneously Acting Controls. e majority of above-mentioned papers on the optimal antiHIV therapy consider models with a single control, which correspond to the therapies where a single drug (a single action) is used. Very few of existing publications consider a possibility of simultaneous use of two drugs [15, 20, 37–40], whereas works exploring interaction of three or more controls are rather exceptional [41, 42]. On the other hand, the current anti-HIV therapy guidelines necessitate the use of the highly active antiretroviral therapy (HAART) that implies administering a combination of several (typically, three or four) antiretrovirals. Moreover, over thirty anti-HIV drugs are currently available, which, by their action, can be classified into six groups that affect six different stages of viral life cycle. (NRTI and NNRTI are in different groups but affect the same stage of viral replication). Such a surplus of drugs’ names and actions implies a large number of their possible combinations. Currently, the drug combinations used in HAART are arranged with the aims of (i) minimizing the adverse side effects and (ii) preventing the development of drug resistance. However, the existing surplus of the drugs and their actions presents an opportunity for designing combinations, as well as schedules of their administering, aimed to optimize the HAART efficacy as well. Apparently, this implies the need to explore interaction and interference of the six (that is, equal to the number of the HIV replication cycle stages that the existing antiretroviral drugs can inhibit) simultaneously applied antiretroviral actions. 1.4. Nonlinearity of Functional Responses. e use of a numerical technique in papers [7–10, 12–22] means that in each case results are obtained for a model with a specific parametrization. However, due to intrinsic complexity of biological systems, the actual forms of functional responses and the parametrization are usually unknown in detail, and one hardly can expect that the information available on a biosystem is, or will be, sufficient. is implies that outcome of a study can depend on factors, which are neither studied, nor understood with a sufficient degree of certainty. is, in turn, leads to a question about the reliability and relevance of the results obtained for models with specific forms of functional responses and parametrization. In the framework of the traditional formulation that was used in these papers, this problem cannot be avoided. is problem was recognized at a very early stage of existence of mathematical biology. One way of dealing with this problem was suggested in 1930s by Kolmogorov [43], who suggested to “reverse” the problem: instead of studying properties of a biological system with specifically defined functional responses and parametrization, Kolmogorov suggested to formulate a model where the functional responses are given by unspecified functions and then establish the properties (such as positivity, monotonicity or concavity/convexity, etc.) that these functions should have in order to ensure the possession of certain properties by the considered system. is Kolmogorov’s concept was successfully used in the global analysis, particularly in combination with the direct Lyapunov method [44–48]. However, it is obvious that this concept hardly can be employed in a combination with numerical analysis. In biomedical applications of the optimal control theory this problem was recognized as well. To the best authors’ knowledge, Horst Behncke was the first who employed the Kolmogorov’s concept to an optimal control problem in mathematical biology [49]. He considered the classic 2-dimensional SIR Kermack–McKendrick epidemic model of an isolated epidemic assuming that the infection rate (which is the most important of all functional responses in epidemic models) is given by an unspecified function constrained by a few conditions. Behncke considered this model under action of one of three controls; namely, vaccination of the susceptible Abstract and Applied Analysis 4and Applied Analysis 4 that our optimal controls are only bang-bang and that no singular arcs are possible. Next, the maximal number of switchings for each control can be analytically estimated and then the specific number of the switching moments for each control can be found numerically as a problem of the finite-dimensional optimization. To deal with this optimal control problem, we employ a following analytical technique: (i) Using Pontryagin Maximum Principle, we found differential equations for the adjoint variables and for the switching functions that entirely determine types of the corresponding optimal controls. (ii) en, analyzing the system of differential equations for switching functions, we prove that the optimal controls are bang-bang. (iii) Our next task is to find accurate estimates for the number of zeros of the switching functions which correspond to switching points of controls. In order to do this, we apply the generalized Rolle’s theorem [53] to the equations of the system for switching functions. However, a direct application of the theorem to these equations is impossible, and we have to convert the matrix of this system into the upper triangular form. In order to do this, we have to find a suitable transformation, which is usually defined by a system of nonautonomous quadratic differential equations, and to prove that solutions of these equations are defined on the entire time interval. One of the results of this paper is that we managed to reduce a 3-dimensional system of nonautonomous quadratic differential equations to two 2-dimensional systems: the analysis of these two 2-dimensional systems is simpler than that of a single 3-dimensional system. (iv) Having the maximal possible number of switchings found, we can immediately find the principal type of the optimal controls. (v) Finally, to solve the problem completely, one can find the specific number of switchings and their precise positions, which can be done numerically. us, our analytical technique allows to reduce the optimal control problem to a considerably simpler problem of finite-dimensional optimization. Despite this possibility, the obtained theoretical results are illustrated by numerical calculations using BOCOP–2.0.5 so ware package, and the corresponding conclusions are made. A crucial for our objectives feature of this approach is that the estimations for the number of switchings and the values of the controls on the subintervals can be found analytically. is implies that there is a significant flexibility with respect to particular forms of the functional responses in the model and, hence, a possibility to use this approach in combination with the Kolmogorov’s “reversion of a problem” concept; we exploit this possibility in the paper. An analytical methods used in this paper can be employed to analyze optimal controls for other possible objectives, as individuals, isolation of the infectious individuals, and education that reduces a probability of infection: these three controls comprise all controls possible for such a simple model. Considering arguments as above, that is taking into consideration the disparity between the costs of a control policy and the costs and losses inflicted by a full-scale epidemic, Behncke employed objective functionals independent of the controls and considered minimization of the cumulative number of the infected over a given time interval. Assuming that in each case only one control is employed and that the controls are bounded, he found qualitative forms of these optimal controls. A similar approach was later applied to more complicated problems. us, the optimal controls for a SIR model of endemically persistent diseases (i.e., a SIR model with demographic processes) where the incidence rate was given by an unspecified function were considered in [33]. In this control model, four simultaneously acting controls were considered and their principal qualitative types established. It is noteworthy that there are more variants of the optimal controls for this problem compared with the single epidemic models considered by Behncke. Moreover, comparison of these results with the controls found for the same model with the standard bilinear transmission rate in [30, 32] shows that the nonlinearity of incidence rate does not affect the principal properties of the control. A more complicated 3-dimensional SEIR model with an unspecified nonlinear incidence rate and five simultaneously acting controls (which comprise all controls possible for this model) was considered in [34]. 1.5. Objectives and Mathematical Technique of this Paper. e objective of this paper is to study and establish properties of the optimal, highly active antiretroviral therapy, at the same time avoiding the above-mentioned drawbacks of the preceding studies. Having this goal in mind, we consider a model of HIV within-host dynamics based on the classical Nowak– May model, which describe interaction of the susceptible and infected target cells and free virus particles [50, 51]. e Nowak–May model postulates that the susceptible target cells can be infected by free virus particles and a er the instance of infection move into the infected population; in their turn, the infected cells produce free virus particles that can infect other susceptible cells. Using the Nowak–May model as basis, we assume that the rate of infection is given by a unspecified function of the concentrations of susceptible cells and free virus particles, constrained by a few conditions. Into the model, we also add a term responsible for the loss of free virus due to infection (see [52] for discussion of this term). In order to mirror the HAART, we introduce six bounded controls that can act simultaneously into this model. ese six controls comprise all controls that are possible for this model and correspond to the six stages of the viral replication cycle that the current antiretrovirals can affect. For this model, we consider a control problem of minimizing the infection level and the viral load at the end point of a finite time interval. e objective functional that we consider is independent of the controls. While for our type of objective functional, the optimal controls are likely to be of the bang-bang type, theoretically, they can contain singular arcs. us in this paper, we prove 5 Abstract and Applied Analysis bilinear and unspecified nonlinear incidence rates is studied in detail both analytically and numerically in [50, 51, 54–56, 44–47]. In particular, and what is of importance for further analysis, it is known that the nonnegative octant of the 3 -dimensional space {(푥, 푦, 푧) ∈ R3 : 푥, 푦, 푧 ≥ 0} is a positive invariant set of system (2.1) and that all solutions of the system initiated in this set are bounded. In this paper, we assume that the incidence rate 푓(푥, 푧) satisfies the following conditions. Assumption 1. Let the function 푓(푥, 푧) satisfy: (i) 푓(푥, 푧) is a continuous differentiable function for all 푥, 푧 ≥ 0; (ii) 푓(푥, 푧) > 0 for all 푥, 푧 > 0, and 푓(0, 푧) = 푓(푥, 0) = 0 for all 푥, 푧 ≥ 0; (iii) Partial derivatives 푓 (푥, 푧) and 푓 (푥, 푧) are positive for all 푥, 푧 > 0 and nonnegative for all 푥 ≥ 0, 푧 = 0 and 푥 = 0, 푧 ≥ 0. Please note that this set of assumptions is less restrictive than these that are usually used in the literature (see, e.g., [44–47, 56]). To make system (2.1) controllable, we assume that on a pre-defined time interval [0, 푇] the system is an object of six independent controls ̃푢(푡), ̃푣(푡), ̃푤(푡), ̃푝(푡), ̃푞(푡) and ̃ 푟(푡), which have the following biological meanings: (i) ̃푢(푡) is a fraction of immune cells among the new target cells entering the system; (ii) ̃푤(푡) is a per capita rate of “immunization” of the susceptible cells that makes the “immunized” cells entirely unsusceptible to the virus; (iii) ̃푣(푡) is a magnitude of indirect measures that inhibit the incidence rate; (iv) ̃푝(푡) is a magnitude of indirect measures aimed that inhibit the free virus particles production by the infected cells; (v) ̃푞(푡) is a per capita rate of immobilization (or removal) of free virus particles from the system; (vi) ̃ 푟(푡) is a per capita rate of killing (or of removal by other means) of the infected cells. We would like to stress that these six controls comprise all interventions that are feasible or potentially possible for model (2.1). We assume that the controls satisfy the following constraints: Here the maximal values of control functions ̃푢 max , ̃푣 max , ̃푤 max , ̃푝 max , ̃푞 max and ̃ 푟 max are pre-defined positive constants. Assuming that these six controls are applied to model (2.1) on the time interval [0, 푇], we obtain the following control model: (2.2) 0 ≤ ̃푢(푡) ≤ ̃푢 max < 1, 0 ≤ ̃푣(푡) ≤ ̃푣 max , 0 ≤ ̃푤(푡) ≤ ̃푤 max , 0 ≤ ̃푝(푡) ≤ ̃푝 max , 0 ≤ ̃푞(푡) ≤ ̃푞 max , 0 ≤ ̃ 푟(푡) ≤ ̃ 푟 max . well as for other models of virus dynamics. ese findings have the immediate practical relevance providing answers to long-existing questions about the interference of drugs with different actions combinations of which are used in HAART and about the conditions and objectives for which the current therapy guidelines are optimal. e paper is organized as follows: In Section 2, we formulate the model and state the optimal control problem. In Section 3, we apply Pontryagin Maximum Principle to find the adjoint system and the system of nonautonomous differential equations for the switching functions. In Sections 4, we prove that the optimal controls are bang-bang. In Section 5, we find the estimates for the number of zeros of the switching functions, and thereby, in Section 6, we describe the principal types of the optimal controls that are possible for this problem. In Section 7, we provide results of computations that illustrate our analytical results. Section 8 contains conclusion and discussion of the results. Appendix A describes details of the technique of the reduction of a 3-dimensional matrix of a system of linear nonautonomous homogeneous differential equations for the switching functions to the upper triangular form. Appendix B is devoted to the verification of one of our auxiliary assumptions. 2. Model and the Optimal Control Problem We consider the following model of the within-host HIV dynamics: Here 푥(푡), 푦(푡) and 푧(푡) are the appropriate concentrations of uninfected (susceptible) target cells, infected target cells and free virus particles at time . e model postulates that the susceptible cells (for HIV these are CD4+ T helper cells) enter into the system (into the patient blood from the thymus, where CD4+ cells mature) at a constant rate , that the susceptible cells are infected by the free virus particles at a rate 푓(푥(푡), 푧(푡)) and the infected cells immediately move to the infected population, that the infected cells produce free virus particles at the rate , and that average life-spans of the susceptible cells, infected cells and free virus particles are, respectively, −1, −1 and −1 (naturally, for HIV ≥ ). e third equation of the model also includes a removal term representing the loss of virus due to infection of cells. ( is term is usually omitted in mathematical models of HIV dynamics as very small compared to the virus death rate 훾푧(푡). However, for low viral loads typical for patients under antiviral therapy, the impact of the removal of free virus due to infection can be of importance). e model (2.1) is nearly classical. In the literature it usually appears with a bilinear incidence rate and without the term representing the loss of free virus due to infection and is referred to as the Nowak–May model [51, 54]. To the best of the authors’ knowledge, a version of the model with the unspecified general nonlinear incidence rate 푓(푥, 푧) firstly appeared in [46]. e dynamics of this model with both (2.1) {{{ {{{ { ⋅ 푥 (푡) = 휆 − 푓(푥(푡), 푧(푡)) − 휇푥(푡), ⋅ 푦 (푡) = 푓(푥(푡), 푧(푡)) − 휎푦(푡), ⋅ 푧 (푡) = 푘푦(푡) − 푓(푥(푡), 푧(푡)) − 훾푧(푡), 푥(0) = 푥0, 푦(0) = 푦0, 푧(0) = 푧0; 푥0, 푦0, 푧0 > 0. Abstract and Applied Analysis 6and Applied Analysis 6


Antiretroviral Drugs and Guidelines.
e World Health Organization and the United Nations Programme on HIV/ AIDS (UNAIDS) estimated that over 36.9 millions of people were living with HIV/AIDS in 2017; other 1.8 million were newly infected with HIV in the same year [1,2]. Of these 36.9 millions of the infected, 21.7 millions were living on antiretroviral therapy [1,2]. ese figures mean that the therapy is currently covers less than 60 percent of the infected. One (and probably the major) of the reasons for the insufficient coverage is the high cost of the antiretroviral therapy, which ranges from almost US $20, 000 per patient per year in the United States (where costs of the therapy are the highest) [3] to US $10, 000-15, 000 per a person per year in other countries [4] (p. 90). Such costs make the treatment inaccessible for the majority of patients in developing countries, where more than 95 percent of infections and HIV-associated deaths occurred.
Starting from 1987, when the first antiretroviral drugs, Nucleoside and Nucleotide Reverse Transcriptase Inhibitors, were approved for HIV treatment, over 30 of antiretroviral drugs were introduced in medical practice. A typical action of these drugs is inhibition of one of the several phases of HIV reproduction cycle. Depending on the stage that is inhibited, the drugs can be broadly classified into one of the following six groups [1, 2, 5]: (1) Entry or Fusion Inhibitors (Enfuvirtide, Maraviroc) Interfere with Binding, Fusion, and Entry of Virus to the Host Cell by Blocking One of Several Targets. (2) Nucleoside reverse transcriptase inhibitors (NRTI) and nucleotide reverse transciptase inhibitors (NtRTI) is the oldest and the largest class of antiretroviral drugs that includes drugs such as Zidovudine, Didanosine, Stavudine, Lamivudine, Abacavir, Tenofovir, and combinations of these. e NRTI and NtRTI inhibit reverse transcription by being incorporated into the newly synthesized viral DNA strand as faulty nucleotides. (3) Nonnucleoside reverse transcriptase inhibitors, or NNRTI, (Nevirapine, Delavirdine, Efavirenz, Etravirine, Rilpivirine) also inhibit reverse transcriptase, but in a different way, by binding to an allosteric site of the enzyme. (4) Integrase inhibitors (Raltegravir, Elvitegravir, Dolutegravir) inhibit the enzyme integrase, which is responsible for integration of the viral DNA into the DNA of the infected cell. (5) Protease inhibitors (PI) inhibit the activity of protease, an enzyme that cut nascent proteins into specific pieces for the final assembly of new functioning virions. ere are ten approved protease inhibitors, namely: Saquinavir, Indinavir, Ritonavir, Nelfinavir, Amprenavir, Lopinavir/Ritonavir, Atazanavir, Fosamprenavir, Tipranavir, and Darunavir. (6) Maturation inhibitors inhibit the last step in gag processing where the viral capsid polyprotein is cleaved, blocking the conversion of the polyprotein into the mature capsid protein. As a result, the virions released have a defective core and are mainly noninfectious.
e guidelines for use of the antiretroviral drugs are a source of controversy within the medical community. us, there is diversity of opinion regarding the stage of infection and the CD4+ T cells count (which are the primary target cells for HIV, and which depletion leads to the immune system collapse and development of AIDS) when the therapy should be started. It is generally accepted that the objectives of the treatment are to maximally suppress the viral load and to maintain CD4+ T cells at an acceptable level. Other factors that have to be also taken into consideration are: (i) a possibility of individual intolerance and adverse side-effects; (ii) a possibility of developing drug resistance for a patient, and (iii) a possibility of emergence of drug-resistant strains, which thereby can be passed onto others. Moreover, the cost of the therapy, which is high, should be also taken into consideration. To reduce the possibility of developing drug resistance and of the emergence of drug-resistant strains and to minimize the side effects, a combination of several antiretroviral drugs (typically, of three or four) is usually used. is approach is known as highly active antiretroviral therapy, or HAART. Also, it is generally accepted that, in order to reduce a possibility of emergence of drug resistant strains, antiretroviral therapy, once started, should never be stopped: as current US Guidelines states, "Patients initiating antiretroviral therapy should be willing and able to commit to lifelong treatment" [6]. e above-mentioned controversy in formulating the guidelines implies that some sort of optimization can be applied, and that the problem should be quantified and explored mathematically. Of course, at the present stage, it would be hardly possible to take into consideration details such as individual sensitivity of a specific patient to a specific drug. However, some common principles that can assist in developing better guidelines can be established. A natural approach for this sort of problems is to use the tools and methods of the classic optimal control theory. e problem of HIV control and the applicability of the optimal control theory to this problem attracted attention of mathematicians fairly early: seminal works by Butler et al. [7] and Kirschner et al. [8] initiated extensive literature in this direction. A mathematical techniques suggested in these publications proved to be successful and was subsequently applied to a variety of mathematical model of HIV (see, e.g., [9][10][11][12][13][14][15][16][17][18][19][20][21][22] and bibliography therein). ese results brought important insight into the problem. However, since this advance utilized essentially the same mathematical techniques, all these results also share some drawbacks, that we have to discuss before we proceed to our study.

Form of the Objective Functional.
Starting from [7,8], all the above-mentioned publication (and, in fact, the vast majority of the literature on the optimal control in mathematical biology and medicine) deal with optimal control problems of maximization of the cumulative (integral) concentration of healthy CD4+ T cells and simultaneous minimization of the cost of the treatment (represented by the weighted squares of the considered controls). Control problems with such functionals are mathematically convenient because they never deal with singular control phenomenon and can be solved numerically as a two-point boundary value problem. A setback of such objective functionals is that the optimal controls, obtained for such objective functionals, are very sensitive to a specific form of the functional and, in particular, to the value of the exponent. at is, the optimal controls obtained for such objective functionals are not robust, and a small variation of the exponent or even the weights leads to significant changes of the optimal controls. Furthermore, the above-mentioned approach does not guarantee the optimality of the obtained solutions. Instead, it ensures only necessary conditions of the optimality and a local extremum of the functional.
At the same time, while the objective functionals that are independent of the controls are free of the mentioned drawbacks, their use involves dealing with considerably more difficult problems, as the control-independent functionals usually lead to bang-bang controls and can have singular arcs.
ere are also other issues of mathematical nature that are related to a form of the objective functional. A comprehensive discussion of mathematical aspects of this problem can be found in recently published book by Schättler and Ledzewicz [23] (p. [48][49][50]. Nevertheless, despite this complexity, there is a certain progress in application of the control-independent functionals in biomedical control problems. Results that affirmatively demonstrate that the use of such functionals is possible and potentially highly beneficial, include publications of Ledzewicz, Schättler and their coauthors on the optimal control in anticancer therapy (see, e.g., [23,24] and bibliography therein), and of the authors on the optimal waste water biotreatment [25,[26][27][28] and on the optimal control of epidemics [29][30][31][32][33][34]. e objective functionals that include squares of the controls naturally occur in problems originated in mechanics and engineering, where square of a control is usually proportional to the energy required for the control action. A typical motivation for the use of the quadratic objective in HIV treatment (and in biomedical applications in general) is minimizing the "costs" in terms of both financial costs and the adverse side effects. at is, this specific form of objective functional implies that the cost and side effects increase as squares of the drug concentrations in patient's bloodstream. However, while there is little doubt that the severity of side effects grows faster than linear with the growth of drug concentration, and, hence, an exponential functional response may be a good representation for this, the actual value of the exponent is uncertain. Moreover, it depends on patient's individual response and varies for different patient and different drugs. is uncertainty combined with the outcome sensitivity to the exponent implies that for clinical practice the value of results, obtained for this particular type of objective functional, is questionable. Furthermore, the quadratic "cost" is hardly able to adequately represent the financial cost of the therapy either.
Taking into consideration the uncertainty of the form of the real-life "cost" function and the sensitivity of the results to this, we consider an optimal control problem for HIV treatment disregarding the "costs". Instead, we assume that the "costs", in terms of both the side effects and the financial cost of the therapy, is limited by the maximal acceptable levels of medication. An argument in support of such a formulation is that the actual objectives of the therapy are: (i) prevention of the development of AIDS by a patient, and (ii) control of the epidemic in a society. For an individual patient, these objectives imply reducing virus load to an acceptable safe level. Considering the severity of the possible consequences of the control failure, individual "costs" of the therapy can be disregarded. We also assume that the both, financial cost of the therapy and the side effects are known and accepted prior the start of the therapy.
An initial advance in this direction was made in [35], where the authors considered an optimal control problem for a reasonably simple 3-dimensional antiretroviral therapy model. e model in this paper was based on the 2-dimensional Wodarz model of HIV dynamics [36], where an equation representing the concentration of a drug inhibiting the incidence rate was added. e drug intake was considered as the single control. e optimal control problem for the model with a control-independent objective functional was state and solved, and the principal properties of the optimal controls were established.

Interaction of Several Simultaneously Acting Controls.
e majority of above-mentioned papers on the optimal anti-HIV therapy consider models with a single control, which correspond to the therapies where a single drug (a single action) is used. Very few of existing publications consider a possibility of simultaneous use of two drugs [15,20,[37][38][39][40], whereas works exploring interaction of three or more controls are rather exceptional [41,42]. On the other hand, the current anti-HIV therapy guidelines necessitate the use of the highly active antiretroviral therapy (HAART) that implies administering a combination of several (typically, three or four) antiretrovirals. Moreover, over thirty anti-HIV drugs are currently available, which, by their action, can be classified into six groups that affect six different stages of viral life cycle. (NRTI and NNRTI are in different groups but affect the same stage of viral replication). Such a surplus of drugs' names and actions implies a large number of their possible combinations. Currently, the drug combinations used in HAART are arranged with the aims of (i) minimizing the adverse side effects and (ii) preventing the development of drug resistance. However, the existing surplus of the drugs and their actions presents an opportunity for designing combinations, as well as schedules of their administering, aimed to optimize the HAART efficacy as well. Apparently, this implies the need to explore interaction and interference of the six (that is, equal to the number of the HIV replication cycle stages that the existing antiretroviral drugs can inhibit) simultaneously applied antiretroviral actions.

Nonlinearity of Functional Responses.
e use of a numerical technique in papers [7][8][9][10][12][13][14][15][16][17][18][19][20][21][22] means that in each case results are obtained for a model with a specific parametrization. However, due to intrinsic complexity of biological systems, the actual forms of functional responses and the parametrization are usually unknown in detail, and one hardly can expect that the information available on a biosystem is, or will be, sufficient. is implies that outcome of a study can depend on factors, which are neither studied, nor understood with a sufficient degree of certainty. is, in turn, leads to a question about the reliability and relevance of the results obtained for models with specific forms of functional responses and parametrization. In the framework of the traditional formulation that was used in these papers, this problem cannot be avoided.
is problem was recognized at a very early stage of existence of mathematical biology. One way of dealing with this problem was suggested in 1930s by Kolmogorov [43], who suggested to "reverse" the problem: instead of studying properties of a biological system with specifically defined functional responses and parametrization, Kolmogorov suggested to formulate a model where the functional responses are given by unspecified functions and then establish the properties (such as positivity, monotonicity or concavity/convexity, etc.) that these functions should have in order to ensure the possession of certain properties by the considered system. is Kolmogorov's concept was successfully used in the global analysis, particularly in combination with the direct Lyapunov method [44][45][46][47][48]. However, it is obvious that this concept hardly can be employed in a combination with numerical analysis.
In biomedical applications of the optimal control theory this problem was recognized as well. To the best authors' knowledge, Horst Behncke was the first who employed the Kolmogorov's concept to an optimal control problem in mathematical biology [49]. He considered the classic 2-dimensional SIR Kermack-McKendrick epidemic model of an isolated epidemic assuming that the infection rate (which is the most important of all functional responses in epidemic models) is given by an unspecified function constrained by a few conditions. Behncke considered this model under action of one of three controls; namely, vaccination of the susceptible that our optimal controls are only bang-bang and that no singular arcs are possible. Next, the maximal number of switchings for each control can be analytically estimated and then the specific number of the switching moments for each control can be found numerically as a problem of the finite-dimensional optimization.
To deal with this optimal control problem, we employ a following analytical technique: (i) Using Pontryagin Maximum Principle, we found differential equations for the adjoint variables and for the switching functions that entirely determine types of the corresponding optimal controls. (ii) en, analyzing the system of differential equations for switching functions, we prove that the optimal controls are bang-bang. (iii) Our next task is to find accurate estimates for the number of zeros of the switching functions which correspond to switching points of controls. In order to do this, we apply the generalized Rolle's theorem [53] to the equations of the system for switching functions. However, a direct application of the theorem to these equations is impossible, and we have to convert the matrix of this system into the upper triangular form. In order to do this, we have to find a suitable transformation, which is usually defined by a system of nonautonomous quadratic differential equations, and to prove that solutions of these equations are defined on the entire time interval. One of the results of this paper is that we managed to reduce a 3-dimensional system of nonautonomous quadratic differential equations to two 2-dimensional systems: the analysis of these two 2-dimensional systems is simpler than that of a single 3-dimensional system. (iv) Having the maximal possible number of switchings found, we can immediately find the principal type of the optimal controls. (v) Finally, to solve the problem completely, one can find the specific number of switchings and their precise positions, which can be done numerically. us, our analytical technique allows to reduce the optimal control problem to a considerably simpler problem of finite-dimensional optimization. Despite this possibility, the obtained theoretical results are illustrated by numerical calculations using BOCOP-2.0.5 soware package, and the corresponding conclusions are made.
A crucial for our objectives feature of this approach is that the estimations for the number of switchings and the values of the controls on the subintervals can be found analytically. is implies that there is a significant flexibility with respect to particular forms of the functional responses in the model and, hence, a possibility to use this approach in combination with the Kolmogorov's "reversion of a problem" concept; we exploit this possibility in the paper. An analytical methods used in this paper can be employed to analyze optimal controls for other possible objectives, as individuals, isolation of the infectious individuals, and education that reduces a probability of infection: these three controls comprise all controls possible for such a simple model. Considering arguments as above, that is taking into consideration the disparity between the costs of a control policy and the costs and losses inflicted by a full-scale epidemic, Behncke employed objective functionals independent of the controls and considered minimization of the cumulative number of the infected over a given time interval. Assuming that in each case only one control is employed and that the controls are bounded, he found qualitative forms of these optimal controls.
A similar approach was later applied to more complicated problems. us, the optimal controls for a SIR model of endemically persistent diseases (i.e., a SIR model with demographic processes) where the incidence rate was given by an unspecified function were considered in [33]. In this control model, four simultaneously acting controls were considered and their principal qualitative types established. It is noteworthy that there are more variants of the optimal controls for this problem compared with the single epidemic models considered by Behncke. Moreover, comparison of these results with the controls found for the same model with the standard bilinear transmission rate in [30,32] shows that the nonlinearity of incidence rate does not affect the principal properties of the control. A more complicated 3-dimensional SEIR model with an unspecified nonlinear incidence rate and five simultaneously acting controls (which comprise all controls possible for this model) was considered in [34].

Objectives and Mathematical Technique of this Paper.
e objective of this paper is to study and establish properties of the optimal, highly active antiretroviral therapy, at the same time avoiding the above-mentioned drawbacks of the preceding studies. Having this goal in mind, we consider a model of HIV within-host dynamics based on the classical Nowak-May model, which describe interaction of the susceptible and infected target cells and free virus particles [50,51]. e Nowak-May model postulates that the susceptible target cells can be infected by free virus particles and a er the instance of infection move into the infected population; in their turn, the infected cells produce free virus particles that can infect other susceptible cells. Using the Nowak-May model as basis, we assume that the rate of infection is given by a unspecified function of the concentrations of susceptible cells and free virus particles, constrained by a few conditions. Into the model, we also add a term responsible for the loss of free virus due to infection (see [52] for discussion of this term).
In order to mirror the HAART, we introduce six bounded controls that can act simultaneously into this model. ese six controls comprise all controls that are possible for this model and correspond to the six stages of the viral replication cycle that the current antiretrovirals can affect. For this model, we consider a control problem of minimizing the infection level and the viral load at the end point of a finite time interval. e objective functional that we consider is independent of the controls.
While for our type of objective functional, the optimal controls are likely to be of the bang-bang type, theoretically, they can contain singular arcs. us in this paper, we prove bilinear and unspecified nonlinear incidence rates is studied in detail both analytically and numerically in [50,51,[54][55][56][44][45][46][47]. In particular, and what is of importance for further analysis, it is known that the nonnegative octant of the 3 -dimensional space 푥, 푦, 푧 ∈ R 3 : 푥, 푦, 푧 ≥ 0 is a positive invariant set of system (2.1) and that all solutions of the system initiated in this set are bounded.
In this paper, we assume that the incidence rate 푓(푥, 푧) satisfies the following conditions. Please note that this set of assumptions is less restrictive than these that are usually used in the literature (see, e.g., [44][45][46][47]56]).
To make system (2.1) controllable, we assume that on a pre-defined time interval [0, 푇] the system is an object of six independent controls 푢(푡), 푣(푡), 푤(푡), 푝(푡), 푞(푡) and 푟(푡), which have the following biological meanings: (i) 푢(푡) is a fraction of immune cells among the new target cells entering the system; (ii) 푤(푡) is a per capita rate of "immunization" of the susceptible cells that makes the "immunized" cells entirely unsusceptible to the virus; (iii) 푣(푡) is a magnitude of indirect measures that inhibit the incidence rate; (iv) 푝(푡) is a magnitude of indirect measures aimed that inhibit the free virus particles production by the infected cells; (v) 푞(푡) is a per capita rate of immobilization (or removal) of free virus particles from the system; (vi) 푟(푡) is a per capita rate of killing (or of removal by other means) of the infected cells.
We would like to stress that these six controls comprise all interventions that are feasible or potentially possible for model (2.1). We assume that the controls satisfy the following constraints: Here the maximal values of control functions 푢 max , 푣 max , 푤 max , 푝 max , 푞 max and 푟 max are pre-defined positive constants. Assuming that these six controls are applied to model (2.1) on the time interval [0, 푇], we obtain the following control model: well as for other models of virus dynamics. ese findings have the immediate practical relevance providing answers to long-existing questions about the interference of drugs with different actions combinations of which are used in HAART and about the conditions and objectives for which the current therapy guidelines are optimal. e paper is organized as follows: In Section 2, we formulate the model and state the optimal control problem. In Section 3, we apply Pontryagin Maximum Principle to find the adjoint system and the system of nonautonomous differential equations for the switching functions. In Sections 4, we prove that the optimal controls are bang-bang. In Section 5, we find the estimates for the number of zeros of the switching functions, and thereby, in Section 6, we describe the principal types of the optimal controls that are possible for this problem. In Section 7, we provide results of computations that illustrate our analytical results. Section 8 contains conclusion and discussion of the results. Appendix A describes details of the technique of the reduction of a 3-dimensional matrix of a system of linear nonautonomous homogeneous differential equations for the switching functions to the upper triangular form. Appendix B is devoted to the verification of one of our auxiliary assumptions.

Model and the Optimal Control Problem
We consider the following model of the within-host HIV dynamics: Here 푥(푡), 푦(푡) and 푧(푡) are the appropriate concentrations of uninfected (susceptible) target cells, infected target cells and free virus particles at time . e model postulates that the susceptible cells (for HIV these are CD4+ T helper cells) enter into the system (into the patient blood from the thymus, where CD4+ cells mature) at a constant rate , that the susceptible cells are infected by the free virus particles at a rate 푓(푥(푡), 푧(푡)) and the infected cells immediately move to the infected population, that the infected cells produce free virus particles at the rate , and that average life-spans of the susceptible cells, infected cells and free virus particles are, respectively, −1 , −1 and −1 (naturally, for HIV ≥ ). e third equation of the model also includes a removal term representing the loss of virus due to infection of cells. ( is term is usually omitted in mathematical models of HIV dynamics as very small compared to the virus death rate 훾푧(푡). However, for low viral loads typical for patients under antiviral therapy, the impact of the removal of free virus due to infection can be of importance). e model (2.1) is nearly classical. In the literature it usually appears with a bilinear incidence rate and without the term representing the loss of free virus due to infection and is referred to as the Nowak-May model [51,54]. To the best of the authors' knowledge, a version of the model with the unspecified general nonlinear incidence rate 푓(푥, 푧) firstly appeared in [46]. e dynamics of this model with both
Proof of this lemma is fairly straightforward and hence we omit it. Proofs of similar statements are given, for example, in [23,[57][58][59].
On the set of admissible controls Ω(푇) for system (2.11), we consider the problem of minimizing a weighted sum of concentrations of the infected cells and free virus particles at the end of time interval [0, 푇]: Such objective functional usually corresponds to the task of complete elimination of the infection. In (2.15), and are positive weights, such that ̸ = . Please note that the objective functional 퐽 푢, 푣, 푤, 푝, 푞, 푟 is independent of the controls.

Lemma 5.
Depending on values of the weighting coefficients α and β, for switching functions 퐿 (푡) and 퐿 (푡), the following statements hold: (3.14) corresponding optimal solution 푥 * (푡), 푦 * (푡), 푧 * (푡) of system (2.11), there exists a vector-function 휓 * (푡) = 휓 * 1 (푡), 휓 * 2 (푡), 휓 * 3 (푡) such that: is a nontrivial solution of the adjoint system In this system, it is easy to note that for almost all 푡 ∈ [0, 푇], and, therefore, the following relationships hold: . Proof. Firstly, we consider switching function 퐿 (푡). Let us suppose the contrary, that is that there is an interval everywhere on interval Δ . Hence, by (3.12) and (3.15), on Δ the controls 푝 * (푡) and 푞 * (푡) simultaneously take either the values min and max , or the values max and min . Moreover, by the inequality from (4.14) and the last equation in (4.9), on the interval Δ the switching function 퐿 (푡) has at most one zero. (For analysis, it is sufficient that it is not equal to zero identically.) Hence, by (3.13), 푟 * ∈ 푟 min ; 푟 max on Δ . Substituting (4.14) and the above-found values of 푝 * (푡) and 푞 * (푡) into the first equality in (4.9), we conclude that either Proof. One of the inequalities for function 퐿 (푡) immediately follows from the arguments similar to those presented in the proof of Lemma 4. Indeed, switching function 퐿 (푡) satisfies the first equation of system (3.16), and, hence, is absolutely continuous and, therefore, is also a continuous function. en, for ̸ = , by eorem on the Stability of the Sign of a Continuous Function, there is such that, depending on the sign of − , either the first inequality in (4.2), or the first inequality in (4.3), holds.
To get the inequality for 퐿 (푡), we integrate the second equation of system (3.16) on the interval [0, 푇]. e integration yields en, taking into consideration the constraints imposed on control 푣 * (푡) by (2.5), Assumption 1 and the already established inequality for 퐿 (푡), we immediately see that the corresponding inequality for function 퐿 (푡) in either (4.2) or (4.3) also necessary hold on the interval 푡 , 푇 . From the absolute continuity of switching function 퐿 (푡) and the eorem on the Stability of the Sign of a Continuous Function it is immediately follows that there is value , < , such that, depending on the sign of − , the second inequality either in (4.2) or in (4.3) holds. is completes the proof. ☐ Corollary 6. Lemmas 4 and 5 combined with formulas (3.9), (3.10), (3.12), (3.13) and (3.15), immediately lead to the following relationships for the optimal controls 푢 * (푡), 푣 * (푡), 푤 * (푡), 푝 * (푡), 푞 * (푡) and 푟 * (푡): e following Lemma is of principal importance for further analysis. (4.6) As all possible cases lead to a contradiction, we have to conclude that the hypothesis was incorrect and that the switching function 퐿 (푡) cannot be identically equal to zero on any subinterval of the interval [0, 푇].
Now we have to consider whether switching function 퐿 (t) can be equal to zero on a subinterval of the interval [0, 푇].
Let us assume that there is an interval Δ ⊂ [0, 푇] on which en, almost everywhere on this interval holds as well. Substituting (4.24) and (4.25) into the second equation of system (3.16) and using the constraints on control 푣 * (푡) from (2.5) and the properties of function 푓(푥, 푧) from Assumption 1, we have to conclude that the contradictory equality 퐿 (푡) = 0 holds for all 푡 ∈ Δ . Hence, our assumption is incorrect and the switching function 퐿 (푡) cannot be identically equal to zero on any subinterval of the interval [0, 푇].
Next, we consider a possibility that the switching function 퐿 (푡) is identically equal to zero on any subinterval of the interval [0, 푇]. Applying the same arguments as those for the function 퐿 (푡), but using the third equation of system (3.16) instead of the second equation, we come to the conclusion that the switching function 퐿 (푡) cannot be identically equal to zero on any subinterval as well.
Finally, we consider the switching function 퐿 (푡). Let us suppose the contrary and assume that there is the interval Δ ⊂ [0, 푇] on which en, almost everywhere on this interval holds. Again, we substitute (4.26) and (4.27) into the fourth equation of system (3.16). Using the constraints on control 푝 * (푡) from (2.5) we find the contradictory equality 퐿 (푡) = 0 for 푡 ∈ Δ . Hence the above-made assumption was incorrect and the switching function 퐿 (푡) is not identically equal to zero on any subinterval of the interval [0, 푇]. is completes the proof. ☐ Lemma 7 combined with equalities (3.8)-(3.13) and (3.15) immediately leads to the following Corollary. In addition to this corollary, we note that the controls 푢 * (푡), 푣 * (푡), 푤 * (푡), 푝 * (푡), 푞 * (푡) and 푟 * (푡) have no singular arcs (see [23,24] for details).    where 0 푟 is a constant. Here, the equality 퐿 0 푟 = 0 is impossible, since, by (3.14), in the case we come to the contradictory equalities (4.13). erefore, 퐿 0 푟 ̸ = 0, and then we have to conclude that on the interval Δ the control 푟 * (푡) takes either the value min or value max . Substituting (4.17) and the abovefound values of control 푤 * (푡) into the first equality of (4.9), we immediately conclude that either or must hold. However, neither of these equalities is possible as contradicting to the control restrictions (2.7) and (2.8) of Assumption 2. erefore, this case is also impossible. Differentiating the first equality of (4.9) on the interval Δ or on its subinterval Δ ὔ , and then using the corresponding differential equations, we get equality where 푤 * ∈ 푤 min ; 푤 max , 푝 * + 푞 * ∈ 푝 min + 푞 max ; 푝 max + 푞 min and 푟 * ∈ 푟 min ; 푟 max . e first equality in (4.9) and (4.22) form a linear homogeneous system with respect to unknown 퐿 (푡) and 퐿 (푡). By (4.21), the solution of this system must be nontrivial, and, hence, the determinant of this system must be equal to zero. at is, the equality must hold. However, by the control restrictions (2.7)-(2.10) of Assumption 2, this equality is impossible. erefore, this case is also impossible.  .7), we write down the following system of differential equations for the switching function 퐿 (푡) and auxiliary functions 퐺 (푡) and 푄 (푡), which is more convenient for the subsequent analysis than system (3.16): Let us rewrite system (5.9) in the following matrix form: Now, using a special change of variables, we reduce the matrix of system (5.10) to the upper-triangular form. is

Estimating the Number of Zeros of the Switching Functions
In previous section, we proved that the optimal controls 푢 * (푡), 푣 * (푡), 푤 * (푡), 푝 * (푡), 푞 * (푡) and 푟 * (푡) are the bang-bang functions with values 푢 min ; 푢 max , 푣 min ; 푣 max , 푤 min ; 푤 max , 푝 min ; 푝 max , 푞 min ; 푞 max and 푟 min ; 푟 max , respectively. In this section, we proceed to estimating the number of switchings of these controls. Equalities (3.8)-(3.13) and (3.15) show that for this task it suffices to estimate the number of zeros of the switching functions 퐿 (푡), 퐿 (푡), 퐿 (푡) and 퐿 (푡) defined by system (3.16). From the analysis of the second and third equations of this system, we conclude that if we have an estimate for the number of zeros of function 퐿 (푡), then, applying the Generalized Rolle's eorem (see, e.g., [53]) to each of these equations, we can find estimates for the number of zeros of functions 퐿 (푡) and 퐿 (푡) as well. e use of this theorem with the fourth equation of system (3.16) gives us an estimate of the number of zeros of function 퐿 (푡). us, the analysis of switching function 퐿 (푡) is the basis for completing our task. e method that allows to estimate the number of zeros of the function 퐿 (푡) consists in the subsequent application of the above-mentioned Generalized Rolle's eorem to the equations of system (3.16).
Let us consider a possibility of application of the Generalized Rolle's eorem to the first equation of system (3.16). e direct application of the theorem to this equation is impossible, because, apart of the term containing the function 퐿 (푡), the equation has two terms independent of 퐿 (푡). . Finally, arguments similar to these used by the authors in [34] allow to conclude that the functions 푤 * (푡) − 푟 * (푡) and 푝 * (푡) + 푞 * (푡) − 푟 * (푡) are also differentiable almost everywhere on the interval [0, 푇]. en, using the second and third equations of system (3.16), we formulate the required differential equation: , we obtain the following systems of equations: Let us apply to systems (5.18) and (5.19) with nonnegative initial conditions the necessary and sufficient condition for solutions 푔 20 (푡), 푔 21 (푡) and 푔 31 (푡), 푔 32 (푡) to be nonnegative for all 푡 ∈ 0, 푡 2 1 and 푡 ∈ 0, 푡 3 1 , respectively [65]. It consists in satisfying the so-called quasi-positivity condition for system (5.18), and for system (5.19).
As a result, we obtain the inequalities: transformation is fairly standard, and is described in detail in Appendix A. Comparing systems (5.10) and (A.1) and applying arguments from Appendix A, we find the desired system in the form Here 푘 푣 (푡) = 푑 푣 (푡) − 0.25푒 2 푣 (푡) is a piecewise-constant function as well, and, as it is not difficult to see, the inequality 푘 (푡) < 0 We have to note that these inequalities are obtained under constraints which follow from the relationships (see Lemma in Section 14, Chapter 4 in [66]). at is, the estimates imply the required continuation of the nonnegative solutions 푔 20 (푡) and 푔 21 (푡) of system (5.18) and the nonnegative solutions 푔 31 (푡) and 푔 32 (푡) of system (5.19) to the entire interval [0, 푇]. Taking into consideration that functions ℎ 20 (푡), ℎ 21 (푡), ℎ 31 (푡) and ℎ 32 (푡) satisfy systems (5.12) and (5.13) and are also defined on the entire interval [0, 푇], we conclude that system (5.11) is defined on this interval as well.
Applying this theorem to the second equation of system (5.11), we have to conclude that the function Denote the zeros of function 퐿 (푡) by * 1 and * 2 ; 휏 * 1 , 휏 * 2 ∈ [0, 푇). Depending on the sign of the difference − , the function 퐿 (푡) is described by one of the following two relationships, namely: We consider the following assumption, the validity of which will be demonstrated in subsequent arguments. Finally, we have to verify the continuation of the nonnegative solutions 푔 20 (푡), 푔 21 (푡), 푔 31 (푡) and 푔 32 (푡) of systems (5.18) and (5.19) to the entire interval [0, 푇] (see eorem 1.6 in [65]). To do this, we define Lyapunov functions: and estimate the derivatives 푑푉 2 /푑푡 and 푑푉 3 /푑푡 of these functions along the trajectories of the systems of differential equations (5.18) and (5.19), respectively. ( e derivatives are the scalar products of the gradients of these functions and the right-hand sides of the corresponding systems.) e Lyapunov functions satisfy Here, 2 and 3 are constants defined as following: (5.28) To estimate the number of zeros of the switching function 퐿 (푡), we take into consideration the boundary condition 퐿 (푇) = 훽 > 0 (see (3.16)). Applying the arguments that are analogous to those used to estimate the number of zeros of function 퐿 (푡), we can prove that the function 퐿 (푡) has at most three distinct zeros on the interval [0, 푇]. We omit the arguments, which are the same as the above-used. Taking into consideration the first inequality from (4.1) in Lemma 4, we have to conclude that the switching function 퐿 (푡) satisfies Here 휂 * 1 , 휂 * 2 , 휂 * 3 ∈ [0, 푇) are zeros of the switching function 퐿 (푡). Finally, we estimate the number of zeros of switching function 퐿 (푡). Taking into consideration the boundary condition 퐿 (푇) = 훼 > 0 (see (3.16)), we show that function 퐿 (푡) has at most four distinct zeros on interval [0, 푇). Let us prove this by contradiction: assume that the function 퐿 (푡) has at least five distinct zeros on the interval [0, 푇). Applying the Generalized Rolle's eorem, we conclude that in this case the function 퐿 (푡) has at least four distinct zeros on the interval (0, 푇). is conclusion contradicts to the above-established estimate, and, therefore, our assumption is incorrect.
Denoting the zeros of function 퐿 (푡) by * 1 , * 2 , * 3 and * 4 ; 휒 * 1 , 휒 * 2 , 휒 * 3 , 휒 * 4 ∈ (0, 푇), and taking into consideration the second inequality from (4.1) in Lemma 4, we get the following relationship describing behavior of the switching function 퐿 (푡): 6. Types of the Optimal Controls e above-found qualitative profiles of the switching functions (5.33)-(5.38) combined with the formulas of the optimal controls (3.9), (3.10) and (3.12), (3.13) enable us to immediately find the possible types of the optimal controls 푣 * (푡), 푤 * (푡), 푞 * (푡) and 푟 * (푡). Specifically, depending on the sign of − , the optimal controls 푣 * (푡) and 푤 * (푡) are of one of the following types: Please note that these formulas completely agree with the corresponding inequalities from (4.2) and (4.3) of Lemma 5. Now we proceed to estimating the number of zeros of the switching function 퐿 (푡). We will show that, taking into consideration that 퐿 (푇) = 0 (see (3.16)), function 퐿 (푡] has at most three distinct zeros on the interval [0, 푇). To prove this preposition, we, again, suppose the opposite: that is, we suppose that function 퐿 (푡) has at least four distinct zeroes on interval [0, 푇]. Applying the Generalized Rolle's eorem, we have to conclude that in this case function 퐿 (푡) has at least three distinct zeroes on the interval (0, 푇), which contradicts to the above-found estimates for the number of zeros of switching function 퐿 (푡). erefore, our assumption was incorrect, and the required estimate for the number of zeros of the switching function 퐿 (푡) is found.
Since 퐿 (푇) = 0 (see (3.16)), one of these three zeros is the end point of the interval [0, 푇]. We denote the other two zeros by * 1 and * 2 ; 휃 * 1 , 휃 * 2 ∈ (0, 푇). en, taking into consideration the corresponding inequalities from Lemma 5, we conclude that, depending on the sign of − , the switching function 퐿 (푡) is described by one of the following two relationships, namely:

Numerical Study
In order to illustrate our analytical results, we consider some numerical examples for problem (2.11), (2.15) using BOCOP-2.0.5 so ware package. BOCOP-2.0.5 is an optimal control interface, which is specifically designed for solving optimal control problems with general path and boundary constraints and with free or fixed final time and implemented in MATLAB [67]. Using discretization of time, problems under consideration are approximated by finite-dimensional optimization problems, which are thereby solved by the well-known opensource so ware package IPOPT designed for large-scale nonlinear optimization.
In our calculations, we use a time grid with 5000 nodes over time interval [0, 푇] that is equivalent to the time step of Δ푡 = 2푇 ⋅ 10 −4 . Since the problem is solved by a direct method using an iterative approach, at each step we impose the acceptable convergence tolerance 휀 푟푒푙 = 10 −14 . In the computations we use the sixth-order Lobatto III C discretization rule (see [67] for details). Parameters and initial conditions of system (2.11), the weights in (2.15), control constraints in (2.5) are as following: ese values are adopted from [68]. We note that > .

Conclusions
In this paper, we used an optimal control theory to model the highly active antiretroviral therapy for the purpose of finding both an optimal scheduler and an optimal drug combination to be used in HAART. To address this issue, we considered a compartmental model of within-host HIV dynamics based on the Nowak-May model. Compared with the classical Nowak-May model, our model assumes that the incidence rate is given by an unspecified nonlinear function 푓(푥, 푧) constrained by (see (2.15)). IPOPT used 50 iterations in order to solve the optimal control problem (2.11), (2.15). Figures 7 and 8 present results for the nonlinear incidence rate 푓(푥, 푧) of (7.3). Figure 7 show the optimal solutions 푥 * (푡), 푦 * (푡) and 푧 * (푡), while Figure 8 display the corresponding optimal controls 푢 * (푡), 푣 * (푡), 푤 * (푡), 푝 * (푡), 푞 * (푡) and 푟 * (푡). 퐽 * = 1.425958 ⋅ 10 −3 is the minimum value of the functional 퐽 푢, 푣, 푤, 푝, 푞, 푟 of (2.15). IPOPT used 55 iterations in order to solve the optimal control problem (2.11) and (2.15). e numerical calculations lead to the following conclusions: (i) A change in the type of the incidence rate 푓(푥, 푧) does not lead to qualitative changes in the behavior of the optimal solutions 푥 * (푡), 푦 * (푡) and 푧 * (푡).  problem of minimizing the viral load and the infected cell concentration at the end of a predetermined time interval is stated and solved with the objective functional independent on the controls. Using Pontryagin Maximum Principle, we investigated optimal dynamic analytically. We proved that the optimal controls are bang-bang with no singular arcs and found all possible variants of the optimal controls. e realization of a specific variant depends on specific initial conditions 푥(0), 푦(0), 푧(0) and the maximal values of derivatives 푓 ὔ (푥, 푧) and 푓 ὔ (푥, 푧) reached in the feasible region of the model. When the maximal number of switchings of the optimal controls is known, the two-point boundary value problem for Pontryagin Maximum Principle can be immediately reduced to a considerably simpler problem of finite-dimensional optimization. For specific parameter values and initial a few biologically feasible conditions and takes account of the loss of free virus due to the target cell infections. e constraints imposed on incidence rate 푓(푥, 푧) summarized in Assumption 1 are biologically motivated and hold for any realistic incidence rate. ese constraints are considerably less restrictive than the conditions for the uniqueness of a positive equilibrium state and the global asymptotic stability of the system [45][46][47][48]. Apart from biologically motivated Assumption 1 and for convenience of the analysis, we also introduce constraints in Assumption 9. Further, our numerical experiments show that constraints in Assumption 9 are sufficient but not necessary conditions.
To mirror HAART, we introduce six simultaneously acting bounded controls into this model. ese six controls comprise all controls possible for this model and represent the actions that the existing anti-HIV drugs exhibit. e optimal control  (i) e first group is composed of controls 푢(푡) ("immunization" of CD4+ T cells before they enter patient's bloodstream) and 푤(푡) ("immunization" of the existing CD4+ T cells) and corresponds to action of the fusion inhibitors. (Drugs that "immunize" CD4+ T cells before they enter patient's blood stream do not exist at the moment, but such drugs can appear in the nearest future). (ii) e second group combines controls 푝(푡) (inhibition of the virus production by the infected cells) and 푞(푡) (immobilization or removal of the free virus particles). e first of these controls correspond to actions of the protease and maturation inhibitors. Drugs represented by control 푞(푡) are not currently available, but they can appear in the nearest future. Action of such drugs is similar to that of antibodies. (iii) e third and the fourth groups contain a single control each; namely, 푟(푡) (killing or immobilizing of the infected cells, or removing them by other means) and It is remarkable that in the absence of control ∼ 푣 (푡) the drugs of the first, second and third groups are completely independent and do not interfere. is implies that in the absence of control ∼ 푣 (푡), the controls of other three groups can be considered separately (which simplifies the analysis), whereas in the presence of this control all four groups should be considered as a complex. What is more important for real-life conditions, the latter problem can be solved numerically. However, in this paper, we did not use this approach; we used, instead, an alternative numerical solver package BOCOP-2.0.5 to verify and illustrate our analytical findings.
It follows from our analysis that for the considered model with a nonlinear incidence rate satisfying few biologically feasible assumptions, the principal qualitative properties of the optimal controls, such as the maximum number of switchings of the bang-bang controls and the order of these switchings, do not depend on a particular form of the nonlinearity. Additionally, our analysis indicates that possible nonlinearities of other functional responses in the model would not affect the principal properties of the optimal controls either, provided that the signs of the derivatives of these responses are preserved. We have to stress, however, that the conclusion that the principal properties of the optimal controls do not depend on a form of the nonlinearity of functional responses of a model holds only for the functionals which are independent of the controls and can be different for the explicitly depending on controls objective functionals.
While our numerical outcomes are presented for treatment duration of 60 days, we proved in this paper that our results do not depend on the length of the treatment and that the optimal schedule and drug combination would be the same for different time intervals.
Another practically relevant conclusion that follows from our analysis is that the six considered controls can be segregated into four groups such that the controls of each group should be considered as a single control and started, carried out and ended simultaneously. Here, ℎ 21 (푡), ℎ 31 (푡) and ℎ 32 (푡) are functions should be chosen later. It is easy to see that such a transformation does not change function 푟 1 (푡); that is, ∼ 푟 1 (푡) = 푟 1 (푡) holds. e other two functions satisfy Moreover, the reverse transformation is also easily determined as e corresponding inverse matrix in equality We apply transformation (A.5) to matrix 퐹(푡). To do this, using variables 푟 (푡), we rewrite system (A.1) as Here we choose matrix 퐻(푡), or that is the same, the functions ℎ 21 (푡), ℎ 31 (푡) and ℎ 32 (푡), to satisfy the equality We re-write equality (A.11) element-wise: e matrix 퐷(푡) that satisfies these equalities has the required upper-triangular form (A.2). On its main diagonal the matrix has values (A.7) practice, this fact also means that a drug or drugs of the fourth control group should be included in HAART. For a continuing uninterrupted HAART, these results also indicate that life-long use of the same drug combination may be not optimal, and that patient-specific alterations of the combinations can be beneficial. However, in this case other considerations, and in particular a possibility of developing of drug-resistant strains, should be taken into consideration as well.
Finally, we note that the qualitative descriptions of the optimal controls (6.1)-(6.7) assuming that constrains (2.7)-(2.10) of Assumption 2 are held. However, while these constrains significantly simplify the analysis, they are not necessary: by arguments that were used in [29,34], it is possible to establish the validity of relationships (6.1)-(6.7) for the cases where conditions (2.7)-(2.10) of Assumption 2 are violated.
Our goal is to reduce the matrix 퐹(푡) to the upper-triangular form: In our earlier papers [28,34], we conducted such a transformation via a series of successive changes of variables. However, such a transformation can be achieved by a simultaneous nondegenerate transformation of a vector 푟(푡) or, in a matrix form,

Abstract and Applied Analysis 20
To formulate the second 2-dimensional subsystem, we introduce variable We would like to emphasize that the 3-dimensional quadratic system (A.18) is equivalent to two 2-dimensional quadratic subsystems (A.21) and (A.26). As analysis of 2 -dimensional systems is usually simpler, further we will deal with the later. In particularly, we have to analyze these quadratic subsystems to find the restrictions on functions 푓 (푡) which guarantee that there exist solutions of subsystems (A.21) and (A.26) defined on the entire interval [0, 푇].

B. Verifying Assumption 9
To verify that Assumption 9 holds for the optimal controls that we found in this paper and the corresponding optimal solutions, we consider system (2.11) for the case of bilinear incidence rate 푓(푥, 푧) = , where is a positive constant: and, due to equalities (A.12) and (A.14), functions ℎ 21 (푡) and ℎ 32 (푡) satisfy the following differential equations: To find a similar equation for function ℎ 31 (푡), we substitute (A.14) into equality (A.13). is yields the equation where ℎ 0 21 , ℎ 0 31 and ℎ 0 32 are constants, can be added to system (A.18). As a result, we have the Cauchy problem, and, as a consequence of the Existence and Uniqueness of Solution eorem [64], solutions ℎ 21 (푡), ℎ 31 (푡), ℎ 32 (푡) to this problem can be defined only locally, in a neighborhood of the value 푡 = 0. However, for this problem it is important to have solutions of the Cauchy problem (A.18), (A.20) defined on the entire interval [0, 푇] (see [28,34]). at is we have to find the restrictions on functions 푓 (푡) such that under these restrictions solutions of problem (A.18) and (A.20) were defined on the entire interval [0, 푇]. en system (A.19) would be also defined on the entire interval and, hence, it would be possible to use the Generalized Rolle's eorem [53] to estimate the number of zeros of switching function 푟 1 (푡) = ∼ 푟 1 (푡). e analysis of the 3-dimensional quadratic system (A.18) can be significantly simplified by splitting this system into two simpler 2-dimensional quadratic subsystems. Specifically, it is easy to see that equations for ℎ 31 (푡) and ℎ 32 (푡) contain no variable ℎ 21 (푡) and, hence, can be considered separately. ese two equations form the first of these quadratic subsystems: e value max is also defined in (2.12). e third inequality in (B.2) is satisfied if the smallest possible values of its le -hand side is greater than or equal to zero. Hence, we have the following inequalities: Combining inequalities (B.3) and (B.10), we obtain the following system: If these inequalities hold, then inequalities (B.2) are satisfied as well. We verified the validity of inequalities (B.11) by direct calculations for the following parameters and initial conditions of system (2.11) and specific values min , min , max , min , max , min , min , max , min and max from (2.5): Other specific values of the parameters and initial conditions of system (2.11) can be found in [69,70].
Data Availability e so ware used to support the findings of this study is free open-access so ware that can be freely downloaded from https://www.bocop.org/. e values of parameters and initial conditions used in computations to support the findings of this study are included within the article. 휅푣 min 푧 min + 푥 min − 0.5 푤 max − 푟 min + 푞 max − 푟 min ≥ 0, 1 − 2푤 max 푤 min − 푟 max 푧 min + 1 − 2푞 max 푝 min + 푞 min − 푟 max 푥 min ≥ 0.

For this case Assumption 9 is
Our goal is to find the restrictions on the parameters and the initial conditions of system (2.11), as well as on values min , max , min , max , min , max , min , max , min , max , min and max defined in (2.5), which ensure that inequalities (B.2) are held. Taking into consideration that all the controls are bangbang and taking either the minimal, or the maximal possible value, we immediately see that the first two inequalities in (B.2) hold if the smallest values of their le -hand sides are greater than or equal to the largest values of the righthand sides. at is, we have the inequalities: Here min and min are the positive lower bounds for solutions 푥(푡) and 푧(푡), respectively. In inequalities (2.12) of Lemma 3 for the general incidence rate 푓(푥, 푧), these values are zero. However, for the specific function 푓(푥, 푧) = the lower bounds min and min can be refined. Specifically, integrating the first equation in (B.1), we get the inequality where To find the lower bound min , we, firstly, have to find the positive lower bound min for 푦(푡) by integrating the second equation of system (B.1). is yields where 푦 min = 푦 0 e −푟 max 푇 . en, using this inequality, we integrate the third equation of system (B.1) and get the inequality 휅푣 min 푤 min − 푟 max 푧 min + 푝 min + 푞 min − 푟 max 푥 min ≥ 푤 max 푞 max , 푤 min 푞 min 휅푣 min 푧 min + 푥 min + 푟 min + 휅푣 min 푤 min 푤 min − 푟 max 푧 min +푞 min 푝 min + 푞 min − 푟 max 푥 min ≥ 푤 max 푞 max 푤 max + 푞 max .