Solving Differential Equation System for HIV Treatment Model

Question

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[0] == 800, Ti[0] == 0.04, 
V[0] == 1.5, l1[20] == 0, l2[20] == 0, l3[20] == 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

Answer 1

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;
equ1 = Thread[D[X, t] == Grad[H, Lambda]];
equ2 = Thread[D[Lambda, t] == -Grad[H, X]];
equ3 = D[H, U] == 0;
solU = Solve[equ3, U][[1]];
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}]

Answer 2

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;
equ1 = Thread[D[X, t] == Grad[H, Lambda]];
equ2 = Thread[D[Lambda, t] == -Grad[H, X]];
(*equ3=D[H,U]==0;
solU=Solve[equ3,U][[1]];*)
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.

2. The boundary conditions for lambda are given as lambda[0] 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;

3. 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.