Application Center - Maplesoft

App Preview:

HIV and antiviral therapy

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

TR-00-22[1].mws

HIV and Antiviral Therapy

Robert M. Corless
Ontario Research Centre for Computer Algebra
http://www.orcca.on.ca

An exploration of "Mathematical Models and the design of public health policy: HIV and Antiviral therapy", by Sunetra Gupta, Roy M. Anderson, and Robert M. May, SIAM Reviw, Vol 35, No. 3, pp. 1-16, March 1993.

Notation of the paper:

p = proportion of cohort

mu*N[0] = sexually active population

V[x] = vaccinated susceptibles

X = unvaccinated susceptibles

s = per capita rate of vaccination

lambda = per partnership probability of loss

c = mean rate of partner change

N[0] = stable population size in the absence of infection

mu = inverse of average duration of sexual activity

Y = untreated infected population

v = per capita rate of development of AIDS

V[y] = number of infected individuals receiving therapy

d = per capita rate of treated infected population development of AIDS

B[1] = probability of transmission from an untreated carrier of HIV

B[2] = probability of transmission from a treated carrier.

N = current size of the sexually active population

> restart;

The assumptions of the exposition in the paper (we can come back and change this later and see how the conclusions change):

> s := 0;

s := 0

> r := 0;

r := 0

The basic equations in differential form (with the dependence on time suppressed for ease of working with equilibria):
the rate of change of unvaccinated susceptibles depends on how many come of age, and how many are infected, and how many stop sexual activity, and how many are vaccinated:

> eq1 := diff(X(t),t) = mu*N[0]*(1-p) - (c*lambda + mu + s)*X;

eq1 := diff(X(t),t) = mu*N[0]*(1-p)-(c*lambda+mu)*X...

Similarly, the rate of change of vaccinated susceptibles depends on how many come of age and are vaccinated, and how many are infected, and how many stop sexual activity, and how many are vaccinated from the X population.

> eq2 := diff(V[x](t),t) = mu*N[0]*p-(c*lambda+mu)*V[x] + s*X;

eq2 := diff(V[x](t),t) = mu*N[0]*p-(c*lambda+mu)*V[...

> eq3 := diff(Y(t),t) = c*lambda*X-(v+r+mu)*Y;

eq3 := diff(Y(t),t) = c*lambda*X-(v+mu)*Y

> eq4 := diff(V[y](t),t) = c*lambda*V[x] + r*Y - (d+mu)*V[y];

eq4 := diff(V[y](t),t) = c*lambda*V[x]-(d+mu)*V[y]

The paper doesn't say so, but the calculation only works out if we say that the total population is just the sum of the parts.

> NN := X + V[x] + Y + V[y];

NN := X+V[x]+Y+V[y]

The per partnership probability of infection depends on the infectiousness and the likelihood that the partner is infected:

> lambda := (B[1]*Y + B[2]*V[y])/NN;

lambda := (B[1]*Y+B[2]*V[y])/(X+V[x]+Y+V[y])

We only care about the ultimate state, when the populations don't change: therefore we set all derivatives to zero.

> equilibria := map(rhs,{eq1,eq2,eq3,eq4});

equilibria := {c*(B[1]*Y+B[2]*V[y])*X/(X+V[x]+Y+V[y...
equilibria := {c*(B[1]*Y+B[2]*V[y])*X/(X+V[x]+Y+V[y...

Maple can solve those population equations, to give us solutions that look a little different than those in the paper but that we can experiment with:

> sols := solve(equilibria, {X,Y,V[x],V[y]});

sols := {V[y] = (-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p...
sols := {V[y] = (-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p...
sols := {V[y] = (-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p...
sols := {V[y] = (-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p...
sols := {V[y] = (-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p...
sols := {V[y] = (-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p...
sols := {V[y] = (-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p...
sols := {V[y] = (-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p...

One solution occurs if no-one is infected to start with, so Y = V[y] = 0. It is the other one we are interested in. It is curious that there is no solution if X=V[x]=0; this is a result of an assumption in the paper that there is always a population of uninfected children "coming of age" to sexual maturity to supply uninfected susceptibles. (This is not pointed out in the paper).

Now let us check that the solutions we get are the same as those in the paper.

> assign(sols[1]);

> Y;

-mu*N[0]*(-mu^2+mu^2*p-c*B[2]*p^2*mu+c*B[1]*p^2*mu-...
-mu*N[0]*(-mu^2+mu^2*p-c*B[2]*p^2*mu+c*B[1]*p^2*mu-...
-mu*N[0]*(-mu^2+mu^2*p-c*B[2]*p^2*mu+c*B[1]*p^2*mu-...
-mu*N[0]*(-mu^2+mu^2*p-c*B[2]*p^2*mu+c*B[1]*p^2*mu-...

This doesn't look like the solution in the paper. We should see if we can simplify the computed Y to see if there is an error.

> Y/c/N[0]/mu/(1-p)*(v+mu); # = lamstar/(c*lamstar+mu)

-(-mu^2+mu^2*p-c*B[2]*p^2*mu+c*B[1]*p^2*mu-2*c*B[1]...
-(-mu^2+mu^2*p-c*B[2]*p^2*mu+c*B[1]*p^2*mu-2*c*B[1]...
-(-mu^2+mu^2*p-c*B[2]*p^2*mu+c*B[1]*p^2*mu-2*c*B[1]...
-(-mu^2+mu^2*p-c*B[2]*p^2*mu+c*B[1]*p^2*mu-2*c*B[1]...

We should get something like lambda/(c*lambda+mu), for a particular value of lambda---let us see if we can get something like the paper does.

> solve(lstar/(c*lstar+mu)=%,lstar);;

-(-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p+c*B[1]*mu*p+mu...

> lstar := %;

lstar := -(-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p+c*B[1...

> lstarf := map(factor,%);

lstarf := -(-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p+c*B[...

> phieq := B[1]*(1-p)/(v+mu) + B[2]*p/(d+mu);

phieq := B[1]*(1-p)/(v+mu)+B[2]*p/(d+mu)

> phi :='phi';

phi := 'phi'

> solve(phieq=phi,B[1]);

-(-B[2]*p+phi*d+phi*mu)*(v+mu)/((d+mu)*(p-1))

> subs(B[1]=%,lstarf);

-(-c*B[2]*p*v-c*(-B[2]*p+phi*d+phi*mu)*(v+mu)*d*p/(...
-(-c*B[2]*p*v-c*(-B[2]*p+phi*d+phi*mu)*(v+mu)*d*p/(...

> collect(%,phi);

-(-c*(v+mu)*mu*p/(p-1)-c*(v+mu)*d*p/(p-1)+c*(v+mu)*...
-(-c*(v+mu)*mu*p/(p-1)-c*(v+mu)*d*p/(p-1)+c*(v+mu)*...

> collect(%,phi,factor);

(d+mu)*(v+mu)*phi/(mu+p*v-d*p+d)-(d+mu)*(v+mu)/((mu...

This agrees with equations (6) and (8) of the paper (if we don't use N=X+V[x]+Y+V[y]) and even if we do it looks as though it is in the right form.

> phi := phieq;

phi := B[1]*(1-p)/(v+mu)+B[2]*p/(d+mu)

> lamstar := %%;

lamstar := (d+mu)*(v+mu)*(B[1]*(1-p)/(v+mu)+B[2]*p/...

Now we compare our computed answer with the formula from the paper by subtracting the one from the other; if the difference is zero, then we know the paper result is correct.

> Y - c*lamstar*mu*N[0]*(1-p)/(c*lamstar+mu)/(v+mu):

> normal(%);

0

It is, so at least the formula for Y is right. But the formula for V[y] also looks different from that of the paper:

> V[y];

(-c*B[2]*p*v+c*B[1]*d*p-c*B[2]*mu*p+c*B[1]*mu*p+mu^...

Again we subtract the formula from the paper and see if it is zero.

> normal( V[y] - c*lamstar*mu*N[0]*p/(c*lamstar+mu)/(d+mu) );

0

It is. Now we see if we can recapitulate their argument to get the equlibrium population in the presence of infection:

> normal( lamstar - c*lamstar*mu*N[0]*phi/(c*lamstar+mu)/Nstar):

> solve(%,Nstar):

> factor(%);

(mu+p*v-d*p+d)*(-mu*B[2]*p+mu*B[1]*p-mu*B[1]-p*B[2]...

> Nstar := %;

Nstar := (mu+p*v-d*p+d)*(-mu*B[2]*p+mu*B[1]*p-mu*B[...

Again, that formula looks different from the one in the paper, which defines things in terms of an auxiliary parameter, Omega.

> Nstar = N[0]*(1-Omega)/(1-(1/c)*Omega/phi);

(mu+p*v-d*p+d)*(-mu*B[2]*p+mu*B[1]*p-mu*B[1]-p*B[2]...
(mu+p*v-d*p+d)*(-mu*B[2]*p+mu*B[1]*p-mu*B[1]-p*B[2]...

> tst := solve(%,Omega);

tst := -(mu*p*v-v*mu-d*mu*p-d*v)/(d*v+d*mu+v*mu+mu^...

That again looks different from the formula in the paper:

> tst - v*(1-p)/(v+mu) - d*p/(d+mu);

-(mu*p*v-v*mu-d*mu*p-d*v)/(d*v+d*mu+v*mu+mu^2)-v*(1...

> normal(%);

0

But it is not. Therefore, all the formulas in the paper are correct.

Now, we can extend the printed results and work out what the uninfected population NA would be:

> eval(Nstar,p=0):

> NA := factor(%);

NA := N[0]*c*B[1]*mu/((c*B[1]-v)*(v+mu))

> effectiveness_ratio := Nstar/NA;

effectiveness_ratio := (mu+p*v-d*p+d)*(-mu*B[2]*p+m...
effectiveness_ratio := (mu+p*v-d*p+d)*(-mu*B[2]*p+m...

> eval(%,p=1);

(v+mu)*(-mu*B[2]-B[2]*v)*(c*B[1]-v)/((d+mu)*(-c*B[2...

> normal(%);

(c*B[1]-v)*B[2]*(v+mu)/(B[1]*(c*B[2]-d)*(d+mu))

This short exploration showed an example of using Maple to verify the results in an area of mathematical biology and public health policy.