# Optimal Controls of the Highly Active Antiretroviral Therapy

Authors: Ellina V. Grigorieva; Evgenii N. Khailov; Andrei Korobeinikov
Journal: Abstract and Applied Analysis
 (2020)

---

## Abstract

In this paper, we study generic properties of the optimal, in a certain sense, highly active antiretroviral therapy (or HAART). To address this problem, we consider a control model based on the3-dimensional Nowak–May within-host HIV dynamics model. Taking into consideration that precise forms of functional responses are usually unknown, we introduce into this model a nonlinear incidence rate of a rather general form given by an unspecified function of the susceptible cells and free virus particles. We also add a term responsible to the loss of free virions due to infection of the target cells. To mirror the idea of highly active anti-HIV therapy, in this model we assume six controls that can act simultaneously. These six controls affecting different stage of virus life cycle comprise all controls possible for this model and account for all feasible actions of the existing anti-HIV drugs. With this control model, we consider an optimal control problem of minimizing the infection level at the end of a given time interval. Using an analytical mathematical technique, we prove that the optimal controls are bang-bang, find accurate estimates for the maximal possible number of switchings of these controls and establish qualitative types of the optimal controls as well as mutual relationships between them. Having the estimate for the number of switchings found, we can reduce the two-point boundary value problem for Pontryagin Maximum Principle to a considerably simpler problem of the finite-dimensional optimization, which can be solved numerically. Despite this possibility, the obtained theoretical results are illustrated by numerical calculations using BOCOP–2.0.5 software package, and the corresponding conclusions are made.


---

## Body

## 1. Introduction

### 1.1. Antiretroviral Drugs and Guidelines

The World Health Organization and the United Nations Programme on HIV/AIDS (UNAIDS) estimated that over36.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]. These 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, over30 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. The 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. There 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.The guidelines for use of the antiretroviral drugs are a source of controversy within the medical community. Thus, 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. This 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].The 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. The 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). These 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. That 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. There 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].The 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. That 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. This 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. The 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. The drug intake was considered as the single control. The 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

The 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–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

The 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 bio-system is, or will be, sufficient. This implies that outcome of a study can depend on factors, which are neither studied, nor understood with a sufficient degree of certainty. This, 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.This 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. This 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 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. Thus, 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

The 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]. The Nowak–May model postulates that the susceptible target cells can be infected by free virus particles and after 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. These 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. The 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. Thus in this paper, we prove 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)
Then, 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. Thus, 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 software 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. This 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 well as for other models of virus dynamics. These 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.The paper is organized as follows: In Section2, 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.


## 1.1. Antiretroviral Drugs and Guidelines

The World Health Organization and the United Nations Programme on HIV/AIDS (UNAIDS) estimated that over36.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]. These 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, over30 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. The 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. There 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.The guidelines for use of the antiretroviral drugs are a source of controversy within the medical community. Thus, 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. This 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].The 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. The 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). These 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. That 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. There 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].The 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. That 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. This 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. The 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. The drug intake was considered as the single control. The 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

The 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–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

The 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 bio-system is, or will be, sufficient. This implies that outcome of a study can depend on factors, which are neither studied, nor understood with a sufficient degree of certainty. This, 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.This 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. This 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 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. Thus, 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

The 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]. The Nowak–May model postulates that the susceptible target cells can be infected by free virus particles and after 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. These 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. The 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. Thus in this paper, we prove 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)
Then, 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. Thus, 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 software 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. This 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 well as for other models of virus dynamics. These 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.The paper is organized as follows: In Section2, 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:(2.1)x˙t=λ−fxt,zt−μxt,y˙t=fxt,zt−σyt,z˙t=kyt−fxt,zt−γzt,x0=x0,y0=y0,z0=z0;x0,y0,z0>0.Herext, yt and zt are the appropriate concentrations of uninfected (susceptible) target cells, infected target cells and free virus particles at time t. The 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 fxt,zt and the infected cells immediately move to the infected population, that the infected cells produce free virus particles at the rate ky, and that average life-spans of the susceptible cells, infected cells and free virus particles are, respectively, μ−1, σ−1 and γ−1 (naturally, for HIV σ≥μ). The third equation of the model also includes a removal term representing the loss of virus due to infection of cells. (This term is usually omitted in mathematical models of HIV dynamics as very small compared to the virus death rate γzt. 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).The 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 fx,z firstly appeared in [46]. The dynamics of this model with both 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 x,y,z∈R3:x,y,z≥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 ratefx,z satisfies the following conditions.Assumption 1.
Let the functionfx,z satisfy:(i)
fx,z is a continuous differentiable function for all x,z≥0;(ii)
fx,z>0 for all x,z>0, and f0,z=fx,0=0 for all x,z≥0;(iii)
Partial derivativesfx′x,z and fz′x,z are positive for all x,z>0 and nonnegative for all x≥0,z=0 and x=0,z≥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,T the system is an object of six independent controls u˜t, v˜t, w˜t, p˜t, q˜t and r˜t, which have the following biological meanings:(i)
u˜t is a fraction of immune cells among the new target cells entering the system;(ii)
w˜t is a per capita rate of “immunization” of the susceptible cells that makes the “immunized” cells entirely unsusceptible to the virus;(iii)
v˜t is a magnitude of indirect measures that inhibit the incidence rate;(iv)
p˜t is a magnitude of indirect measures aimed that inhibit the free virus particles production by the infected cells;(v)
q˜t is a per capita rate of immobilization (or removal) of free virus particles from the system;(vi)
r˜t 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:(2.2)0≤u˜t≤u˜max<1,0≤v˜t≤v˜max,0≤w˜t≤w˜max,0≤p˜t≤p˜max,0≤q˜t≤q˜max,0≤r˜t≤r˜max.Here the maximal values of control functionsu˜max, v˜max, w˜max, p˜max, q˜max and r˜max are pre-defined positive constants.Assuming that these six controls are applied to model (2.1) on the time interval 0,T, we obtain the following control model:(2.3)x˙t=λ1−u˜t−1+v˜t−1fxt,zt−μ+w˜txt,y˙t=1+v˜t−1fxt,zt−σ+r˜tyt,z˙t=k1+p˜t−1yt−1+v˜t−1fxt,zt−γ+q˜tzt,x0=x0,y0=y0,z0=z0;x0,y0,z0>0.For convenience of notation, we introduce modified controls:(2.4)ut=λ1−u˜t,vt=1+v˜t−1,wt=μ+w˜t,pt=1+p˜t−1,qt=γ+q˜t,rt=σ+r˜t.These modified controls satisfy the following constraints:(2.5)umin≤ut≤umax,vmin≤vt≤vmax,wmin≤wt≤wmax,pmin≤pt≤pmax,qmin≤qt≤qmax,rmin≤rt≤rmax,where(2.6)umin=λ1−u˜max,umax=λ,vmin=1+v˜max−1,vmax=1,wmax=μ+w˜max,wmin=μ,pmin=k1+p˜max−1,pmax=k,qmax=γ+q˜maxqmin=γ,rmax=σ+r˜max,rmin=σ.In addition to Assumption1, we make the following assumption.Assumption 2.
Maximal and minimal levelswmin, wmax, pmin, pmax, qmin, qmax, rmin, rmax, of the controls satisfy the conditions:(2.7)wmin≠qmin,qmax,rmin,rmax,
(2.8)
w
m
a
x
≠
q
min
,
q
m
a
x
,
r
min
,
r
m
a
x
,
(2.9)
r
min
≠
p
min
+
q
m
a
x
,
p
m
a
x
+
q
min
,
(2.10)
r
m
a
x
≠
p
min
+
q
m
a
x
,
p
m
a
x
+
q
min
.As a set of admissible controls,ΩT, we consider the set of all Lebesque measurable functions ut, vt, wt, pt, qt and rt satisfying inequalities (2.5) for almost all t∈0,T.Using modified controls (2.4), we rewrite control system (2.3) as following:(2.11)x˙t=ut−vtfxt,zt−wtxt,y˙t=vtfxt,zt−rtyt,z˙t=ptyt−vtfxt,zt−qtzt,x0=x0,y0=y0,z0=z0;x0,y0,z0>0.The following lemma ensures the boundedness, positiveness, and continuity of solutions for system (2.11).Lemma 3.
For any admissible controlsut, vt, wt,pt, qt, and rt, the corresponding absolutely continuous solution xt, yt, zt to system (2.11) are defined on the entire interval 0,T, and for its components xt, yt, zt the inequalities:(2.12)0<xt<xmax,0<yt<ymax,0<zt<zmax,
where(2.13)xmax,ymax=x0+y0+λminμ;σ,
(2.14)
z
m
a
x
=
y
0
+
z
0
+
k
y
m
a
x
min
σ
;
γ
,
hold for allt∈0,T.Proof of this lemma is fairly straightforward and hence we omit it. Proofs of similar statements are given, for example, in [23, 57–59].On the set of admissible controlsΩT 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,T:(2.15)Ju,v,w,p,q,r=αyT+βzT→minu·,v·,w·,p·,q·,r·.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 Ju,v,w,p,q,r is independent of the controls.For optimal control problem (2.11), (2.15), inequalities (2.12) in Lemma 3 and Theorem 4 of Chapter 4 in [60] ensure the existence of the optimal controls u∗t, v∗t, w∗t, p∗t, q∗t and r∗t and the corresponding optimal solution x∗t,y∗t,z∗t of system (2.11). In the following sections we establish the principal properties of these controls.


## 3. Pontryagin Maximum Principle

To analyze the optimal control problem (2.11), (2.15), we use Pontryagin Maximum Principle [61]. The Hamiltonian of this problem is(3.1)Hx,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=u−vfx,z−wxψ1+vfx,z−ryψ2+py−vfx,z−qzψ3,whereψ1, ψ2 and ψ3 are adjoint variables. The Hamiltonian satisfies:(3.2)∂H∂xx,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=vfx′x,zψ2−ψ1−ψ3−wψ1,∂H∂yx,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=−rψ2+pψ3,∂H∂zx,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=vfz′x,zψ2−ψ1−ψ3−qψ3,∂H∂ux,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=ψ1,∂H∂vx,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=fx,zψ2−ψ1−ψ3,∂H∂wx,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=−xψ1,∂H∂px,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=yψ3,∂H∂qx,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=−zψ3,∂H∂rx,y,z,ψ1,ψ2,ψ3,u,v,w,p,q,r=−yψ2.Herefx′x,z and fz′x,z are the partial derivatives of function fx,z with respect to variables x and z, respectively.According to Pontryagin Maximum Principle, for the optimal controlsu∗t, v∗t, w∗t, p∗t, q∗t, r∗t and the corresponding optimal solution x∗t,y∗t,z∗t of system (2.11), there exists a vector-function ψ∗t=ψ1∗t,ψ2∗t,ψ3∗t such that:(i)
The vector-functionψ∗t=ψ1∗t,ψ2∗t,ψ3∗t is a nontrivial solution of the adjoint system(3.3)
ψ
˙
1
∗
t
=
w
∗
t
ψ
1
∗
t
−
v
∗
t
f
x
′
x
∗
t
,
z
∗
t
ψ
2
∗
t
−
ψ
1
∗
t
−
ψ
3
∗
t
,
ψ
˙
2
∗
t
=
r
∗
t
ψ
2
∗
t
−
p
∗
t
ψ
3
∗
t
,
ψ
˙
3
∗
t
=
q
∗
t
ψ
3
∗
t
−
v
∗
t
f
z
′
x
∗
t
,
z
∗
t
ψ
2
∗
t
−
ψ
1
∗
t
−
ψ
3
∗
t
,
ψ
1
∗
T
=
0
,
ψ
2
∗
T
=
−
α
,
ψ
3
∗
T
=
−
β
.In this system, it is easy to note that(3.4)ψ˙1∗=−∂H∂xx∗,y∗,z∗,ψ1∗,ψ2∗,ψ3∗,u∗,v∗,w∗,p∗,q∗,r∗,ψ˙2∗=−∂H∂yx∗,y∗,z∗,ψ1∗,ψ2∗,ψ3∗,u∗,v∗,w∗,p∗,q∗,r∗,ψ˙3∗=−∂H∂zx∗,y∗,z∗,ψ1∗,ψ2∗,ψ3∗,u∗,v∗,w∗,p∗,q∗,r∗,and(3.5)ψ1∗T=−∂J∂xT,ψ2∗T=−∂J∂yT,ψ3∗T=−∂J∂zT.(ii)
The controlsu∗t, v∗t, w∗t, p∗t, q∗t, r∗t maximize the Hamiltonian(3.6)
H
x
∗
t
,
y
∗
t
,
z
∗
t
,
ψ
1
∗
t
,
ψ
2
∗
t
,
ψ
3
∗
t
,
u
,
v
,
w
,
p
,
q
,
r
,with respect to variables:
(3.7)
u
∈
u
min
,
u
max
,
v
∈
v
min
,
v
max
,
w
∈
w
min
,
w
max
,
p
∈
p
min
,
p
max
,
q
∈
q
min
,
q
max
,
r
∈
r
min
,
r
max
,for almost allt∈0,T, and, therefore, the following relationships hold:
(3.8)
u
∗
t
=
u
max
,
if
L
u
t
>
0
,
∀
u
∈
u
min
,
u
max
,
if
L
u
t
=
0
,
u
min
,
if
L
u
t
<
0
,
(3.9)
v
∗
t
=
v
max
,
if
L
v
t
>
0
,
∀
v
∈
v
min
,
v
max
,
if
L
v
t
=
0
,
v
min
,
if
L
v
t
<
0
,
(3.10)
w
∗
t
=
w
max
,
if
L
w
t
>
0
,
∀
w
∈
w
min
,
w
max
,
if
L
w
t
=
0
,
w
min
,
if
L
w
t
<
0
,
(3.11)
p
∗
t
=
p
max
,
if
L
p
t
>
0
,
∀
p
∈
p
min
,
p
max
,
if
L
p
t
=
0
,
p
min
,
if
L
p
t
<
0
,
(3.12)
q
∗
t
=
q
max
,
if
L
q
t
>
0
,
∀
q
∈
q
min
,
q
max
,
if
L
q
t
=
0
,
q
min
,
if
L
q
t
<
0
,
(3.13)
r
∗
t
=
r
max
,
if
L
r
t
>
0
,
∀
r
∈
r
min
,
r
max
,
if
L
r
t
=
0
,
r
min
,
if
L
r
t
<
0
.In these relationships, by Lemma3, functions(3.14)Lut=ψ1∗t,Lvt=ψ2∗t−ψ1∗t−ψ3∗t,Lwt=−ψ1∗t,Lpt=ψ3∗t,Lqt=−ψ3∗t,Lrt=−ψ2∗t,are the switching functions, which determine a type of the optimal controlsu∗t, v∗t, w∗t, p∗t, q∗t and r∗t.Please note that in (3.14)(3.15)Lut=−Lwt,Lpt=−Lqt.Hence, it suffices to study the properties of only two of these four functions, sayLwt and Lqt. Therefore, further we focus on the analysis of the switching functions Lvt, Lwt, Lqt and Lrt. By (3.3) and (3.14), we obtain the following Cauchy problem for the switching functions Lvt, Lwt, Lqt and Lrt:(3.16)L˙vt=v∗tfx′x∗t,z∗t+fz′x∗t,z∗t+r∗tLvt+w∗t−r∗tLwt+p∗t+q∗t−r∗tLqt,L˙wt=w∗tLwt+v∗tfx′x∗t,z∗tLvt,L˙qt=q∗tLqt+v∗tfx′x∗t,z∗tLvt,L˙rt=r∗tLrt−p∗tLqt,LvT=β−α,LwT=0,LqT=β,LrT=α.We actively use this Cauchy problem in subsequent arguments.


## 4. Properties of the Switching Functions

We start the analysis by establishing some properties of switching functionsLvt, Lwt, Lqt, Lrt.Lemma 4.
There exist valuestq,tr∈0,T such that
(4.1)
L
q
t
>
0
f
o
r
t
∈
t
q
,
T
;
L
r
t
>
0
f
o
r
t
∈
t
r
,
T
.
hold.Proof.
The switching functionsLqt and Lrt are the solutions of the third and fourth equations of system (3.16), respectively, and, hence, these functions are continuous functions. Since LqT=β, LrT=α and α,β>0, by Theorem on the Stability of the Sign of a Continuous Function (see, e.g., [62]) there are intervals tq,T and tr,T where inequalities (4.1) necessary hold.Lemma 5.
Depending on values of the weighting coefficients α and β, for switching functionsLvt and Lwt, the following statements hold:(i)
ifα<β, then there exist tv,tw∈0,T, tw<tv, such that(4.2)Lvt>0fort∈tv,T;Lwt<0fort∈tw,T;(ii)
ifα>β, then there exist tv,tw∈0,T, tw<tv, such that
(4.3)
L
v
t
<
0
f
o
r
t
∈
t
v
,
T
;
L
w
t
>
0
f
o
r
t
∈
t
w
,
T
.Proof.
One of the inequalities for functionLvt immediately follows from the arguments similar to those presented in the proof of Lemma 4. Indeed, switching function Lvt satisfies the first equation of system (3.16), and, hence, is absolutely continuous and, therefore, is also a continuous function. Then, for α≠β, by Theorem on the Stability of the Sign of a Continuous Function, there is tv 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 forLwt, we integrate the second equation of system (3.16) on the interval 0,T. The integration yields(4.4)Lwt=−∫tTe−∫tsw∗ξdξv∗sfx′x∗s,z∗sLvsds.
Then, taking into consideration the constraints imposed on controlv∗t by (2.5), Assumption 1 and the already established inequality for Lvt, we immediately see that the corresponding inequality for function Lwt in either (4.2) or (4.3) also necessary hold on the interval tv,T. From the absolute continuity of switching function Lwt and the Theorem on the Stability of the Sign of a Continuous Function it is immediately follows that there is value tw, tw<tv, such that, depending on the sign of α−β, the second inequality either in (4.2) or in (4.3) holds. This completes the proof.Corollary 6.
Lemmas4 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 u∗t, v∗t, w∗t, p∗t, q∗t and r∗t:(i)
ifα<β, then(4.5)u∗t=umax,w∗t=wminfort∈tw,T;v∗t=vmaxfort∈tv,T;p∗t=pmin,q∗t=qmaxfort∈tq,T;r∗t=rmaxfort∈tr,T;(ii)
ifα>β, then(4.6)u∗t=umin,w∗t=wmaxfort∈tw,T;v∗t=vminfort∈tv,T;p∗t=pmin,q∗t=qmaxfort∈tq,T;r∗t=rmaxfort∈tr,T.The following Lemma is of principal importance for further analysis.Lemma 7.
The switching functionsLvt, Lwt, Lqt and Lrt are not identically equal to zero on any subinterval of interval 0,T.Proof.
Firstly, we consider switching functionLvt. Let us suppose the contrary, that is that there is an interval Δv=t0,t1⊂0,T on which(4.7)Lvt=0.Then(4.8)L˙vt=0,almost everywhere on this interval as well and fort∈Δv equations of system (3.16) are(4.9)0=w∗t−r∗tLwt+p∗t+q∗t−r∗tLqt,L˙wt=w∗tLwt,L˙qt=q∗tLqt,L˙rt=r∗tLrt−p∗tLqt.Hence,(4.10)Lwt=Lw0e∫t0tw∗ξdξ,and(4.11)Lqt=Lq0e∫t0tq∗ξdξ,whereLw0, Lq0 are some constants.Now we have to consider the possible cases depending on the particular values ofLw0 and Lq0.Case 1. If Lw0=0 and Lq0=0, then, by (4.10) and (4.11),(4.12)Lwt=0,Lqt=0everywhere on intervalΔv. Substituting (4.7) and (4.12) into (3.14), we get(4.13)ψ1∗t=0,ψ2∗t=0,ψ3∗t=0everywhere on intervalΔv. System (3.3) is linear and homogeneous, and hence, equalities (4.13) necessary hold on the entire interval 0,T. This contradicts to the the boundary conditions at t=T and, hence, this case is impossible.Case 2. If Lw0=0 and Lq0≠0, then, from (4.10) and (4.11),(4.14)Lwt=0,Lqt≠0everywhere on intervalΔv. Hence, by (3.12) and (3.15), on Δv the controls p∗t and q∗t simultaneously take either the values pmin and qmax, or the values pmax and qmin. Moreover, by the inequality from (4.14) and the last equation in (4.9), on the interval Δv the switching function Lrt has at most one zero. (For analysis, it is sufficient that it is not equal to zero identically.) Hence, by (3.13), r∗∈rmin;rmax on Δv. Substituting (4.14) and the above-found values of p∗t and q∗t into the first equality in (4.9), we conclude that either(4.15)pmin+qmax−r∗=0,or(4.16)pmax+qmin−r∗=0necessary holds. However, either of these equalities contradicts to the control restrictions (2.9) and (2.10) of Assumption 2. Therefore, this case is also impossible.Case 3. If Lw0≠0 and Lq0=0, then, by (4.10) and (4.11),(4.17)Lwt≠0,Lqt=0everywhere on the intervalΔv. Hence, by (3.10), for all t∈Δv, the control w∗t takes either value wmin, or value wmax. Equality in (4.17) substituted into the last equation in (4.9) yields(4.18)Lrt=Lr0e∫t0tr∗ξdξ,whereLr0 is a constant. Here, the equality Lr0=0 is impossible, since, by (3.14), in the case we come to the contradictory equalities (4.13). Therefore, Lr0≠0, and then we have to conclude that on the interval Δv the control r∗t takes either the value rmin or value rmax. Substituting (4.17) and the above-found values of control w∗t into the first equality of (4.9), we immediately conclude that either(4.19)wmin=r∗,or(4.20)wmax=r∗must hold. However, neither of these equalities is possible as contradicting to the control restrictions (2.7) and (2.8) of Assumption 2. Therefore, this case is also impossible.Case 4. If Lw0≠0 and Lq0≠0, then, by (4.10) and (4.11),(4.21)Lwt≠0,Lqt≠0everywhere on intervalΔv. Hence, by (3.10) and (3.12), for all t∈Δv, control w∗t takes either value wmin, or value wmax, and controls p∗t and q∗t simultaneously take either values pmin and qmax, or values pmax and qmin. By the same arguments as in Case 2, we have to conclude that control r∗t takes either value rmin or value rmax on the entire interval Δv or on its subinterval Δv′.Differentiating the first equality of (4.9) on the interval Δv or on its subinterval Δv′, and then using the corresponding differential equations, we get equality(4.22)w∗w∗−r∗Lwt+q∗p∗+q∗−r∗Lqt=0,t∈Δvt∈Δv′,wherew∗∈wmin;wmax, p∗+q∗∈pmin+qmax;pmax+qmin and r∗∈rmin;rmax. The first equality in (4.9) and (4.22) form a linear homogeneous system with respect to unknown Lwt and Lqt. By (4.21), the solution of this system must be nontrivial, and, hence, the determinant of this system must be equal to zero. That is, the equality(4.23)w∗−r∗p∗+q∗−r∗w∗−q∗=0,must hold. However, by the control restrictions (2.7)–(2.10) of Assumption 2, this equality is impossible. Therefore, this case is also impossible.As all possible cases lead to a contradiction, we have to conclude that the hypothesis was incorrect and that the switching functionLvt cannot be identically equal to zero on any subinterval of the interval 0,T.Now we have to consider whether switching functionLwt can be equal to zero on a subinterval of the interval 0,T. Let us assume that there is an interval Δw⊂0,T on which(4.24)Lwt=0.Then, almost everywhere on this interval(4.25)L˙wt=0,holds as well. Substituting (4.24) and (4.25) into the second equation of system (3.16) and using the constraints on control v∗t from (2.5) and the properties of function fx,z from Assumption 1, we have to conclude that the contradictory equality Lvt=0 holds for all t∈Δw. Hence, our assumption is incorrect and the switching function Lwt cannot be identically equal to zero on any subinterval of the interval 0,T.Next, we consider a possibility that the switching functionLqt is identically equal to zero on any subinterval of the interval 0,T. Applying the same arguments as those for the function Lwt, but using the third equation of system (3.16) instead of the second equation, we come to the conclusion that the switching function Lqt cannot be identically equal to zero on any subinterval as well.Finally, we consider the switching functionLrt. Let us suppose the contrary and assume that there is the interval Δr⊂0,T on which(4.26)Lrt=0.Then, almost everywhere on this interval(4.27)L˙rt=0,holds. Again, we substitute (4.26) and (4.27) into the fourth equation of system (3.16). Using the constraints on control p∗t from (2.5) we find the contradictory equality Lqt=0 for t∈Δr. Hence the above-made assumption was incorrect and the switching function Lrt is not identically equal to zero on any subinterval of the interval 0,T. This completes the proof.Lemma7 combined with equalities (3.8)–(3.13) and (3.15) immediately leads to the following Corollary.Corollary 8.
The optimal controlsu∗t, v∗t, w∗t, p∗t, q∗t and r∗t are the bang-bang functions taking values umin;umax, vmin;vmax, wmin;wmax, pmin;pmax, qmin;qmax and rmin;rmax, respectively. Moreover, each of the pairs composed of controls u∗t and w∗t, and p∗t and q∗t switches from the minimum values to the maximum values and vice versa at the same switching moments, which are specific for each of these pairs.In addition to this corollary, we note that the controlsu∗t, v∗t, w∗t, p∗t, q∗t and r∗t have no singular arcs (see [23, 24] for details).


## 5. Estimating the Number of Zeros of the Switching Functions

In previous section, we proved that the optimal controlsu∗t, v∗t, w∗t, p∗t, q∗t and r∗t are the bang-bang functions with values umin;umax, vmin;vmax, wmin;wmax, pmin;pmax, qmin;qmax and rmin;rmax, 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 Lvt, Lwt, Lqt and Lrt 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 Lvt, then, applying the Generalized Rolle’s Theorem (see, e.g., [53]) to each of these equations, we can find estimates for the number of zeros of functions Lwt and Lqt as well. The use of this theorem with the fourth equation of system (3.16) gives us an estimate of the number of zeros of function Lrt. Thus, the analysis of switching function Lvt is the basis for completing our task. The method that allows to estimate the number of zeros of the function Lvt consists in the subsequent application of the above-mentioned Generalized Rolle’s Theorem to the equations of system (3.16).Let us consider a possibility of application of the Generalized Rolle’s Theorem to the first equation of system (3.16). The direct application of the theorem to this equation is impossible, because, apart of the term containing the function Lvt, the equation has two terms independent of Lvt. Therefore, we introduce an auxiliary function(5.1)Gvt=w∗t−r∗tLwt+p∗t+q∗t−r∗tLqt,which, by Lemma7 and relationships (2.7)–(2.10) of Assumption 2, is not equal to zero on any subinterval of 0,T. Then, the first equation of system (3.16) can be rewritten as(5.2)L˙vt=v∗tfx′x∗t,z∗t+fz′x∗t,z∗t+r∗tLvt+Gvt.In order to formulate the differential equation for functionGvt, we note that the switching functions Lwt and Lqt are absolutely continuous (as the solutions of system (3.16) and, therefore, are almost everywhere differentiable functions on the interval 0,T [63]. Theorem 3 of Chapter 4 in [60], gives us absolute continuity and hence, also almost everywhere differentiability of functions w∗tLwt and p∗t+q∗tLqt. Finally, arguments similar to these used by the authors in [34] allow to conclude that the functions w∗t−r∗t and p∗t+q∗t−r∗t are also differentiable almost everywhere on the interval 0,T. Then, using the second and third equations of system (3.16), we formulate the required differential equation:(5.3)
G
˙
v
t
=
v
∗
t
w
∗
t
−
r
∗
t
f
x
′
x
∗
t
,
z
∗
t
+
p
∗
t
+
q
∗
t
−
r
∗
t
f
z
′
x
∗
t
,
z
∗
t
L
v
t
+
w
∗
t
w
∗
t
−
r
∗
t
L
w
t
+
q
∗
t
p
∗
t
+
q
∗
t
−
r
∗
t
L
q
t
.Let us introduce another auxiliary function(5.4)Qvt=w∗tw∗t−r∗tLwt+q∗tp∗t+q∗t−r∗tLqt.Then, equation (5.3) can be rewritten as follows:(5.5)G˙vt=v∗tw∗t−r∗tfx′x∗t,z∗t+p∗t+q∗t−r∗tfz′x∗t,z∗tLvt+Qvt.By the above-presented arguments, functionQvt is also almost everywhere differentiable and satisfies(5.6)Q˙vt=v∗tw∗tw∗t−r∗tfx′x∗t,z∗t+q∗tp∗t+q∗t−r∗tfz′x∗t,z∗tLvt+w∗2tw∗t−r∗tLwt+q∗2tp∗t+q∗t−r∗tLqt.The last two terms of this equation can be expressed in terms of the auxiliary functionsGvt and Qvt. Then,(5.7)Q˙vt=v∗tw∗tw∗t−r∗tfx′x∗t,z∗t+q∗tp∗t+q∗t−r∗tfz′x∗t,z∗tLvt−w∗tq∗tGvt+w∗t+q∗tQvt.Let us denote(5.8)avt=v∗tfx′x∗t,z∗t+fz′x∗t,z∗t+r∗t,bvt=v∗tw∗t−r∗tfx′x∗t,z∗t+p∗t+q∗t−r∗tfz′x∗t,z∗t,cvt=v∗tw∗tw∗t−r∗tfx′x∗t,z∗t+q∗tp∗t+q∗t−r∗tfz′x∗t,z∗t,dvt=w∗tq∗t,evt=w∗t+q∗t.Due to the constraints on the controlsv∗t, w∗t, q∗t and r∗t in (2.5) and the properties of function fx,z in Assumption 1, we conclude that functions avt, dvt and evt are positive almost everywhere on the interval 0,T. Moreover, functions dvt and evt are piecewise-constant and functions avt, bvt and cvt are piecewise-continuous. Combining equations (5.2), (5.5) and (5.7), we write down the following system of differential equations for the switching function Lvt and auxiliary functions Gvt and Qvt, which is more convenient for the subsequent analysis than system (3.16):(5.9)L˙vt=avtLvt+Gvt,t∈0,T,G˙vt=bvtLvt+Qvt,Q˙vt=cvtLvt−dvtGvt+evtQvt.Let us rewrite system (5.9) in the following matrix form:(5.10)L˙vtG˙vtQ˙vt=avt10bvt01cvt−dvtevtLvtGvtQvt.Now, using a special change of variables, we reduce the matrix of system (5.10) to the upper-triangular form. This 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(5.11)L˜˙vt=avt−h21tL˜vt+G˜vt,G˜˙vt=h21t−h32tG˜vt+Q˜vt,Q˜˙vt=h32t+evtQ˜vt,whereL˜vt=Lvt is the switching function, and G˜vt, Q˜vt are the new auxiliary functions corresponding to L˜vt. According to the arguments of Appendix A, the functions h20t, h21t, h31t and h32t perform such a transformation and, by (A.21) and (A.26), satisfy the following systems of differential equations:(5.12)h˙20t=cvt+evt−avth20t+dvth21t+h20th21t,h˙21t=−bvt−h20t−avth21t+h212t,h200=h200=h320h210−h310,h210=h210,(5.13)h˙31t=−cvt+evt−avth31t−bvth32t+h31th32t,h˙32t=dvt−h31t+evth32t+h322t,h310=h310,h320=h320.Solutionsh20t, h21t of system (5.12) and solutions h31t, h32t of system (5.13) are determined on the maximal intervals of existence of such solutions [60, 64]: specifically, we assume that functions h20t and h21t are determined on the interval Δ2=t02,t12, and functions h31t and h32t are determined on the interval Δ3=t03,t13. Without loss of generality, we suppose that t12≤T and t13≤T.As we noted in AppendixA, we need to find the restrictions on the functions avt, bvt, cvt, dvt, evt and initial conditions h210, h310, h320, under which systems (5.12) and (5.13) have solutions defined on the entire interval 0,T. To do this, using the fact that dvt and evt are piecewise-constant functions, we transform the equations of systems (5.12) and (5.13) as follows:(5.14)ddth20t+dvt=−h20t+dvt−h21t+evt−avth20t+dvt+cvt−evt−avtdvt,ddt−h21t=−−h21t2+h20t+dvt−avt−h21t+bvt−dvt,(5.15)ddth31t−kvt=−h31t−kvt−h32t−0.5evt+0.5evt−avth31t−kvt+bvt−kvt−h32t−0.5evt−cvt−0.5bvt+kvtevt+avtkvt,ddt−h32t−0.5evt=−−h32t−0.5evt2+h31t−kvt.Herekvt=dvt−0.25ev2t is a piecewise-constant function as well, and, as it is not difficult to see, the inequality kvt<0 holds for almost all t∈0,T. Performing in these systems the change of variables(5.16)g20t=h20t+dvt,g21t=−h21t,g31t=h31t−kvt,g32t=−h32t−0.5evt,and introducing values(5.17)g210=−h210,g310=h310−kv+0,g320=−h320−0.5ev+0,whereev+0=limt→+0evt and kv+0=limt→+0kvt, we obtain the following systems of equations:(5.18)g˙20t=Φ20t,g20t,g21t=−g20tg21t+evt−avtg20t+cvt−evt−avtdvt,g˙21t=Φ21t,g20t,g21t=−g212t+g20t−avtg21t+bvt−dvt,g200=g320+ev+0g210−g310,g210=g210,(5.19)g˙31t=Φ31t,g31t,g32t=−g31tg32t+0.5evt−avtg31t+bvt−kvtg32t−cvt−0.5bvt+kvtevt+avtkvt,g˙32t=Φ32t,g31t,g32t=−g322t+g31t,g310=g310,g320=g320.Let us apply to systems (5.18) and (5.19) with nonnegative initial conditions(5.20)g320+ev+0g210−g310≥0,g210≥0,g310≥0,g320≥0,the necessary and sufficient condition for solutionsg20t, g21t and g31t, g32t to be nonnegative for all t∈0,t12 and t∈0,t13, respectively [65]. It consists in satisfying the so-called quasi-positivity condition(5.21)Φ20t,0,g21≥0,forallg21≥0,Φ21t,g20,0≥0,forallg20≥0,for system (5.18), and(5.22)Φ31t,0,g32≥0,forallg32≥0,Φ32t,g31,0≥0,forallg31≥0,for system (5.19).As a result, we obtain the inequalities:(5.23)bvt≥dvt,avtdvt+cvt≥dvtevt,avtkvt+cvt≤0.5bvt+kvtevt,which hold for almost allt∈0,T. These inequalities define the required restrictions on the functions avt, bvt, cvt, dvt and evt. Moreover, as it is easy to see choosing appropriate values g210, g310 and g320, inequalities (5.20) are satisfied as well.We consider the following assumption, the validity of which will be demonstrated in subsequent arguments.Assumption 9.
Let the system of inequalities(5.24)vtwt−rtfx′xt,zt+pt+qt−rtfz′xt,zt≥wtqt,wtqtvtfx′xt,zt+fz′xt,zt+rt+vtwtwt−rtfx′xt,zt+qtpt+qt−rtfz′xt,zt≥wtqtwt+qt,wt−qt2vtfx′xt,zt+fz′xt,zt−0.5wt−rt+qt−rt+4vt1−2wtwt−rtfx′xt,zt+1−2qtpt+qt−rtfz′xt,zt≥0
hold almost everywhere on the interval0,T for arbitrary admissible controls ut, vt, wt, pt, qt and rt and the components xt, yt and zt of the corresponding solution xt,yt,zt to system (2.11).By definitions of functionsavt, bvt, cvt, dvt, evt and kvt, it is easy to see that under Assumption 9 the inequalities (5.23) hold. The verification of this assumption for the case of a bilinear function fx,z=κxz, where κ is a positive constant, is discussed in Appendix B.Finally, we have to verify the continuation of the nonnegative solutionsg20t, g21t, g31t and g32t of systems (5.18) and (5.19) to the entire interval 0,T (see Theorem 1.6 in [65]). To do this, we define Lyapunov functions:(5.25)V2g20,g21=0.5g202+g212,V3g31,g32=0.5g312+g322,and estimate the derivativesdV2/dt and dV3/dt of these functions along the trajectories of the systems of differential equations (5.18) and (5.19), respectively. (The derivatives are the scalar products of the gradients of these functions and the right-hand sides of the corresponding systems.) The Lyapunov functions satisfy(5.26)dV2dtt=−g21tg202t+g212t+evt−avtg202t+g20tg21t−avtg212t+cvt−evt−avtdvtg20t+bvt−dvtg21t+cvt−evt−avtdvtg20t+bvt−dvtg21t≤K2V2g20t,g21t,(5.27)dV3dtt=−g32tg312t+g322t+0.5evt−avtg312t+bvt−kvt+1g31tg32t−cvt−0.5bvt+kvtevt+avtkvtg31t≤K3V3g31t,g32t.Here,K2 and K3 are constants defined as following:(5.28)K2=2maxt∈0,Tavt+evt;1;bvt+dvt;cvt+avt+evtdvt,(5.29)K3=2maxt∈0,Tavt+0.5evt;bvt+kvt+1;cvt+0.5bvt+kvtevt+avtkvt.We have to note that these inequalities are obtained under constraints(5.30)g202t+g212t≥1,g312t+g322t≥1,which follow from the relationships(5.31)limt→t12−0g20t,g21t⊤=+∞,limt→t13−0g31t,g32t⊤=+∞(see Lemma in Section 14, Chapter 4 in [66]). That is, the estimates(5.32)dV2dtt≤K2V2g20t,g21t,dV3dtt≤K3V3g31t,g32t,imply the required continuation of the nonnegative solutionsg20t and g21t of system (5.18) and the nonnegative solutions g31t and g32t of system (5.19) to the entire interval 0,T. Taking into consideration that functions h20t, h21t, h31t and h32t satisfy systems (5.12) and (5.13) and are also defined on the entire interval 0,T, we conclude that system (5.11) is defined on this interval as well.Using system (5.11), we can estimate the number of zeros of the switching function Lvt=L˜vt. Taking into consideration that LvT=β−α≠0 (see (3.16)), we will show that this switching function has at most two distinct zeros on the interval 0,T.Let us suppose the contrary, that is that the functionLvt=L˜vt has at least three distinct zeroes on this interval. Then, by virtue of the Generalized Rolle’s Theorem [53] applied to the first equation of system (5.11), the function G˜vt has at least two distinct zeroes on the interval 0,T. Applying this theorem to the second equation of system (5.11), we have to conclude that the function Q˜vt has at least one zero on the interval 0,T. However, this function satisfies the linear homogeneous equation, and, therefore, if Q˜vt=0 at a single point of interval 0,T, then Q˜vt=0 on the entire interval 0,T. Moreover, in this case the second equation of system (5.11) is also the linear homogeneous equation, and, hence, if function G˜vt has at least two zeros on interval 0,T, than G˜vt=0 everywhere on the interval 0,T as well. Finally, applying the same line of arguments to the first equation of system (5.11), we have to conclude that in this case L˜vt=0 everywhere on the entire interval 0,T. However, this contradicts to Lemma 7 and, hence, the assumption is incorrect. Thus, the required estimate for the number of zeros of the switching function Lvt=L˜vt is established.Denote the zeros of functionLvt by τ1∗ and τ2∗; τ1∗,τ2∗∈0,T. Depending on the sign of the difference β−α, the function Lvt is described by one of the following two relationships, namely:(i)
Ifβ−α>0, then(5.33)Lvt=>0,if0≤t<τ1∗,=0,ift=τ1∗,<0,ifτ1∗<t<τ2∗,=0,ift=τ2∗,>0,ifτ2∗<t≤T;(ii)
Ifβ−α<0, then(5.34)
L
v
t
=
<
0
,
if
0
≤
t
<
τ
1
∗
,
=
0
,
if
t
=
τ
1
∗
,
>
0
,
if
τ
1
∗
<
t
<
τ
2
∗
,
=
0
,
if
t
=
τ
2
∗
,
<
0
,
if
τ
2
∗
<
t
≤
T
.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 functionLwt. We will show that, taking into consideration that LwT=0 (see (3.16)), function Lwt has at most three distinct zeros on the interval 0,T. To prove this preposition, we, again, suppose the opposite: that is, we suppose that function Lwt has at least four distinct zeroes on interval 0,T. Applying the Generalized Rolle’s Theorem, we have to conclude that in this case function Lvt has at least three distinct zeroes on the interval 0,T, which contradicts to the above-found estimates for the number of zeros of switching function Lvt. Therefore, our assumption was incorrect, and the required estimate for the number of zeros of the switching function Lwt is found.SinceLwT=0 (see (3.16)), one of these three zeros is the end point of the interval 0,T. We denote the other two zeros by θ1∗ and θ2∗; θ1∗,θ2∗∈0,T. Then, taking into consideration the corresponding inequalities from Lemma 5, we conclude that, depending on the sign of β−α, the switching function Lwt is described by one of the following two relationships, namely:(ii)
ifβ−α>0, then(5.35)Lwt=<0,if0≤t<θ1∗,=0,ift=θ1∗,>0,ifθ1∗<t<θ2∗,=0,ift=θ2∗,<0,ifθ2∗<t<T,=0,ift=T;(ii)
Ifβ−α<0, then(5.36)
L
w
t
=
>
0
,
if
0
≤
t
<
θ
1
∗
,
=
0
,
if
t
=
θ
1
∗
,
<
0
,
if
θ
1
∗
<
t
<
θ
2
∗
,
=
0
,
if
t
=
θ
2
∗
,
>
0
,
if
θ
2
∗
<
t
<
T
,
=
0
,
if
t
=
T
.To estimate the number of zeros of the switching functionLqt, we take into consideration the boundary condition LqT=β>0 (see (3.16)). Applying the arguments that are analogous to those used to estimate the number of zeros of function Lwt, we can prove that the function Lqt has at most three distinct zeros on the interval 0,T. 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 Lqt satisfies(5.37)Lqt=<0,if0≤t<η1∗,=0,ift=η1∗,>0,ifη1∗<t<η2∗,=0,ift=η2∗,<0,ifη2∗<t<η3∗,=0,ift=η3∗,>0,ifη3∗<t≤T.Hereη1∗,η2∗,η3∗∈0,T are zeros of the switching function Lqt.Please note that, by (3.15), relationships (5.35)–(5.37) also describe the behavior of the switching functions Lut and Lpt.Finally, we estimate the number of zeros of switching functionLrt. Taking into consideration the boundary condition LrT=α>0 (see (3.16)), we show that function Lrt has at most four distinct zeros on interval 0,T. Let us prove this by contradiction: assume that the function Lrt has at least five distinct zeros on the interval 0,T. Applying the Generalized Rolle’s Theorem, we conclude that in this case the function Lqt has at least four distinct zeros on the interval 0,T. This conclusion contradicts to the above-established estimate, and, therefore, our assumption is incorrect.Denoting the zeros of functionLrt by χ1∗, χ2∗, χ3∗ and χ4∗; χ1∗,χ2∗,χ3∗,χ4∗∈0,T, and taking into consideration the second inequality from (4.1) in Lemma 4, we get the following relationship describing behavior of the switching function Lrt:(5.38)Lrt=>0,if0≤t<χ1∗,=0,ift=χ1∗,<0,ifχ1∗<t<χ2∗,=0,ift=χ2∗,>0,ifχ2∗<t<χ3∗,=0,ift=χ3∗,<0,ifχ3∗<t<χ4∗,=0,ift=χ4∗,>0,ifχ4∗<t≤T.


## 6. Types of the Optimal Controls

The 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 v∗t, w∗t, q∗t and r∗t. Specifically, depending on the sign of β−α, the optimal controls v∗t and w∗t are of one of the following types:(i)
Ifβ−α>0, then(6.1)v∗t=vmax,if0≤t≤τ1∗,vmin,ifτ1∗<t≤τ2∗,vmax,ifτ2∗<t≤T,w∗t=wmin,if0≤t≤θ1∗,wmax,ifθ1∗<t≤θ2∗,wmin,ifθ2∗<t≤T.(ii)
Ifβ−α<0, then(6.2)
v
∗
t
=
v
min
,
if
0
≤
t
≤
τ
1
∗
,
v
max
,
if
τ
1
∗
<
t
≤
τ
2
∗
,
v
min
,
if
τ
2
∗
<
t
≤
T
,
w
∗
t
=
w
max
,
if
0
≤
t
≤
θ
1
∗
,
w
min
,
if
θ
1
∗
<
t
≤
θ
2
∗
,
w
max
,
if
θ
2
∗
<
t
≤
T
.Here,τ1∗,τ2∗∈0,T and θ1∗,θ2∗∈0,T are the switchings of the optimal controls v∗t and w∗t, respectively.The optimal controlq∗t is of the following type:(6.3)q∗t=qmin,if0≤t≤η1∗,qmax,ifη1∗<t≤η2∗,qmin,ifη2∗<t≤η3∗,qmax,ifη3∗<t≤T,whereη1∗,η2∗,η3∗∈0,T are the moments of switching.The optimal controlr∗t is described as following:(6.4)r∗t=rmax,if0≤t≤χ1∗,rmin,ifχ1∗<t≤χ2∗,rmax,ifχ2∗<t≤χ3∗,rmin,ifχ3∗<t≤χ4∗,rmax,ifχ4∗<t≤T,whereχ1∗,χ2∗,χ3∗,χ4∗∈0,T are the moments of switching.Taking into consideration (3.15), we can immediately describe the type of controls u∗t and p∗t as well. Specifically, depending on the sign of the difference β−α, the optimal control u∗t belongs to one of the following two types:(i)
Ifβ−α>0, then(6.5)u∗t=umax,if0≤t≤θ1∗,umin,ifθ1∗<t≤θ2∗,umax,ifθ2∗<t≤T.(ii)
Ifβ−α<0, then(6.6)
u
∗
t
=
u
min
,
if
0
≤
t
≤
θ
1
∗
,
u
max
,
if
θ
1
∗
<
t
≤
θ
2
∗
,
u
min
,
if
θ
2
∗
<
t
≤
T
.The optimal controlp∗t is of the following qualitative type:(6.7)p∗t=pmax,if0≤t≤η1∗,pmin,ifη1∗<t≤η2∗,pmax,ifη2∗<t≤η3∗,pmin,ifη3∗<t≤T.


## 7. 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 software 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 open-source software package IPOPT designed for large-scale nonlinear optimization.In our calculations, we use a time grid with5000 nodes over time interval 0,T that is equivalent to the time step of Δt=2T·10−4. Since the problem is solved by a direct method using an iterative approach, at each step we impose the acceptable convergence tolerance εrel=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:(7.1)λ=0.13μ=0.013σ=0.4k=0.4γ=0.05x0=5.0y0=0.15z0=20.0α=1.0β=0.01umin=0.10umax=0.13vmin=0.8vmax=1.0wmin=0.013wmax=0.018pmin=0.1pmax=0.4qmin=0.05qmax=0.08rmin=0.4rmax=0.6These values are adopted from [68]. We note that α>β.Results of the numerical calculations forT=60.0 are presented in Figures 1–4. Figures 1 and 2 correspond to the outcomes for the following bilinear incidence rate(7.2)fx,z=κxz.Figure 1.
Optimal solutions forT=60.0: x∗t, y∗t, z∗t.Figure 2.
Optimal controls forT=60.0: upper row: u∗t, v∗t, w∗t; lower row: p∗t, q∗t, r∗t.withκ=0.00216. Figure 1 show the optimal solutions x∗t, y∗t and z∗t, while Figure 2 display the corresponding optimal controls u∗t, v∗t, w∗t, p∗t, q∗t and r∗t. J∗=0.270579 is the minimum value of the functional Ju,v,w,p,q,r (see (2.15)). IPOPT used 40 iterations in order to solve the optimal control problem (2.11) and (2.15).Figures3 and 4 present results for a nonlinear incidence rate fx,z of the Michaelis–Menten kinetics type with saturation in both variables:(7.3)fx,z=κx1+mxz1+lz,Figure 3.
Optimal solutions forT=60.0: x∗t, y∗t, z∗t.Figure 4.
Optimal controls forT=60.0: upper row: u∗t, v∗t, w∗t; lower row: p∗t, q∗t, r∗t.whereκ=0.00216, m=0.1, l=0.001. Figure 3 show the optimal solutions x∗t, y∗t and z∗t, while Figure 4 display the corresponding optimal controls u∗t, v∗t, w∗t, p∗t, q∗t and r∗t. J∗=0.250366 is the minimum value of the functional Ju,v,w,p,q,r of (2.15). IPOPT used 45 iterations in order to solve the optimal control problem (2.11) and (2.15).Results of the numerical calculations forT=120.0 are presented in Figures 5–8. Figures 5 and 6 correspond to the outcomes for the bilinear incidence rate fx,z of (7.2). Figure 5 shows the optimal solutions x∗t, y∗t and z∗t, while Figure 5 displays the corresponding optimal controls u∗t, v∗t, w∗t, p∗t, q∗t and r∗t. J∗=1.249370·10−3 is the minimum value of the functional Ju,v,w,p,q,r (see (2.15)). IPOPT used 50 iterations in order to solve the optimal control problem (2.11), (2.15).Figure 5.
Optimal solutions forT=120.0: x∗t, y∗t, z∗t.Figure 6.
Optimal controls forT=120.0: upper row: u∗t, v∗t, w∗t; lower row: p∗t, q∗t, r∗t.Figures7 and 8 present results for the nonlinear incidence rate fx,z of (7.3). Figure 7 show the optimal solutions x∗t, y∗t and z∗t, while Figure 8 display the corresponding optimal controls u∗t, v∗t, w∗t, p∗t, q∗t and r∗t. J∗=1.425958·10−3 is the minimum value of the functional Ju,v,w,p,q,r of (2.15). IPOPT used 55 iterations in order to solve the optimal control problem (2.11) and (2.15).Figure 7.
Optimal solutions forT=120.0: x∗t, y∗t, z∗t.Figure 8.
Optimal controls forT=120.0: upper row: u∗t, v∗t, w∗t; lower row: p∗t, q∗t, r∗t.The numerical calculations lead to the following conclusions:(i)
A change in the type of the incidence ratefx,z does not lead to qualitative changes in the behavior of the optimal solutions x∗t, y∗t and z∗t.(ii)
Taking into account the relationship (3.15) between the switching functions Lut, Lwt and Lpt, Lqt, it can be seen that the optimal controls u∗t, v∗t and w∗t are piecewise constant functions with switchings. At the same time, the optimal controls p∗t, q∗t and r∗t are the constant functions over the whole interval 0,T. Moreover, all optimal controls behave at the end of the time interval precisely in accordance with formulas (6.2)–(6.4), (6.6) and (6.7).


## 8. 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 functionfx,z constrained by a few biologically feasible conditions and takes account of the loss of free virus due to the target cell infections. The constraints imposed on incidence rate fx,z summarized in Assumption 1 are biologically motivated and hold for any realistic incidence rate. These 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–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. These six controls comprise all controls possible for this model and represent the actions that the existing anti-HIV drugs exhibit. The optimal control 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. The realization of a specific variant depends on specific initial conditionsx0,y0,z0 and the maximal values of derivatives fx′x,z and fz′x,z 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 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 of60 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.(i)
The first group is composed of controlsu˜t (“immunization” of CD4+ T cells before they enter patient’s bloodstream) and w˜t (“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)
The second group combines controlsp˜t (inhibition of the virus production by the infected cells) and q˜t (immobilization or removal of the free virus particles). The first of these controls correspond to actions of the protease and maturation inhibitors. Drugs represented by control q˜t are not currently available, but they can appear in the nearest future. Action of such drugs is similar to that of antibodies.(iii)
The third and the fourth groups contain a single control each; namely,r˜t (killing or immobilizing of the infected cells, or removing them by other means) and v˜t (inhibiting the infection rate). Control v˜t corresponds to action of NRTI, NtRTI, NNRTI and the integrase inhibitors. Drugs represented by control r˜t do not exist, but such drugs can be developed in the next few years.It is remarkable that in the absence of controlv˜t the drugs of the first, second and third groups are completely independent and do not interfere. This implies that in the absence of control v˜t, 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 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.



---

*Source: 8107106-2020-02-01.xml*