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 the -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.
1. Introduction
1.1. Antiretroviral Drugs and Guidelines
The World Health Organization and the United Nations Programme on HIV/AIDS (UNAIDS) estimated that over millions of people were living with HIV/AIDS in 2017; other million were newly infected with HIV in the same year [1, 2]. Of these millions of the infected, millions were living on antiretroviral therapy [1, 2]. These figures mean that the therapy is currently covers less than 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 $ per patient per year in the United States (where costs of the therapy are the highest) [3] to US $– 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 percent of infections and HIV-associated deaths occurred.
Starting from 1987, when the first antiretroviral drugs, Nucleoside and Nucleotide Reverse Transcriptase Inhibitors, were approved for HIV treatment, over 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 -dimensional antiretroviral therapy model. The model in this paper was based on the -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 -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 -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 -dimensional system of nonautonomous quadratic differential equations to two -dimensional systems: the analysis of these two -dimensional systems is simpler than that of a single -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 Section 2, we formulate the model and state the optimal control problem. In Section 3, we apply Pontryagin Maximum Principle to find the adjoint system and the system of nonautonomous differential equations for the switching functions. In Sections 4, we prove that the optimal controls are bang-bang. In Section 5, we find the estimates for the number of zeros of the switching functions, and thereby, in Section 6, we describe the principal types of the optimal controls that are possible for this problem. In Section 7, we provide results of computations that illustrate our analytical results. Section 8 contains conclusion and discussion of the results. Appendix A describes details of the technique of the reduction of a -dimensional matrix of a system of linear nonautonomous homogeneous differential equations for the switching functions to the upper triangular form. Appendix B is devoted to the verification of one of our auxiliary assumptions.
2. Model and the Optimal Control Problem
We consider the following model of the within-host HIV dynamics:
Here , and are the appropriate concentrations of uninfected (susceptible) target cells, infected target cells and free virus particles at time . 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 and the infected cells immediately move to the infected population, that the infected cells produce free virus particles at the rate , and that average life-spans of the susceptible cells, infected cells and free virus particles are, respectively, , and (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 . 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 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 -dimensional space is a positive invariant set of system (2.1) and that all solutions of the system initiated in this set are bounded.
In this paper, we assume that the incidence rate satisfies the following conditions.
Assumption 1. Let the function satisfy:(i) is a continuous differentiable function for all ;(ii) for all , and for all ;(iii)Partial derivatives and are positive for all and nonnegative for all and .
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 the system is an object of six independent controls , , , , and , which have the following biological meanings:(i) is a fraction of immune cells among the new target cells entering the system;(ii) is a per capita rate of “immunization” of the susceptible cells that makes the “immunized” cells entirely unsusceptible to the virus;(iii) is a magnitude of indirect measures that inhibit the incidence rate;(iv) is a magnitude of indirect measures aimed that inhibit the free virus particles production by the infected cells;(v) is a per capita rate of immobilization (or removal) of free virus particles from the system;(vi) is a per capita rate of killing (or of removal by other means) of the infected cells.
We would like to stress that these six controls comprise all interventions that are feasible or potentially possible for model (2.1).
We assume that the controls satisfy the following constraints:
Here the maximal values of control functions , , , , and are pre-defined positive constants.
Assuming that these six controls are applied to model (2.1) on the time interval , we obtain the following control model:
For convenience of notation, we introduce modified controls:
These modified controls satisfy the following constraints:
where
In addition to Assumption 1, we make the following assumption.
Assumption 2. Maximal and minimal levels , , , , , , , , of the controls satisfy the conditions:
As a set of admissible controls, , we consider the set of all Lebesque measurable functions , , , , and satisfying inequalities (2.5) for almost all .
Using modified controls (2.4), we rewrite control system (2.3) as following:
The following lemma ensures the boundedness, positiveness, and continuity of solutions for system (2.11).
Lemma 3. For any admissible controls , , ,, , and , the corresponding absolutely continuous solution , , to system (2.11) are defined on the entire interval , and for its components , , the inequalities:wherehold for all .
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 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 :
Such objective functional usually corresponds to the task of complete elimination of the infection. In (2.15), and are positive weights, such that . Please note that the objective functional is independent of the controls.
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 , , , , and and the corresponding optimal solution 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
where , and are adjoint variables. The Hamiltonian satisfies:
Here and are the partial derivatives of function with respect to variables and , respectively.
According to Pontryagin Maximum Principle, for the optimal controls , , , , , and the corresponding optimal solution of system (2.11), there exists a vector-function such that:(i)The vector-function is a nontrivial solution of the adjoint system
In this system, it is easy to note that
and(ii)The controls , , , , , maximize the Hamiltonian
with respect to variables:
for almost all , and, therefore, the following relationships hold:
In these relationships, by Lemma 3, functions
are the switching functions, which determine a type of the optimal controls , , , , and .
Please note that in (3.14)
Hence, it suffices to study the properties of only two of these four functions, say and . Therefore, further we focus on the analysis of the switching functions , , and . By (3.3) and (3.14), we obtain the following Cauchy problem for the switching functions , , and :
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 functions , , , .
Lemma 4. There exist values such that
hold.
Proof. The switching functions and are the solutions of the third and fourth equations of system (3.16), respectively, and, hence, these functions are continuous functions. Since , and , by Theorem on the Stability of the Sign of a Continuous Function (see, e.g., [62]) there are intervals and where inequalities (4.1) necessary hold.
Lemma 5. Depending on values of the weighting coefficients α and β, for switching functions and , the following statements hold:(i)if , then there exist , , such that(ii)if , then there exist , , such that
Proof. One of the inequalities for function immediately follows from the arguments similar to those presented in the proof of Lemma 4. Indeed, switching function satisfies the first equation of system (3.16), and, hence, is absolutely continuous and, therefore, is also a continuous function. Then, for , by Theorem on the Stability of the Sign of a Continuous Function, there is such that, depending on the sign of , either the first inequality in (4.2), or the first inequality in (4.3), holds.
To get the inequality for , we integrate the second equation of system (3.16) on the interval . The integration yieldsThen, taking into consideration the constraints imposed on control by (2.5), Assumption 1 and the already established inequality for , we immediately see that the corresponding inequality for function in either (4.2) or (4.3) also necessary hold on the interval . From the absolute continuity of switching function and the Theorem on the Stability of the Sign of a Continuous Function it is immediately follows that there is value , , such that, depending on the sign of , the second inequality either in (4.2) or in (4.3) holds. This completes the proof.
Corollary 6. Lemmas 4 and 5 combined with formulas (3.9), (3.10), (3.12), (3.13) and (3.15), immediately lead to the following relationships for the optimal controls , , , , and :(i)if , then(ii)if , then
The following Lemma is of principal importance for further analysis.
Lemma 7. The switching functions , , and are not identically equal to zero on any subinterval of interval .
Proof. Firstly, we consider switching function . Let us suppose the contrary, that is that there is an interval on which
Then
almost everywhere on this interval as well and for equations of system (3.16) are
Hence,
and
where , are some constants.
Now we have to consider the possible cases depending on the particular values of and .
Case 1. If and , then, by (4.10) and (4.11),
everywhere on interval . Substituting (4.7) and (4.12) into (3.14), we get
everywhere on interval . System (3.3) is linear and homogeneous, and hence, equalities (4.13) necessary hold on the entire interval . This contradicts to the the boundary conditions at and, hence, this case is impossible.
Case 2. If and , then, from (4.10) and (4.11),
everywhere on interval . Hence, by (3.12) and (3.15), on the controls and simultaneously take either the values and , or the values and . Moreover, by the inequality from (4.14) and the last equation in (4.9), on the interval the switching function has at most one zero. (For analysis, it is sufficient that it is not equal to zero identically.) Hence, by (3.13), on . Substituting (4.14) and the above-found values of and into the first equality in (4.9), we conclude that either
or
necessary 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 and , then, by (4.10) and (4.11),
everywhere on the interval . Hence, by (3.10), for all , the control takes either value , or value . Equality in (4.17) substituted into the last equation in (4.9) yields
where is a constant. Here, the equality is impossible, since, by (3.14), in the case we come to the contradictory equalities (4.13). Therefore, , and then we have to conclude that on the interval the control takes either the value or value . Substituting (4.17) and the above-found values of control into the first equality of (4.9), we immediately conclude that either
or
must hold. However, neither of these equalities is possible as contradicting to the control restrictions (2.7) and (2.8) of Assumption 2. Therefore, this case is also impossible.
Case 4. If and , then, by (4.10) and (4.11),
everywhere on interval . Hence, by (3.10) and (3.12), for all , control takes either value , or value , and controls and simultaneously take either values and , or values and . By the same arguments as in Case 2, we have to conclude that control takes either value or value on the entire interval or on its subinterval .
Differentiating the first equality of (4.9) on the interval or on its subinterval , and then using the corresponding differential equations, we get equality
where , and . The first equality in (4.9) and (4.22) form a linear homogeneous system with respect to unknown and . By (4.21), the solution of this system must be nontrivial, and, hence, the determinant of this system must be equal to zero. That is, the equality
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 function cannot be identically equal to zero on any subinterval of the interval .
Now we have to consider whether switching function can be equal to zero on a subinterval of the interval . Let us assume that there is an interval on which
Then, almost everywhere on this interval
holds as well. Substituting (4.24) and (4.25) into the second equation of system (3.16) and using the constraints on control from (2.5) and the properties of function from Assumption 1, we have to conclude that the contradictory equality holds for all . Hence, our assumption is incorrect and the switching function cannot be identically equal to zero on any subinterval of the interval .
Next, we consider a possibility that the switching function is identically equal to zero on any subinterval of the interval . Applying the same arguments as those for the function , but using the third equation of system (3.16) instead of the second equation, we come to the conclusion that the switching function cannot be identically equal to zero on any subinterval as well.
Finally, we consider the switching function . Let us suppose the contrary and assume that there is the interval on which
Then, almost everywhere on this interval
holds. Again, we substitute (4.26) and (4.27) into the fourth equation of system (3.16). Using the constraints on control from (2.5) we find the contradictory equality for . Hence the above-made assumption was incorrect and the switching function is not identically equal to zero on any subinterval of the interval . This completes the proof.
Lemma 7 combined with equalities (3.8)–(3.13) and (3.15) immediately leads to the following Corollary.
Corollary 8. The optimal controls , , , , and are the bang-bang functions taking values , , , , and , respectively. Moreover, each of the pairs composed of controls and , and and 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 controls , , , , and 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 controls , , , , and are the bang-bang functions with values , , , , and , respectively. In this section, we proceed to estimating the number of switchings of these controls. Equalities (3.8)–(3.13) and (3.15) show that for this task it suffices to estimate the number of zeros of the switching functions , , and defined by system (3.16). From the analysis of the second and third equations of this system, we conclude that if we have an estimate for the number of zeros of function , then, applying the Generalized Rolle’s Theorem (see, e.g., [53]) to each of these equations, we can find estimates for the number of zeros of functions and 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 . Thus, the analysis of switching function is the basis for completing our task. The method that allows to estimate the number of zeros of the function 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 , the equation has two terms independent of . Therefore, we introduce an auxiliary function
which, by Lemma 7 and relationships (2.7)–(2.10) of Assumption 2, is not equal to zero on any subinterval of . Then, the first equation of system (3.16) can be rewritten as
In order to formulate the differential equation for function , we note that the switching functions and are absolutely continuous (as the solutions of system (3.16) and, therefore, are almost everywhere differentiable functions on the interval [63]. Theorem 3 of Chapter 4 in [60], gives us absolute continuity and hence, also almost everywhere differentiability of functions and . Finally, arguments similar to these used by the authors in [34] allow to conclude that the functions and are also differentiable almost everywhere on the interval . Then, using the second and third equations of system (3.16), we formulate the required differential equation:
Let us introduce another auxiliary function
Then, equation (5.3) can be rewritten as follows:
By the above-presented arguments, function is also almost everywhere differentiable and satisfies
The last two terms of this equation can be expressed in terms of the auxiliary functions and . Then,
Let us denote
Due to the constraints on the controls , , and in (2.5) and the properties of function in Assumption 1, we conclude that functions , and are positive almost everywhere on the interval . Moreover, functions and are piecewise-constant and functions , and 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 and auxiliary functions and , which is more convenient for the subsequent analysis than system (3.16):
Let us rewrite system (5.9) in the following matrix form:
Now, using a special change of variables, we reduce the matrix of system (5.10) to the upper-triangular form. 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
where is the switching function, and , are the new auxiliary functions corresponding to . According to the arguments of Appendix A, the functions , , and perform such a transformation and, by (A.21) and (A.26), satisfy the following systems of differential equations:
Solutions , of system (5.12) and solutions , of system (5.13) are determined on the maximal intervals of existence of such solutions [60, 64]: specifically, we assume that functions and are determined on the interval , and functions and are determined on the interval . Without loss of generality, we suppose that and .
As we noted in Appendix A, we need to find the restrictions on the functions , , , , and initial conditions , , , under which systems (5.12) and (5.13) have solutions defined on the entire interval . To do this, using the fact that and are piecewise-constant functions, we transform the equations of systems (5.12) and (5.13) as follows:
Here is a piecewise-constant function as well, and, as it is not difficult to see, the inequality holds for almost all . Performing in these systems the change of variables
and introducing values
where and , we obtain the following systems of equations:
Let us apply to systems (5.18) and (5.19) with nonnegative initial conditions
the necessary and sufficient condition for solutions , and , to be nonnegative for all and , respectively [65]. It consists in satisfying the so-called quasi-positivity condition
for system (5.18), and
for system (5.19).
As a result, we obtain the inequalities:
which hold for almost all . These inequalities define the required restrictions on the functions , , , and . Moreover, as it is easy to see choosing appropriate values , and , 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 inequalitieshold almost everywhere on the interval for arbitrary admissible controls , , , , and and the components , and of the corresponding solution to system (2.11).
By definitions of functions , , , , and , 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 , where is a positive constant, is discussed in Appendix B.
Finally, we have to verify the continuation of the nonnegative solutions , , and of systems (5.18) and (5.19) to the entire interval (see Theorem 1.6 in [65]). To do this, we define Lyapunov functions:
and estimate the derivatives and 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
Here, and are constants defined as following:
We have to note that these inequalities are obtained under constraints
which follow from the relationships
(see Lemma in Section 14, Chapter 4 in [66]). That is, the estimates
imply the required continuation of the nonnegative solutions and of system (5.18) and the nonnegative solutions and of system (5.19) to the entire interval . Taking into consideration that functions , , and satisfy systems (5.12) and (5.13) and are also defined on the entire interval , 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 . Taking into consideration that (see (3.16)), we will show that this switching function has at most two distinct zeros on the interval .
Let us suppose the contrary, that is that the function 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 has at least two distinct zeroes on the interval . Applying this theorem to the second equation of system (5.11), we have to conclude that the function has at least one zero on the interval . However, this function satisfies the linear homogeneous equation, and, therefore, if at a single point of interval , then on the entire interval . Moreover, in this case the second equation of system (5.11) is also the linear homogeneous equation, and, hence, if function has at least two zeros on interval , than everywhere on the interval 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 everywhere on the entire interval . 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 is established.
Denote the zeros of function by and ; . Depending on the sign of the difference , the function is described by one of the following two relationships, namely:(i)If , then(ii)If , then
Please note that these formulas completely agree with the corresponding inequalities from (4.2) and (4.3) of Lemma 5.
Now we proceed to estimating the number of zeros of the switching function . We will show that, taking into consideration that (see (3.16)), function has at most three distinct zeros on the interval . To prove this preposition, we, again, suppose the opposite: that is, we suppose that function has at least four distinct zeroes on interval . Applying the Generalized Rolle’s Theorem, we have to conclude that in this case function has at least three distinct zeroes on the interval , which contradicts to the above-found estimates for the number of zeros of switching function . Therefore, our assumption was incorrect, and the required estimate for the number of zeros of the switching function is found.
Since (see (3.16)), one of these three zeros is the end point of the interval . We denote the other two zeros by and ; . Then, taking into consideration the corresponding inequalities from Lemma 5, we conclude that, depending on the sign of , the switching function is described by one of the following two relationships, namely:(ii)if , then(ii)If , then
To estimate the number of zeros of the switching function , we take into consideration the boundary condition (see (3.16)). Applying the arguments that are analogous to those used to estimate the number of zeros of function , we can prove that the function has at most three distinct zeros on the interval . We omit the arguments, which are the same as the above-used. Taking into consideration the first inequality from (4.1) in Lemma 4, we have to conclude that the switching function satisfies
Here are zeros of the switching function .
Please note that, by (3.15), relationships (5.35)–(5.37) also describe the behavior of the switching functions and .
Finally, we estimate the number of zeros of switching function . Taking into consideration the boundary condition (see (3.16)), we show that function has at most four distinct zeros on interval . Let us prove this by contradiction: assume that the function has at least five distinct zeros on the interval . Applying the Generalized Rolle’s Theorem, we conclude that in this case the function has at least four distinct zeros on the interval . This conclusion contradicts to the above-established estimate, and, therefore, our assumption is incorrect.
Denoting the zeros of function by , , and ; , and taking into consideration the second inequality from (4.1) in Lemma 4, we get the following relationship describing behavior of the switching function :
6. Types of the Optimal Controls
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 , , and . Specifically, depending on the sign of , the optimal controls and are of one of the following types:(i)If , then(ii)If , then
Here, and are the switchings of the optimal controls and , respectively.
The optimal control is of the following type:
where are the moments of switching.
The optimal control is described as following:
where are the moments of switching.
Taking into consideration (3.15), we can immediately describe the type of controls and as well. Specifically, depending on the sign of the difference , the optimal control belongs to one of the following two types:(i)If , then(ii)If , then
The optimal control is of the following qualitative type:
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 with nodes over time interval that is equivalent to the time step of . Since the problem is solved by a direct method using an iterative approach, at each step we impose the acceptable convergence tolerance . 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:
These values are adopted from [68]. We note that .
Results of the numerical calculations for are presented in Figures 1–4. Figures 1 and 2 correspond to the outcomes for the following bilinear incidence rate
with . Figure 1 show the optimal solutions , and , while Figure 2 display the corresponding optimal controls , , , , and . is the minimum value of the functional (see (2.15)). IPOPT used iterations in order to solve the optimal control problem (2.11) and (2.15).
Figures 3 and 4 present results for a nonlinear incidence rate of the Michaelis–Menten kinetics type with saturation in both variables:
where , , . Figure 3 show the optimal solutions , and , while Figure 4 display the corresponding optimal controls , , , , and . is the minimum value of the functional of (2.15). IPOPT used iterations in order to solve the optimal control problem (2.11) and (2.15).
Results of the numerical calculations for are presented in Figures 5–8. Figures 5 and 6 correspond to the outcomes for the bilinear incidence rate of (7.2). Figure 5 shows the optimal solutions , and , while Figure 5 displays the corresponding optimal controls , , , , and . is the minimum value of the functional (see (2.15)). IPOPT used iterations in order to solve the optimal control problem (2.11), (2.15).
Figures 7 and 8 present results for the nonlinear incidence rate of (7.3). Figure 7 show the optimal solutions , and , while Figure 8 display the corresponding optimal controls , , , , and . is the minimum value of the functional of (2.15). IPOPT used iterations in order to solve the optimal control problem (2.11) and (2.15).
The numerical calculations lead to the following conclusions:(i)A change in the type of the incidence rate does not lead to qualitative changes in the behavior of the optimal solutions , and .(ii)Taking into account the relationship (3.15) between the switching functions , and , , it can be seen that the optimal controls , and are piecewise constant functions with switchings. At the same time, the optimal controls , and are the constant functions over the whole interval . 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 function 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 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 conditions and the maximal values of derivatives and reached in the feasible region of the model. When the maximal number of switchings of the optimal controls is known, the two-point boundary value problem for Pontryagin Maximum Principle can be immediately reduced to a considerably simpler problem of finite-dimensional optimization. For specific parameter values and initial conditions, the latter problem can be solved numerically. However, in this paper, we did not use this approach; we used, instead, an alternative numerical solver package BOCOP–2.0.5 to verify and illustrate our analytical findings.
It follows from our analysis that for the considered model with a nonlinear incidence rate satisfying few biologically feasible assumptions, the principal qualitative properties of the optimal controls, such as the maximum number of switchings of the bang-bang controls and the order of these switchings, do not depend on a particular form of the nonlinearity. Additionally, our analysis indicates that possible nonlinearities of other functional responses in the model would not affect the principal properties of the optimal controls either, provided that the signs of the derivatives of these responses are preserved. We have to stress, however, that the conclusion that the principal properties of the optimal controls do not depend on a form of the nonlinearity of functional responses of a model holds only for the functionals which are independent of the controls and can be different for the explicitly depending on controls objective functionals.
While our numerical outcomes are presented for treatment duration of 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 controls (“immunization” of CD4+ T cells before they enter patient’s bloodstream) and (“immunization” of the existing CD4+ T cells) and corresponds to action of the fusion inhibitors. (Drugs that “immunize” CD4+ T cells before they enter patient’s blood stream do not exist at the moment, but such drugs can appear in the nearest future).(ii)The second group combines controls (inhibition of the virus production by the infected cells) and (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 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, (killing or immobilizing of the infected cells, or removing them by other means) and (inhibiting the infection rate). Control corresponds to action of NRTI, NtRTI, NNRTI and the integrase inhibitors. Drugs represented by control do not exist, but such drugs can be developed in the next few years.
It is remarkable that in the absence of control the drugs of the first, second and third groups are completely independent and do not interfere. This implies that in the absence of control , the controls of other three groups can be considered separately (which simplifies the analysis), whereas in the presence of this control all four groups should be considered as a complex. What is more important for real-life 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.
Appendices
A. Reduction of the Matrix of the Linear System
On time interval we consider the following third-order system of linear nonautonomous homogeneous differential equations:
where
Here are given Lebesgue measurable bounded functions; and are sign-definite almost everywhere on the interval functions; , are absolutely continuous solutions to system (A.1) satisfying the boundary conditions
where , are constants. Usually in this problems is the switching function and the object of investigation, and and are the corresponding auxiliary functions.
Our goal is to reduce the matrix to the upper-triangular form:
In our earlier papers [28, 34], we conducted such a transformation via a series of successive changes of variables. However, such a transformation can be achieved by a simultaneous nondegenerate transformation of a vector
or, in a matrix form, , where
Here, , and are functions should be chosen later. It is easy to see that such a transformation does not change function ; that is, holds. The other two functions satisfy
Moreover, the reverse transformation is also easily determined as
The corresponding inverse matrix in equality is
We apply transformation (A.5) to matrix . To do this, using variables , we rewrite system (A.1) as
Here we choose matrix , or that is the same, the functions , and , to satisfy the equality
We re-write equality (A.11) element-wise:
The matrix that satisfies these equalities has the required upper-triangular form (A.2). On its main diagonal the matrix has values
and, due to equalities (A.12) and (A.14), functions and satisfy the following differential equations:
To find a similar equation for function , we substitute (A.14) into equality (A.13). This yields the equation
Combining equations (A.16) and (A.17), we finally get a system of nonautonomous nonhomogeneous quadratic differential equations for functions , and :
The transformed system (A.10) in this case has the following form:
Initial conditions
where , and are constants, can be added to system (A.18).
As a result, we have the Cauchy problem, and, as a consequence of the Existence and Uniqueness of Solution Theorem [64], solutions , , to this problem can be defined only locally, in a neighborhood of the value . However, for this problem it is important to have solutions of the Cauchy problem (A.18), (A.20) defined on the entire interval (see [28, 34]). That is we have to find the restrictions on functions such that under these restrictions solutions of problem (A.18) and (A.20) were defined on the entire interval . Then system (A.19) would be also defined on the entire interval and, hence, it would be possible to use the Generalized Rolle’s Theorem [53] to estimate the number of zeros of switching function .
The analysis of the 3-dimensional quadratic system (A.18) can be significantly simplified by splitting this system into two simpler -dimensional quadratic subsystems. Specifically, it is easy to see that equations for and contain no variable and, hence, can be considered separately. These two equations form the first of these quadratic subsystems:
To formulate the second -dimensional subsystem, we introduce variable
Substituting (A.22) into the fist equation of system (A.18), we get equation
Besides, function satisfies the following differential equation:
Initial condition for variable immediately follows from (A.22):
Combining (A.23)–(A.25) and adding the initial condition for function from (A.20), we get the second -dimensional quadratic subsystem
We would like to emphasize that the -dimensional quadratic system (A.18) is equivalent to two -dimensional quadratic subsystems (A.21) and (A.26). As analysis of -dimensional systems is usually simpler, further we will deal with the later. In particularly, we have to analyze these quadratic subsystems to find the restrictions on functions which guarantee that there exist solutions of subsystems (A.21) and (A.26) defined on the entire interval .
B. Verifying Assumption 9
To verify that Assumption 9 holds for the optimal controls that we found in this paper and the corresponding optimal solutions, we consider system (2.11) for the case of bilinear incidence rate , where is a positive constant:
For this case Assumption 9 is
Our goal is to find the restrictions on the parameters and the initial conditions of system (2.11), as well as on values , , , , , , , , , , and defined in (2.5), which ensure that inequalities (B.2) are held.
Taking into consideration that all the controls are bang-bang and taking either the minimal, or the maximal possible value, we immediately see that the first two inequalities in (B.2) hold if the smallest values of their left-hand sides are greater than or equal to the largest values of the right-hand sides. That is, we have the inequalities:
Here and are the positive lower bounds for solutions and , respectively. In inequalities (2.12) of Lemma 3 for the general incidence rate , these values are zero. However, for the specific function the lower bounds and can be refined. Specifically, integrating the first equation in (B.1), we get the inequality
where
To find the lower bound , we, firstly, have to find the positive lower bound for by integrating the second equation of system (B.1). This yields
where . Then, using this inequality, we integrate the third equation of system (B.1) and get the inequality
where
We would like to note that the values and in inequalities (B.4) and (B.7) can be also found by the integration of the first and third equations of system (B.1) and a subsequent evaluation of the results. These values are found in (2.12), Lemma 3:
The value is also defined in (2.12).
The third inequality in (B.2) is satisfied if the smallest possible values of its left-hand side is greater than or equal to zero. Hence, we have the following inequalities:
Combining inequalities (B.3) and (B.10), we obtain the following system:
If these inequalities hold, then inequalities (B.2) are satisfied as well.
We verified the validity of inequalities (B.11) by direct calculations for the following parameters and initial conditions of system (2.11) and specific values , , , , , , , , and from (2.5):
Other specific values of the parameters and initial conditions of system (2.11) can be found in [69, 70].
Data Availability
The software used to support the findings of this study is free open-access software that can be freely downloaded from https://www.bocop.org/. The values of parameters and initial conditions used in computations to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
E. Khailov is funded by the Russian Foundation for Basic Research according to the research project 18-51-45003 IND_a.