In pursuant of some vital models for HIV dynamics and treatment progression, we identified and formulated as penultimate model, a set of 7–Dimension classical mathematical model, which accounted for the dynamical interplay of dual HIV – parasitoid pathogen infections on dual immune systems, studied using multiple chemotherapy cocktail in the presence of enhanced immune effectors response. The model was considered as a continuous multiple chemotherapy treatment (MCT) and as periodic multiple chemotherapy treatment (PMC), transformed to an optimal control problem. The positivity of the model state variables and stability properties was conducted. Deploying classical optimal control theory, the model used Pontryagin’s maximum principle to investigate the existence of optimal control strategy, established the optimality control system and justified the uniqueness of the system solutions. Numerical methods were explored to numerically solve the existing model via Runge–Kutter – 4 in a Mathcad surface. The result of the numerical analysis did not only identified PMC treatment as possible technique for the reduction of drug side–effects and suppression of dual HIV –pathogen infection by enhanced immune effectors response but largely established continuous MCT, which indicated complete elimination of dual HIV – pathogen viruses and provided window for quantification of minimized systemic cost as a more formidable approach in tackling the menace of the of the deadly dual infectivity. Thus, a broader verification and application of the model to related infectious disease is therefore suggested.

Keywords: multiple–chemotherapy–treatment, periodic–multiple–chemotherapy, dual–hiv–pathogen–infection, immune–effectors–response, optimal–weight–factor, optimal–control–measures, penalty–condition

MCT: multiple chemotherapy treatment; PMC: periodic multiple chemotherapy treatment; HIV: human immune deficiency virus; AIDS: acquired immune deficiency syndrome; STI: structured treatment interruptions; ODEs: ordinary differential equations;

In spite of lack of complete eradication of the world most acclaimed viral disease – the human immune deficiency virus (HIV), which often snowball into the deadly – acquired immune deficiency syndrome (AIDS), a lot has been tinkered by research scientists in the area of suppression and elimination of these infectious diseases and its affiliated pathogenic infections. Furthermore, following the multiplicity of HIV infections, with closest limitations, which inevitably assume the motivating factors of this present study, is the methodological application of multiple chemotherapies. A factor that comes along with the impairment of drugs resistance on continuous prolong treatment. Not left out are the consequences of optimal cost benefit from these chemotherapies as well as the determinative factors of CD8 immune effector cells and other immune responders considered as key players in the established viral load and pathogenic infection set–points.^1,2 From the study on dynamic multidrug therapies for HIV with optimal and structured treatment interruptions (STI) control approaches, the model.¹ described the interaction of the immune system with HIV and permitted drug “cocktail” therapies in the presence of structured treatment interruptions. The result showed how STI therapies can lead to long control of HIV by the immune response system after discontinuation of therapy.

Moreso, considering prevalence and strength of CD8 effector cells and other immune response, which affects disease progression rate, the models.^3,4 formulated treatment strategies that aimed at boosting adaptive cellular immune responses. Other models that accounted for structured treatment interruption strategies for single viral load infection treatment as well as the use of STI strategies includes.^5–10 Interesting, all the above models focuses on STI for single viral load infections, while the models.^2,11 had considered multiple treatment of dual HIV infectivity. Model.² studied quantitative approach to parametric identifiability of dual HIV – parasitoid infection; while the model.¹¹ investigated the analysis of parametric HIV infectivity and optimal control for the treatment of dual HIV–parasitoid pathogen infection. The results of both investigations were of immense contribution to the evaluation and suppression of viral load and pathogenic infections.

Therefore, in deploying ordinary differential equations (ODEs), the present study is formulated as a 7–Dimensional dual HIV–pathogen dynamic differential model. We present the model as an optimal control problem and the method of analysis explores classical numerical method – the Pontryagin maximum principle. The method leads to the establishment of existence and uniqueness of optimality control solution. The methodological application of the study involves clinical chemotherapy cocktail from the class of antiretroviral therapies known as reverse transcriptase inhibitors (RTIs) and protease inhibitors (PIs), collectively called highly active antiretroviral therapy (HAART). Inclusively, the combination of RTI and PIs drugs increases CD4⁺ T cell count, which are the target cells of viral load and pathogeneses.¹² Thus, the present model, which utilizes drug cocktail from RTI and PIs, is formulated as a dynamic intermittent (periodic) treatment and methodological application of periodic multiple chemotherapies (PMC) for the prevention and suppression of dual viral load and pathogenic infections on host vector – human immune system. Furthermore, the model is task with the maximization of healthy CD4⁺ T cells and minimization of optimal cost, a procedure that adopts the application of optimal control technique (approach). For more details on the above, we refer readers to the models.^13–17 Therefore, the present paper as against those mentioned above, is aimed at resolving the pertinent issues in the range of methodological application of periodic multiple chemotherapies, complete zero application of optimal control measures on multiple drug cocktail, periodic alternate zero application of multiple drug cocktail, maximization of healthy CD4⁺ T cell count and evaluation of the strength of immune responses; and as well, the minimization of systemic cost. Two more important factors apprehensive of this model are drug validity period, which often span within 500–750 days before the onset of drug side–effect and drug resistivity.^{7,9,13,15,17–19} and initiation of treatment time from onset of infection, the detail which can be found in.²

Thus, the structural content of this work is characterized by eight sections, with section 1, covering the introductory aspect. Devoted in section 2, is the material and method used for the formulation and schematic presentation of the model as a classical 7–Dimensional mathematical differential model. The corresponding optimal control strategy and optimality system characterized by optimal control solution forms the fulcrum of section 3. Section 4 is devoted to the formulation and derivation of optimal multiple periodic chemotherapy (MPC) theoretical approach for the control of viral load and parasitoid–pathogen. In section 5 we conducted numerical computations for the continuous optimal multiple chemotherapy without control measures on treatment factors. The model is stepped further with the numerical simulations for the continuous multiple chemotherapy treatment with control measures on treatment optimal weight factors in section 6. Constituting section 7 is an explicit discussion of the results of conducted experiment. Finally, incisive model conclusion and recommendations are given in section 8. It is anticipated that this present study will offer a brighter hope to the eradication of the deadly disease.

The material and methods of this study is subdivided into two sections (2 and 3) respectively. Embedded in section 2, are the schematic representation of system model followed by derivation of model optimal control problem. We also show here that the state variables are non–negative and as well, ascertain the stability behavior of the system. Section 3, is devoted to the formulation of optimal control strategy and the optimality system of the model, which explores classical numerical methods known as Pontryagin maximum principle.

Schematic representation and model formulation

In structuring the progression of this model, we condone the fact that there are two virions (viral load and pathogen), which attack the immune system composing of two target cells (T–lymphocytes and macrophages). The scheme adopts the introduction of two chemotherapy cocktail (RTI and PIs) and taken into cognizance, natural interplay of adaptive immune effectors response. Thus, the model is constituted by 7 – subgroups, schematically represented as in Figure 1 below:

Figure 1 Schematic representation of dual HIV-pathogen infection with PMC treatment.

Obviously, from Figure 1, if these subgroups represent the population variables, with unit volume $m{m}^3$ , then we define as follows: $U_i$  i=1, 2 as uninfected T–lymphocytes and macrophages cells, $U_i^{\ast }$ , i=1, 2 as infected T–lymphocytes and macrophages cells, V – viral load, P – pathogen and M – immune effectors. Therefore, using classical ODEs, the simulative interactions of the above biological structure, is the formulation of a 7–Dimensional dynamic optimal control problem having its physiological derivation as follows:

$U_1{}^{\prime }=b_1-{\alpha }_1{U}_1-(1-{\rho }_1)g_1(V+P)U_1$   (1)

$U_2{}^{\prime }=b_2-{\alpha }_2{U}_2-(1-\frac{{\rho }_1}{r_1})g_2(V+P)U_2$   (2)

$U_1^{\ast }{}^{\prime }=(1-{\rho }_1)g_1(V+P)U_1-\mu U_1^{\ast }-h_{1}M{U}_1^{\ast }$   (3)

$U_2^{\ast }{}^{\prime }=(1-\frac{{\rho }_1}{r_1})g_2(V+P)U_2-\mu U_2^{\ast }-h_{2}M{U}_2^{\ast }$   (4)

$V^{\prime }=(1-{\rho }_2)p_{p}V-(\mu +n)V-[(1-{\rho }_1){\gamma }_1{g}_1{U}_1+(1-\frac{{\rho }_1}{r_1}){\gamma }_2{g}_2{U}_2]V$   (5)

$P^{\prime }=(1-\frac{{\rho }_2}{r_2})p_{p}P-(\mu +n)P-[(1-{\rho }_1){\gamma }_1{g}_1{U}_1+(1-\frac{{\rho }_1}{r_1}){\gamma }_2{g}_2{U}_2]P$   (6)

$M^{\prime }=b_M+\frac{w_M(U_1^{\ast }+U_2^{\ast })}{(U_1^{\ast }+U_2^{\ast })+H_w}M-\frac{q_M(U_1^{\ast }+U_2^{\ast })}{(U_1^{\ast }+U_2^{\ast })+H_q}M-{\mu }_{M}M$   (7)

with initial values $U_1(0)=U_{(1)}{}_0,U_2(0)=U_{(2)}{}_0,U_1^{\ast }(0)=U_{{}_{(1)}}^{\ast }{}_0,U_2^{\ast }(0)=U_{{}_{(2)}}^{\ast }{}_0,V(0)=V_0,P(0)=P_0$  and $M(0)=M_0$  for all $t=t_0$  and which satisfies the biological variables and parameter values as defined in Tables 1 & 2 below. The present model is in tune with the investigation by.²⁰ where the author accounted for only single drug–RTI on single infection – HIV; and the study.¹ accounted for chemotherapy “cocktail” (RTI and PIs) used for the study of dynamics of single infection–HIV.

Laconically, equations (1)–(7) present the model of the system and the epidemiological terms are carefully defined thus: in equation (1),$b_1$  – the source term of uninfected CD4⁺ T–lymphocytes and ${\alpha }_1{U}_1$  – the death rate of uninfected T–lymphocytes cells. The term $(1-{\rho }_1)g_1(V+P)U_1$  is the interplay of viruses $(V,P)$  with infectivity rate $g_1$ , and ${\rho }_1(t)$  – drug efficacy, which model the RTI that block new infections, in $U_1,U_2,U_1^{*}$  and $U_2^{\ast }$  of the population. From equation (2),$b_2$  – source term of uninfected macrophages cells, ${\alpha }_2{U}_2$  – the death rate of uninfected macrophages cells. Here, the term $(1-{\rho }_1/r_1)g_2(V+P)U_2$  – represents the balance of treatment of viruses with infectivity rate $g_2$ , on $U_2$  and ${\rho }_1/r_1$  – as treatment efficacy reduction in population $U_2$ . The implication is that $g_1,g_2$  and the ratio $1/r_1$  constitutes the difference between the two cell populations. For these two cells, we assume that $0\le x_1\le {\rho }_1(t)\le y_1<1$ , such that $x_1$ and $y_1$ denotes minimal and maximal drug efficacy respectively with $r_1\in [0,1]$ .

 In equations (3) and (4), the first terms $(1-{\rho }_1)g_1(V+P)U_1$  and $(1-{\rho }_1/r_1)g_2(V+P)U_2$ are the results of the interactions of free CD4⁺ T–lymphocytes and macrophages with free virions $(V,P)$  infectivity and the activation of proportionate (RTI) drug on infected $U_1^{\ast }$  and $U_2^{\ast }$ . The terms on infected $\mu U_1^{\ast }$  and $\mu U_2^{\ast }$ are the death rates of infected $U_1^{\ast }$  and $U_2^{\ast }$ . At the onset of infection, the third terms $h_{i}M{U}_i^{\ast }$ , i=1, 2, represents the cytotoxic T–lymphatic (CTL) also known as CD8⁺ T cells, which exhibit the tendency of detecting and killing infected cells as much as possible. Therefore, in equations (3) and (4), terms $h_{i}M{U}_i^{\ast }$ , i=1, 2, defines death of infected cells at the rate $h_{i}M$ , which depends on the density of immune effectors M.

$J=\left(\begin{array}{ccccccc}-{\alpha }_1-g_1(V+P)& 0& 0& 0& -g_1{U}_1& -g_1{U}_1& 0\\ 0& -{\alpha }_2-g_2(V+P)& 0& 0& -g_2{U}_2& -g_2{U}_2& 0\\ g_1(V+P)& 0& -\mu -h_{1}M& 0& g_1{U}_1& g_1{U}_1& -h_1{U}_1^{\ast }\\ 0& g_2(V+P)& 0& -\mu -h_{2}M& g_2{U}_2& g_2{U}_2& -h_2{U}_2^{\ast }\\ -{\gamma }_1{g}_{1}V& -{\gamma }_2{g}_{2}V& 0& 0& \begin{array}{l}{p}_v-(\mu +n)-({\gamma }_1{g}_1{U}_1\\ +{\gamma }_2{g}_2{U}_2)\end{array}& 0& 0\\ -{\gamma }_1{g}_{1}P& -{\gamma }_2{g}_{2}P& 0& 0& 0& \begin{array}{l}{p}_p-(\mu +n)-({\gamma }_1{g}_1{U}_1\\ +{\gamma }_2{g}_2{U}_2)\end{array}& 0\\ 0& 0& B_{7,3}& B_{7,4}& 0& 0& B_{7,7}\end{array}\right)$   (10)

From equation (5), the first term $(1-{\rho }_2)p_{v}V$ defines the control function with ${\rho }_2(t)$ , representing the efficacy of protease inhibitors and having viral load multiplicative capacity $p_v$ . The second term $(\mu +n)V$ defines reduction of infected cells due to clearance rate and the natural death rate of viral load. Thus, the productivity $p_v$  is reduced to $(1-{\rho }_2)p_v$ , where $0\le x_2\le {\rho }_2(t)\le y_2<1$ . The cumulative effect of viral load on both targeted cells in $U_1$  and $U_2$  is represented by the term $-[(1-{\rho }_1){\gamma }_1{g}_1{U}_1+(1-{\rho }_1/r_1){\gamma }_2{g}_2{U}_2]V$ , which reveal ${\gamma }_i$ , i=1, 2 as probability of multiple virions infections on target cells $U_1$  and $U_2$ . In a chronological manner, equation (6) can be interpreted as in equation (5) but with respect to pathogenic infection and having differential effect of $1/r_2$ as treatment reduction efficacy on $U_1^{\ast }$  and $U_2^{\ast }$ .

Remarkably, the clearance rate of infected cells.¹ is conspicuously visible through the immune effector cells (CTL), designated by M, of equation (7). The dynamics of $M^{\prime }$ is adopted from the study of.²¹ The term $b_M$  – denotes the replication (or proliferation) of more effector cells from the co–existence of infected cell and immune effector cells. The second term $w_M(U_1^{\ast }+U_2^{\ast })M/(U_1^{\ast }+U_2^{\ast })+H_w$ – defines the maximum birth rate of immune effectors in the presence of infected cells$U_1^{\ast }$ , $U_2^{\ast }$  and $M$ ; and having saturation birth constant $H_w$ . The third term $q_M(U_1^{\ast }+U_2^{\ast })M/(U_1^{\ast }+U_2^{\ast })+H_q$  is attributed to immune effector maximum death rate, where $H_q$  is the effector death saturation constant. Lastly, the term ${\mu }_{M}M$ is the immune effectors natural death rate. Thus, the aforementioned descriptions are relevant in the administration of chemotherapy cocktail and PMC scenarios. Furthermore, the inclusion of immune effectors is a direct consequence of the sensitive role in the practical sense of PMC, which we will illustrate in the course of our numerical simulations.

From model (1)–(7) we observe that the function ${\rho }_i(t)$ , $i=1,2$ , determines the percentage drug efficacy and the maximal use of chemotherapies. Thus, ${\rho }_i(t)$ is a measurable function having limit $t\in [t_0,t_f]$  on the domain $0\le {\rho }_i(t)\le 1$ . This time interval accounts for the design of periodic treatment and as well, satisfies the objectives of this study as clearly mentioned in the introductory section. Furthermore, to ensure that our objective functional satisfies the control system, we take $t\in [t_0,t_f]$  such that $t\le 30$  months of drugs allowable validity period.^2,13,18

Therefore, the validity of the system equations (1)–(7) comes to bear, if we can genuinely generate compactible numerical values for the variables and parameters of the model. Attaining this, we adopt to the closest parameters values from those existing and valid studies as contain in our literature.^1,14,17 Thus, the simulative parameter values that satisfy the biological variables and parameters of our model are explicitly defined as in Tables 1 & 2 below:

Next, to satisfactorily establish the optimality and methodological application of our chemotherapy cocktail, it become necessary to first show that model state variables are non–negative and as well, discuss the disease stability analysis.

Positivity and compatibility of model variables

It is obvious that since maximal cost of chemotherapy cocktail is given by

${({\rho }_1(t)+{\rho }_2(t))}^2$ , then possible drug severities is accounted for by the following proposition.

 Proposition 2.1: If drug severities (or drug hazardous side–effects) emerges in the course of treatment, then the introduction of optimal weight factors ${\left\{K_i\ge 0,L_i\ge 0,\delta \ge 0\right\}}_{i=1,2}$ is justified. Also, if ${\left\{x_i,y_i\right\}}_{i=1,2}$  represents minimal and maximal drugs efficacies, then the inequality $0\le x_i\le \left\{K_i(t),L_i(t),\delta (t)\right\}\le y_i<1$ , i =1, 2 holds. Therefore, positivity for which model variables are compactible is studied using the following theorem.

 Theorem 2.1: If equations (1)–(7) represents the model equations, such that for $\varphi =\left\{\begin{array}{l}{U}_1,U_2,U_1^{\ast },U_2^{\ast },V,P,M\in {ℝ}_{+}^7:U_1(0)>0,U_2(0)>0,U_1^{\ast }(0)>0,U_2^{\ast }(0)>0,\\ V(0)>0,P(0)>0,M(0)>0\end{array}\right\}$ , then the solution of the system $\left\{U_1(t),U_2(t),U_1^{\ast }(t),U_2^{\ast }(t),V(t),P(t),M(t)\right\}$ are non–negative and compactible for all $t\ge 0$ .^22–24

Proof: We show that the state variables are life varying integers. Thus, invoking equations (1)–(7) step–wisely, we differentiating each of each them to justify their non–negativity.

From equation (1), we have $U_1{}^{\prime }=b_1-{\alpha }_1{U}_1-(1-{\rho }_1)g_1(V+P)U_1$ . Differentiating with respect to $U_1$ , we obtain $\frac{d{U}_1}{dt}\ge -{\alpha }_1{U}_1$  $\Rightarrow$   $\frac{d{U}_1}{dt}+{\alpha }_1{U}_1\ge 0$ . Since $t\ge 0$ , then $U_1=+\infty$  and is positive. The integrating factor IF is given by $IF=e^{{\displaystyle \int {\alpha }_{1}dt}}=e^{{\alpha }_{1}t}$ . Multiplying the equation by the integrating factor, we have $$e^{{\alpha }_{1}t}\{U_1{}^{\prime }+({\alpha }_1)U_1\}\ge 0$$ . Rewriting the left hand side of the equation, we obtain $$\frac{d}{dt}(e^{{\alpha }_{1}t}{U}_1)\ge 0$$ .

Integrating both sides, leads to $\frac{d}{dt}(e^{{\alpha }_{1}t}{U}_1)\ge 0$ . Dividing by the integrating factor, we have $$U_1(t)=C{e}^{{\alpha }_{1}t}$$ . Applying the initial condition i.e. $t=0$ , $U_1(t)=U_1(0).U_1(0)\ge C$ , then $U_1(t)\ge U_1(0)e^{{\alpha }_{1}t}\ge 0$ ,$t\ge 0$ . Therefore, $U_1(t)\ge 0$  is non–negative and compactible. Next, taking equation (2), we set

$\frac{d{U}_2}{dt}=b_2-{\alpha }_2{U}_2-(1-\frac{{\rho }_1}{r_1})g_2(V+P)U_2$

Differentiating, we have, $\frac{d{U}_2}{dt}\ge -{\alpha }_2{U}_2$  and taking the integral gives ${\displaystyle \int \frac{d{U}_2}{U_2}}\ge -{\displaystyle \int {\alpha }_{2}dt}$ . Then applying the integrating factor $IF=e^{-{\alpha }_{2}t}$ , we have $U_2(t)\ge U_2(0)e^{-{\alpha }_{2}t}\ge 0$ , $t\ge 0$ . Applying initial condition at $t=0$ , $U_2(t)\ge U_2(0)$ . Therefore, for $t\to \infty$ , $U_2(t)\ge 0$ . Similarly, from equation (3), we have,

$\frac{d{U}_1^{\ast }}{dt}=(1-{\rho }_1)g_1(V+P)U_1-\mu U_1^{\ast }-h_{1}M{U}_1^{\ast }$ .

Differentiating, we see that $\frac{d{U}_1^{\ast }}{dt}\ge -\mu U_1^{\ast }$  $\Rightarrow$  $\frac{d{U}_1^{\ast }}{U_1^{\ast }}\ge -\mu dt$ and taking the integral gives ${\displaystyle \int \frac{d{U}_1^{\ast }}{U_1^{\ast }}}\ge -{\displaystyle \int \mu dt}$ . Applying the integrating factor $IF=e^{-\mu t}$ , we have $U_1^{\ast }(t)\ge U_1^{\ast }(0)e^{-\mu t}\ge 0$ , $t\ge 0$ . Hence, $t\to \infty$ , $U_1^{\ast }(t)\ge 0$ .

From equation (4), we have, $U_2^{\ast }{}^{\prime }=(1-\frac{{\rho }_1}{r_1})g_2(V+P)U_2-\mu U_2^{\ast }-h_{2}M{U}_2^{\ast }$ . Differentiating with respect to $U_2^{\ast }$ , we obtain, $\frac{d{U}_2^{\ast }}{dt}\ge -\mu U_2^{\ast }$  $\Rightarrow$  $\frac{d{U}_2^{\ast }}{U_2^{\ast }}\ge -\mu dt$  and taking the integral gives ${\displaystyle \int \frac{d{U}_2^{\ast }}{U_2^{\ast }}}\ge -{\displaystyle \int \mu dt}$ . Applying the integrating factor $IF=e^{-\mu t}$ , we have $U_2^{\ast }(t)\ge U_2^{\ast }(0)e^{-\mu t}\ge 0$ , $t\ge 0$ . When $t\to \infty$ , . Furthermore, taking equation (5), we set $V^{\prime }=(1-{\rho }_2)p_{p}V-(\mu +n)V-[(1-{\rho }_1){\gamma }_1{g}_1{U}_1+(1-\frac{{\rho }_1}{r_1}){\gamma }_2{g}_2{U}_2]V.$

Differentiating with respect to $V$ , we obtain, $\frac{dV}{dt}\ge [(1-{\rho }_2)p_v-(\mu +n)]V$ $\Rightarrow$ $\frac{dV}{V}\ge [(1-{\rho }_2)p_v-(\mu +n)]dt$ .

Taking the integral gives the expression ${\displaystyle \int \frac{dV}{V}}\ge {\displaystyle \int [(1-{\rho }_2)p_v-(\mu +n)]dt}$ . Now, applying the integrating factor $IF=e^{[(1-{\rho }_2)p_v-(\mu +n)]t}$ , we have $V(t)\ge V(0)e^{[(1-{\rho }_2)p_v-(\mu +n)]t}\ge 0$ , $t\ge 0$ . Therefore, when $t\to \infty$ , $V(t)\ge 0$ .

Taking equation (6), we have $P^{\prime }=(1-\frac{{\rho }_2}{r_2})p_{p}P-(\mu +n)P-[(1-{\rho }_1){\gamma }_1{g}_1{U}_1+(1-\frac{{\rho }_1}{r_1}){\gamma }_2{g}_2{U}_2]P$ . Differentiating with respect to P, we obtain, $\frac{dP}{dt}\ge [(1-\frac{{\rho }_2}{r_2})p_p-(\mu +n)]P$  $\Rightarrow$  $\frac{dP}{P}\ge [(1-\frac{{\rho }_2}{r_2})p_p-(\mu +n)]dt.$  Taking the integral gives the expression ${\displaystyle \int \frac{dP}{P}}\ge {\displaystyle \int [(1-\frac{{\rho }_2}{r_2})p_p-(\mu +n)]dt}$ . Now, applying the integrating factor $IF=e^{[(1-\frac{{\rho }_2}{r_2})p_p-(\mu +n)]t}$ , we have $P(t)\ge P(0)e^{[(1-\frac{{\rho }_2}{r_2})p_p-(\mu +n)]t}\ge 0$ , $t\ge 0$ . Therefore, when $t\to \infty$ , $P(t)\ge 0$ .

Finally, from equation (7), we have, $M^{\prime }=b_M+\frac{w_M(U_1^{\ast }+U_2^{\ast })}{(U_1^{\ast }+U_2^{\ast })+H_w}M-\frac{q_M(U_1^{\ast }+U_2^{\ast })}{(U_1^{\ast }+U_2^{\ast })+H_q}M-{\mu }_{M}M$ .

Differentiating with respect to M, we obtain, $\frac{dM}{dt}\ge -{\mu }_{M}M$  $\Rightarrow$  $\frac{dM}{M}\ge -{\mu }_{M}dt.$  Taking the integral gives the expression ${\displaystyle \int \frac{dM}{M}}\ge {\displaystyle \int -{\mu }_{M}dt}$ . Now, applying the integrating factor $IF=e^{-{\mu }_{M}t}$ , we have $M(t)\ge M(0)e^{-{\mu }_{M}t}\ge 0$ , $t\ge 0$ . Therefore, when $t\to \infty$ , . Hence, all the model variables are non–negative and this completes the proof. We complete this section by highlighting the properties for which the PMC model is locally asymptotically stable.

Model stability analysis

Though the system stability analysis seems somewhat complex due to the nature of nonlinear systems involved, we’ll discuss the existing stability properties of the model in a simpler manner. Suppose, $\nu =(U_1,U_2,U_1^{\ast },U_2^{\ast },V,P,M)$  defines the vectorial capacity of the model, then the system (1)–(7) satisfies the equation

$\frac{d\nu (t)}{dt}=f(t,\nu ;z)$   (8)

Where $f(t,\nu ;z)$  is the right–side of the ODE system and , the vector parameters as in Table 2. Then, we implore the Runge–Kutter of order 4, to solve the equation $f(t,\nu ;z)=0$ for the steady states ${\overline{\nu }}_k$ . Next, we calculate the Jacobian matrix of the partial derivatives of the right side of the differential equations with respect to the state variables, i.e. $\frac{\partial f(t,\nu ;z)}{\partial \nu }=\left[\frac{\partial f_i(t,\nu ;z)}{\partial {\nu }_j}\right]$   (9).

From equation (1)–(7), if we let ${\rho }_1={\rho }_2=0$ for the simple reason that we are poised to establish the system stability behavior when off treatment is applied, then the Jacobian matrix is derive as:

$J=\left(\begin{array}{ccccccc}-{\alpha }_1-g_1(V+P)& 0& 0& 0& -g_1{U}_1& -g_1{U}_1& 0\\ 0& -{\alpha }_2-g_2(V+P)& 0& 0& -g_2{U}_2& -g_2{U}_2& 0\\ g_1(V+P)& 0& -\mu -h_{1}M& 0& g_1{U}_1& g_1{U}_1& -h_1{U}_1^{\ast }\\ 0& g_2(V+P)& 0& -\mu -h_{2}M& g_2{U}_2& g_2{U}_2& -h_2{U}_2^{\ast }\\ -{\gamma }_1{g}_{1}V& -{\gamma }_2{g}_{2}V& 0& 0& \begin{array}{l}{p}_v-(\mu +n)-({\gamma }_1{g}_1{U}_1\\ +{\gamma }_2{g}_2{U}_2)\end{array}& 0& 0\\ -{\gamma }_1{g}_{1}P& -{\gamma }_2{g}_{2}P& 0& 0& 0& \begin{array}{l}{p}_p-(\mu +n)-({\gamma }_1{g}_1{U}_1\\ +{\gamma }_2{g}_2{U}_2)\end{array}& 0\\ 0& 0& B_{7,3}& B_{7,4}& 0& 0& B_{7,7}\end{array}\right)$ (10)

Where

$B_{7,3}=B_{7,4}=\frac{w_M{H}_{w}M}{{(U_1^{\ast }+U_2^{\ast }+H_w)}^2}-\frac{q_M{H}_{q}M}{{(U_1^{\ast }+U_2^{\ast }+H_q)}^2}$

and

$B_{7,7}=\left(\frac{w_M}{U_1^{\ast }+U_2^{\ast }+H_w}-\frac{q_M}{U_1^{\ast }+U_2^{\ast }+H_q}\right)U_1^{\ast }+U_2^{\ast }-{\mu }_M$ .

So equation (10) exhibits non–singularity behavior since the diagonal of the Jacobian matrix is non–zero.

Then if we substitute the resulting computation of ${\overline{\nu }}_k$  steady state for v, in equation (10), we obtain the ODE system dynamics that is linearized about the equilibrium ${\overline{\nu }}_k$ . So we see from here that linearized ODE theory ascertain the fact that if the eigenvalues of the matrix have all negative real parts, then the equilibrium ${\overline{\nu }}_k$ is locally asymptotically stable. Therefore, given specific parameter values as in Tables 1 & 2, the model (1)–(7) exhibits three physical steady states and several non–physical steady states (omitted here for brevity). Thus, it is of interest to note that detail analysis of these stability behaviors is left to readers with related analysis as in model.^1,8,25 This thought is in line with the focus of the present study – optimization control of PMC treatment and the maximization of CD4⁺ T cells and macrophages.

Unlike model.¹ which considered single infection (HIV) in two immune systems with dual chemotherapy cocktail under STI program and which the method of analysis explored stability and linearization technique, the present model consider dual infectious variables under PMC program and focuses on utilizing classical numerical method known as Pontryagin’s maximum principle. This method allows the verification of existence of model as a function and the uniqueness of the system solution. To accomplish this task, we establish the model optimality control strategy from the derived optimal control problem.

Optimal control strategy and optimality system

Having shown that the model state variables are non–negative with known stability behavioral pattern, we then establish in this section, the optimal control strategy for continuous chemotherapy cocktail and the consequences following the introduction of optimal weight factors. This will lead to the derivation of existence of the model, establish the model optimality control system and lastly, prove the uniqueness of the solution of the system.

Optimal control strategy for continuous chemotherapy cocktail

We recall that optimal control strategy is the derivation of mathematical model (as in this case of model (1)–(7)) with which we define the objective functional as a function of maximization of the control variables. Indeed the objective functional of any optimal control problem is an integral equation, which model the trade–off between virions and pathogen concentration, organ health and use of therapeutics.^26,27 Therefore, for a HIV–pathogen dynamics having optimal control problem as in equations (1)–(7), the objective functional that maximizes the control system is derive as:

$R({\rho }_1,{\rho }_2)={\displaystyle \underset{t_0}{\overset{t_f}{\int }}[K_{1}V(t)+K_{2}P(t)+L_1{\rho }_1^2+L_2{\rho }_2^2}-\delta M(t)]dt$   (12)

where ${\rho }_1(t),{\rho }_2(t)$ are the control variables for RTIs and PIs respectively. The quantities $K_{i=1,2},L_{i=1,2}$  and δ are the optimal weight factors on the virions, control treatment inputs and immune effectors respectively. These control constants maximizes the system benefits quantified along the level of CD4⁺ T cells concentration as represented by the first and second terms of equation (12). The third and fourth terms of the equations accounts for the minimization of systemic cost of drugs treatment cocktail. So that, if $0\le {\rho }_{i=1,2}\le 1$  represents maxima drug usage, then periodic zero alternate application of multiple drug cocktail is visible. Therefore, the maximal cost of drug usage is given by ${({\rho }_1(t),{\rho }_2(t))}^2$ .^28,29

The introduction of optimal weight factors ${\left\{K_i\ge 0,L_i\ge 0,\delta \ge 0\right\}}_{i=1,2}$  is a consequence of the fact that cost benefits are nonlinear. Thus, these serve as simple nonlinear controls on the system variables. And we can satisfactorily say that the objective functional (12) completely meets the aims of the present investigation, which is primed by maximization of immune systems and immune effector concentration in the presence of minimal application of systemic cost as well as maximal suppression of both viral load and pathogenic infections. Therefore, we compatibly opt for an optimal control pair $({\rho }_1^{\ast },{\rho }_2^{\ast })$  that is ascribes by the expression ${}_{0\le {\rho }_i\le 1}{}^{\max }R({\rho }_i^{\ast })=\min {\left\{R({\rho }_1,{\rho }_2)/({\rho }_1,{\rho }_2)\in Q\right\}}_{i=1,2}$  subject to the system of ODEs (1)–(7) and such that $Q=\{R({\rho }_1,{\rho }_2)/{\rho }_i,Q={\rho }_i/{\rho }_i$  is measurable with $x_i\le {\rho }_i\le y_i$  for all $t\in [t_0,t_f]$ , for $i=1,2\}$ a measurable control set. This is an innovative optimal control theory, which is in line with several varying control theories formulated; see for example.^7,9,13,17,19 Next, we verify the existence of the model optimality control for a PMC treatment.

 Existence of optimality control strategy

The existence of an optimality control for PMC treatment of a dual HIV–pathogen infections can be proved from the point of.¹⁰ where we possibly show that the right sides of equations (1)–(7) are bounded by a linear function of the state and control variables. Also, we show that the integrand of the objective functional (12) is concave on Q and bounded below, which again affirm the compatibility required of the model. Therefore, there exists on the basis of system boundedness of solution, super–solutions of the system ${\overline{U}}_1,{\overline{U}}_2,{\overline{U}}_1^{\ast },{\overline{U}}_2^{\ast },\overline{V},\overline{P},$ and , satisfying the equation:

$\{\begin{array}{l}\frac{d{\overline{U}}_1}{dt}=b_1-[{\rho }_1(t)]{\overline{U}}_1\\ \frac{d{\overline{U}}_2}{dt}=b_2-[{\rho }_1(t)]{\overline{U}}_2\\ \frac{d{\overline{U}}_1^{\ast }}{dt}=\frac{h_1{\overline{U}}_1^{\ast }}{C_1+{\overline{U}}_1^{\ast }}\\ \frac{d{\overline{U}}_2^{\ast }}{dt}=\frac{h_2{\overline{U}}_2^{\ast }}{C_1+{\overline{U}}_2^{\ast }}\\ \frac{d\overline{V}}{dt}=\frac{{\rho }_2(t)\overline{V}}{C_2+\overline{V}}\\ \frac{d\overline{P}}{dt}=\frac{{\rho }_2(t)\overline{P}}{C_2+\overline{P}}\\ \frac{d\overline{M}}{dt}=\frac{b_M-({\rho }_1(t)+{\rho }_2(t))\overline{M}}{C_3+\overline{M}}\end{array}$  ,(13)

Where $C_{i=1,2,3}$  are half saturation constants on $U_1^{*},U_2^{*};V$ and P with equation (13) been compactly bounded on a finite time interval. Thus, determining the existence of the optimal control to the model, we invoke theorem 4.1, p. 68–69.¹⁰

Theorem 3.1: Given an optimal control system with model equations (1)–(7) and having proposition 2.1, there exists a PMC optimal control ${\overrightarrow{\rho }}^{\ast }=({\rho }_1^{\ast },{\rho }_2^{\ast })\in Q$ such that

$\begin{array}{c}\max \\ ({\rho }_1,{\rho }_2)\in Q\end{array}R({\rho }_1,{\rho }_2)=R({\rho }_1^{\ast },{\rho }_2^{\ast })$

Proof: Taking proceeding from Thm. 4.1.¹⁰ we state and show that the following conditions are justified:

1. Then control set ${\left\{{\rho }_i\right\}}_{i=1,2}$ is Lebesgue–integrable function on interval $[t_0,t_f]$ and the corresponding state variable is nonempty.
2. The measurable control set Q, is convex and closed.
3. The right side (RHS) of the state system is continuous and bounded above by sum of the bounded control and state variables and can by written as a linear function of ${\rho }_i$ , i=1,2 , with coefficients define by proposition 2.1 and on the state variables.
4. The integrand of the objective functional is concave on the measurable set Q.
5. There exist consistence $b_1,b_2>0$ and $\beta >1$ , such that the integrand $L(U_1,U_2,{\rho }_1,{\rho }_2)$ of the objective functional satisfies

$L(U_1,U_2,{\rho }_1,{\rho }_2)\le b_2-b_1{({\left|{\rho }_1\right|}^2+{\left|{\rho }_2\right|}^2)}^{\beta /2}$  .

Then, we at once resort to the result of (.³⁰ Thm. 9.2.1, p.182) for the existence of the state system (1)–(7) with bounded coefficients, which satisfies condition (i). It is obvious that solution here is bounded. It follows by definition that our control set is closed and convex, thus satisfying condition (ii). Since, our state system is bilinear in ${\rho }_1,{\rho }_2$ , the RHS of (1)–(7) satisfies condition (iii) in the sense of boundedness of solutions. Furthermore, the integrand of the objective functional $\left(K_{1}V(t)+K_{2}P(t)+L_1{\rho }_1^2+L_2{\rho }_2^2-\delta M(t)\right)$ is concave on the measurable control set Q. Finally, the completeness of the existence of solution of the optimal control lies in the established fact that $\left[K_{1}V(t)+K_{2}P(t)+L_1{\rho }_1^2+L_2{\rho }_2^2-\delta M(t)\right]\le b_2-b_1({\left|{\rho }_1\right|}^2+{\left|{\rho }_2\right|}^2)$  where $b_2$  depends on the upper bound on V and P; and $b_1>0$  since ${\{L_i,K_i,\delta \}}_{i=1,2}>0$ . This completes the proof.

Optimality control system

Here, relying on the fact that we were able to prove the existence of the model equations, then we can establish the model optimality system. This is achieved by first defining the necessary conditions for an optimal control for periodic multiple treatments under dual infectious variables. Then, for equation (12), the penalty term on the constraints of the objective functional is the Hamiltonian arguments define by the Lagrangian:$L(U_1,U_2,U_1^{\ast },U_2^{\ast },V,P,M,{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4,{\tau }_5,{\tau }_6,{\tau }_7)$  =

$\begin{array}{l}{K}_{1}V+K_{2}P+L_1{\rho }_1^2+L_2{\rho }_2^2+\delta M\\ +{\tau }_1[b_1-{\alpha }_1{U}_1-(1-{\rho }_1)g_1(V+P)U_1]\\ +{\tau }_2[b_1-{\alpha }_2{U}_2-(1-\frac{{\rho }_1}{r_1})g_2(V+P)U_2]\\ +{\tau }_3[(1-{\rho }_1)g_1(V+P)U_1-\mu U_1^{\ast }-h_{1}M{U}_1^{\ast }]\\ +{\tau }_4[(1-\frac{{\rho }_1}{r_1})g_2(V+P)U_2-\mu U_2^{\ast }-h_{2}M{U}_2^{\ast }]\\ +{\tau }_5\{(1-{\rho }_1)p_{v}V-(\mu +n)V-[(1-{\rho }_1){\gamma }_1{g}_1{U}_1+(1-\frac{{\rho }_1}{r_1}){\gamma }_2{g}_2{U}_2]V\}\\ +{\tau }_6\{(1-\frac{{\rho }_2}{r_2})p_{P}P-(\mu +n)P-[(1-{\rho }_1){\gamma }_1{g}_1{U}_1+(1-\frac{{\rho }_1}{r_1}){\gamma }_2{g}_2{U}_2]P\}\\ +{\tau }_7\left[b_M+\frac{w_M(U_1^{\ast }+U_2^{\ast })}{(U_1^{\ast }+U_2^{\ast })+H_w}M-\frac{q_M(U_1^{\ast }+U_2^{\ast })}{(U_1^{\ast }+U_2^{\ast })+H_q}M-{\mu }_{M}M\right]\\ -q_{11}({\rho }_1-x_1)-q_{12}(y_1-{\rho }_1)-q_{21}({\rho }_2-x_2)-q_{22}(y_2-{\rho }_2),\end{array}$  

Where $q_{ij}(t)\ge 0$  are the penalty multipliers satisfying $q_{11}(t)({\rho }_1(t)-x_1)=q_{12}(t)(y_1-{\rho }_1(t))=0$  at the optimal ${\rho }_1={\rho }_1^{\ast }$ and $q_{21}(t)({\rho }_2(t)-x_2)=q_{22}(t)(y_2-{\rho }_2(t))=0$  at the optimal ${\rho }_2={\rho }_2^{\ast }$ . The optimal control $({\rho }_1^{\ast },{\rho }_2^{\ast })$  to be determined leads to the following theorem.

 Theorem 3.2: Let ${\rho }_1^{\ast },{\rho }_2^{\ast }$ be the given optimal control to be determine, such that, if the system $U_1^{*},U_2^{*},U_1^{*}{}^{*},U_2^{*}{}^{*},V^{*},P^{*}$ and $M^{\ast }$  are the solutions for the corresponding state variables (1)–( 7), then there exists adjoint variables ${\tau }_i$ , i=1, 2... 7 satisfying

${\tau }_1{}^{\prime }=-1\left\{\begin{array}{l}{\tau }_1[-{\alpha }_1-(1-{\rho }_1^{\ast })g_1(V^{\ast }+P^{\ast })]+{\tau }_3[(1-{\rho }_1^{\ast })g_1(V^{\ast }+P^{\ast })]\\ +{\tau }_5[(1-{\rho }_1^{\ast }){\gamma }_1{g}_1{V}^{\ast }]+{\tau }_6[(1-{\rho }_1^{\ast }){\gamma }_1{g}_1{P}^{\ast }]\end{array}\right\}$  

${\tau }_2{}^{\prime }=-1\left\{\begin{array}{l}{\tau }_2[-{\alpha }_2-(1-\frac{{\rho }_1^{\ast }}{r_1})g_2(V^{\ast }+P^{\ast })]+{\tau }_4[(1-\frac{{\rho }_1^{\ast }}{r_1})g_2(V^{\ast }+P^{\ast })]\\ +{\tau }_5[(1-\frac{{\rho }_1^{\ast }}{r_1}){\gamma }_2{g}_2{V}^{\ast }]+{\tau }_6[(1-\frac{{\rho }_1^{\ast }}{r_1}){\gamma }_2{g}_2{P}^{\ast }]\end{array}\right\}$

${\tau }_3{}^{\prime }=-1\left\{{\tau }_3(-\mu -h_1{M}^{\ast })+{\tau }_7\left(\frac{w_{M}M{H}_w}{{(U_1^{\ast }+U_2^{\ast }+H_w)}^2}-\frac{q_{M}M{H}_q}{{(U_1^{\ast }+U_2^{\ast }+H_q)}^2}\right)\right\}$

${\tau }_4{}^{\prime }=-1\left\{{\tau }_4(-\mu -h_2{M}^{*})+{\tau }_7\left(\frac{w_{M}M{H}_w}{{(U_1^{*}+U_2^{*}+H_w)}^2}-\frac{q_{M}M{H}_q}{{(U_1^{*}+U_2^{*}+H_q)}^2}\right)\right\}$

${\tau }_5{}^{\prime }=-1\left\{\begin{array}{l}{K}_1+{\tau }_5[(1-{\rho }_2^{*})p_v-(\mu +n)-(1-{\rho }_1^{*}){\gamma }_1{g}_1{U}_1^{*}+(1-\frac{{\rho }_1^{*}}{r_1}){\gamma }_2{g}_2{U}_2^{*}]\\ -{\tau }_1[(1-{\rho }_1^{*})g_1{U}_1^{*}]-{\tau }_2[(1-\frac{{\rho }_1^{*}}{r_1})g_1{U}_2^{*}]+{\tau }_3[(1-{\rho }_1^8)g_1{U}_1^{*}]+{\tau }_4[(1-\frac{{\rho }_1^{*}}{r_1})g_2{U}_2^{*}]\end{array}\right\}$

${\tau }_6{}^{\prime }=-1\left\{\begin{array}{l}{K}_2+{\tau }_6[(1-\frac{{\rho }_2^{*}}{r_2})p_p-(\mu +n)-(1-{\rho }_1^{*}){\gamma }_1{g}_1{U}_1^{*}+(1-\frac{{\rho }_1^{*}}{r_1}){\gamma }_2{g}_2{U}_2^{*}]\\ -{\tau }_1[(1-{\rho }_1^{*})g_1{U}_1^{*}]-{\tau }_2[(1-\frac{{\rho }_1^{*}}{r_1})g_2{U}_2^{*}]+{\tau }_3[(1-{\rho }_1^{*})g_1{U}_1^{*}]+{\tau }_4[(1-\frac{{\rho }_1^{*}}{r_1})g_2{U}_2^{*}]\end{array}\right\}$

${\tau }_7=-1\left\{-\delta -{\tau }_3{h}_1{U}_1^{**}-{\tau }_4{h}_2{U}_2^{**}+{\tau }_7\left[\frac{w_M(U_1^{**}+U_2^{**})}{(U_1^{**}+U_2^{**})+H_w}-\frac{q_M(U_1^{**}+U_2^{**})}{(U_1^{**}+U_2^{**})+H_q}-{\mu }_M\right]\right\}$

and having ${\tau }_i(t_f)=0,i=1,2,\mathrm{..},7$  as transversality conditions. Moreover,

${\rho }_1^{\ast }(t)=\min \left\{\max \left\{x_1,\frac{1}{2L_1}({\tau }_1{U}_1^{*}(t))+({\tau }_2{U}_2^{*}(t))\right\},y_1\right\}$  

and ${\rho }_2^{\ast }(t)=\min \left\{\max \left\{x_2,-\frac{{\tau }_3+{\tau }_4+{\tau }_5+{\tau }_6+{\tau }_7}{2L_2}\frac{U_1^{**}(t)+U_2^{**}(t)+V^{*}(t)+P^{*}(t)+M^{*}(t)}{(C_1+C_2+C_3)[U_3^{**}(t)+U_2^{**}(t)+V^{*}(t)+P^{*}(t)+M^{*}(t)}\right\},y_2\right\}$ .