# Solving Differential Equation System for HIV Treatment Model

I was working on a project about optimal strategies for HIV treatment, models used from [Butler, Kirschner, and Lenhart] 1997. This model explains the spread of HIV viruses in the human body, where there is one control function u(t). My work is following pontryagin maximum principle. But i have a problem solving the differential equation system, where there are 6 differential equations with 6 initial conditions. Here i use software Mathematica 11 when I execute the code there is no error that appears and he just returning DSolve code that i write.

ClearAll["Global*"];
dT = s/(1 + V[t]) - m1 T[t] + r T[t] (1 - (T[t] + Ti[t])/Tm) -
u[t] k V[t] T[t];
dTi = u[t] k V[t] T[t] - m2 Ti[t];
dV = M m2 Ti[t] - m3 V[t];
H = A T[t] - (1 - u[t])^2 + l1[t] dT + l2[t] dTi + l3[t] dV ;
u[t] = 1/2 (2 - k l1[t] T[t] V[t] + k l2[t] T[t] V[t])
eq1 = D[l1[t], t] == D[H, l1[t]];
eq2 = D[l2[t], t] == D[H, l2[t]];
eq3 = D[l3[t], t] == D[H, l3[t]];
eq4 = D[T[t], t] == -D[H, T[t]];
eq5 = D[Ti[t], t] == -D[H, Ti[t]];
eq6 = D[V[t], t] == -D[H, V[t]];
s = 10; m1 = 0.02; m2 = 0.5; m3 = 4.4; r = 0.03; Tm = 1500; k = \
0.000024; M = 300; A = 1;

DSolve[{eq1, eq2, eq3, eq4, eq5, eq6, T == 800, Ti == 0.04,
V == 1.5, l1 == 0, l2 == 0, l3 == 0}, {l1[t], l2[t],
l3[t], T[t], Ti[t], V[t]}, t]


I cannot guarantee that the analytical solution exists, but is there something wrong with code I wrote? Or is there an alternative solution to complete the differential equation system?

I hope you are pleased to check the file that I attached. Control-Optimal.nb

Many thanks • Thanks for your comment, it give this error: NDSolve::deqn: Equation or list of equations expected instead of eq1 in the first argument {eq1,eq2,eq3,eq4,eq5,eq6,T==800,Ti==0.04,V==1.5,l1==0,l2==0,l3==0}. Dec 24, 2018 at 13:25
• I really appreciate it, thank you Dec 24, 2018 at 13:32
• I'm sure there is no analytical solution, so you'll need to use NDSolve. I noticed that the equations for T'[t], Ti'[t] and V'[t] in your code do not match the system you want to solve. In any case, you might want to include all your code in the question in case people do not want to download a notebook. Dec 24, 2018 at 14:42
• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Dec 24, 2018 at 14:46
• Oh i see, i will focus on NDSolve, thanks for your advice @ChrisK Dec 24, 2018 at 23:27

Structuring the script.

parms = {s -> 10, m1 -> 0.02, m2 -> 0.5, m3 -> 4.4, r -> 0.03, Tm -> 1500, k -> 0.000024, M -> 300, A -> 1, T0 -> 800, Ti0 -> 0.04, V0 -> 1.5};
x1 = T[t];
x2 = Ti[t];
x3 = V[t];
X = {x1, x2, x3};
Lambda = {la1[t], la2[t], la3[t]};
U = u[t];
f1 = s/(1 + x3) - m1 x1 + r x1 (1 - (x1 + x2)/Tm) - U k x1 x3;
f2 = U k x3 x1 - m2 x2;
f3 = M m2 x2 - m3 x3;
F = {f1, f2, f3};
L = A x1 - (1 - U)^2;
H = L + Lambda.F;
equ3 = D[H, U] == 0;
solU = Solve[equ3, U][];
equs = Join[equ1, equ2];
cinitsX = Thread[X == {T0, Ti0, V0}] /. {t -> 0};
cinitsLambda = Thread[Lambda == 0] /. {t -> 0};
cinits = Join[cinitsX, cinitsLambda];
vars = Join[X, Lambda];
DE = Join[equs, cinits] /. solU /. parms;
solDE = NDSolve[DE, vars, {t, 0, 0.01}];
Plot[Evaluate[vars /. solDE], {t, 0, 0.01}]

• thankyou very much, its help me to learn this algorithm a lot Jan 3, 2019 at 17:52

Cesario's answer is almost there, but I'd like to add a couple points of correction/improvement.

1. The problem formulation has constraints on the control (0<=u<=1) and positivity constraints on the states. I would suggest rather than solving for root where the derivative is zero (equ3 above), use Maximize and add the constraints there:
parms = {s -> 10, m1 -> 0.02, m2 -> 0.5, m3 -> 4.4, r -> 0.03,
Tm -> 1500, k -> 0.000024, M -> 300, A -> 1, T0 -> 800,
Ti0 -> 0.04, V0 -> 1.5};
x1 = T[t];
x2 = Ti[t];
x3 = V[t];
X = {x1, x2, x3};
Lambda = {la1[t], la2[t], la3[t]};
U = u[t];
f1 = s/(1 + x3) - m1 x1 + r x1 (1 - (x1 + x2)/Tm) - U k x1 x3;
f2 = U k x3 x1 - m2 x2;
f3 = M m2 x2 - m3 x3;
F = {f1, f2, f3};
L = A x1 - (1 - U)^2;
H = L + Lambda . F;
(*equ3=D[H,U]==0;
solU=Solve[equ3,U][];*)
cons = {0 <= u[t] <= 1 && T[t] > 0 && Ti[t] > 0 && V[t] > 0, k > 0};
equ3 = Last@Maximize[{H, cons}, U];
solU = FullSimplify[equ3, cons]


This (nearly) yields the piecewise function reported in the original work. Please note that the third ODE in the SE problem is missing a term found in the original problem, but the idea is the same. 1. The boundary conditions for lambda are given as lambda in Cesario's answer, but these should apply to the end time tf, which is the gradient of the end cost (Mayer term) and here {0,0,0}. Programmatically:
M=0; (*Mayer cost*)
tf=100;  (* from paper *)
cinitsLambda = Thread[(Lambda /. {t -> tf}) == Grad[M, X /. t -> tf]]; Now the setup for NDSolve can be completed:

equs = Join[equ1, equ2];
cinits = Join[cinitsX, cinitsLambda];
vars = Join[X, Lambda];
DE = Join[equs, cinits] /. solU /. parms;

1. And you might hope that the following NDSolve statement would work:
solDE = NDSolve[DE, vars, {t, 0, tf}];
Plot[Evaluate[vars /. solDE], {t, 0, tf}]
`

But alas it does not seem to converge in an overnight run, nor does the original equations from the paper. (vide supra)

I know that BVPs are tough problems to solve, so I tried an alternative implementation in my personal fork of PSOPT pseudospectral code, as well as the trapezoidal method in Mathematica. These are direct solvers (vs. the Pontryagin indirect method) and also have problems with these equations. I have come to conclusion that there's something fishy or extreme stiffness.

I have tried this code on numerous multivariate optimal control problems and it seems to yield answers from Ross, Bryson and Ho, et al.