Generating the Figures in Chapter 1 of Otto and Day

Clear["Global`*"]

Figure 1.5 - Number of HIV particles within the bloodstream (based on Phillips (1996) Science 271:497-499)

Clear["Global`*"]

"The model is defined by four equations describing the interrelated changes 
over time in the number of activated, uninfected CD4 lymphocytes (R), 
latently infected cells (L), actively infected cells (E), and free virions (V)."

"These equations can be explained as follows."

"Activated, uninfected CD4 lymphocytes (R) arise at a constant rate Γτ,
where Γ is the rate at which new, uninfected CD4 lymphocytes arise and
τ is the proportion that are activated, and they are removed by
HIV-independent cell death at rate μ or by infection at rate βV"

dr[t]/dt ==  Γ*τ - μ*r[t] - β*r[t]*v[t]

"Upon infection, a proportion p of cells become latently infected (L),
and these are removed by either HIV-independent cell death or by
activation at rate α"

dl[t]/dt == p*β*r[t]*v[t] - μ*l[t] - α*l[t]

"Actively infected cells are generated immediately after infection or 
from the activation of latently infected cells before they die at rate δ"

de[t]/dt == (1 - p)*β*r[t]*v[t] + α*l[t] - δ*e[t]

"Free virions are produced at rate π by actively infected cells and
removed at rate σ"

dv[t]/dt ==πP*e[t] - σ*v[t]

A couple of notes on notation:  We have converted the variables to lower case letters, reserving capital letters for functions that are previously defined in Mathematica.  Also, π is defined Mathematica as 3.14, so to represent the parameter in Phillips' model we use πP.

?NDSolve

Here is where we plug in the parameter values from Phillips' Table 1.

VarSet={Γ->1.36,μ->1.36*10^(-3),τ->0.2,β->0.00027,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->100};

Solution[t_]=
NDSolve[{
	D[r[t],t]==Γ*τ-μ*r[t]-β*r[t]*v[t],
	D[l[t],t]==p*β*r[t]*v[t]-μ*l[t]-α*l[t],
	D[e[t],t]==(1-p)*β*r[t]*v[t]+α*l[t]-δ*e[t],
	D[v[t],t]==πP*e[t]-σ*v[t],
	
	r[0]==1000*τ,
	l[0]==0,
	e[0]==0,
	v[0]==4*10^(-7)}
/.VarSet,
{r[t],l[t],e[t],v[t]},{t,0,120}]

LogPlot[(v[t]/.Solution[t])[[1]],{t,0,120},PlotRange->{Automatic,{0.1,10000}}]

Figure 1.6 - Spread of HIV within the San Francisco gay community (based on Blower et al. (1996) Science 287:650-654)

Clear["Global`*"]

Analysis with no treatment (parameter values from note 15):

VarSetNull={bsu->1/10,vsu->1/12,c->17/10,pi->2133,mu->1/30,sig->0};

This leads to an equilibrium HIV prevalence of about 30% (within rounding errors of the parameters):

NSolve[{0==pi-x[t]*(c*(bsu*ysu[t]+bst*yst[t])/(x[t]+ysu[t]+yst[t])+mu),
0==x[t]*c*(bsu*ysu[t]+bst*yst[t])/(x[t]+ysu[t]+yst[t])+yst[t]*gs-ysu[t]*(sig+vsu+mu),
0==ysu[t]*sig-yst[t]*(gs+vst+mu)}/.VarSetNull/.x[t]->x/.ysu[t]->ysu/.yst[t]->yst,{x,ysu,yst}]//N

ysu/(x+ysu)/.%[[1]]

x0=24611.5;
ysu0=11251;
yst0=0;

To numerically solve the differential equations, we now set the rest of the parameters to zero:

VarSetNull={bsu->1/10,vsu->1/12,c->17/10,pi->2133,mu->1/30,sig->0,gs->0,bst->0,vst->0};

Solution[t_]:=
NDSolve[{
	D[x[t],t]==pi-x[t]*(c*(bsu*ysu[t]+bst*yst[t])/(x[t]+ysu[t]+yst[t])+mu),
	D[ysu[t],t]==x[t]*c*(bsu*ysu[t]+bst*yst[t])/(x[t]+ysu[t]+yst[t])+yst[t]*gs-ysu[t]*(sig+vsu+mu),
	D[yst[t],t]==ysu[t]*sig-yst[t]*(gs+vst+mu),
	x[0]==10000,
	ysu[0]==2000,
	yst[0]==0}
/.VarSetNull,
{x[t],ysu[t],yst[t]},{t,0,250}]

Number of uninfected individuals (pre-treatment era):

Plot[Evaluate[x[t]/.Solution[t],{t,0,250}]]

(Evaluate is used in the above so that the general Solution[t] is used before t is set to 0, 1, 2, ... j.)

Number of infected individuals (pre-treatment era):

Plot[Evaluate[ysu[t]/.Solution[t],{t,0,250}]]

Effects of treatment:  Using median values for each parameter (in addition to note 15, these are found in note 21, 23, 25 (sigma is discussed on p. 651 and the increase in risky behavior, ie c is discussed in  note17.)

Clear[VarSet,Solution,i,tablewithout,tablewith,AIDSaverted]

VarSet[i_]={bsu->1/10,vsu->1/12,c->(1+i)*17/10,pi->2133,mu->1/30,gs-> 1/20, bst->1/4*1/10, vst->1/27,sig->1/2};

Solution[t_,i_]:=
NDSolve[{
	D[x[t],t]==pi-x[t]*(c*(bsu*ysu[t]+bst*yst[t])/(x[t]+ysu[t]+yst[t])+mu),
	D[ysu[t],t]==x[t]*c*(bsu*ysu[t]+bst*yst[t])/(x[t]+ysu[t]+yst[t])+yst[t]*gs-ysu[t]*(sig+vsu+mu),
	D[yst[t],t]==ysu[t]*sig-yst[t]*(gs+vst+mu),
	x[0]==x0,
	ysu[0]==ysu0,
	yst[0]==yst0}
/.VarSet[i],
{x[t],ysu[t],yst[t]},{t,0,250}]

Numbers of HIV- (thin curve), untreated HIV+ (short dashed curve), treated HIV+ (long dashed curve), total HIV+ (thick curve) if behavior remains unchanged:

Show[
Plot[Evaluate[x[t]/.Solution[t,0],{t,0,250}],PlotRange->{{0,20},{0,70000}}],
Plot[Evaluate[ysu[t]/.Solution[t,0],{t,0,250}],PlotStyle->Dashing[{0.01}],PlotRange->{{0,20},{0,70000}}],
Plot[Evaluate[yst[t]/.Solution[t,0],{t,0,250}],PlotStyle->Dashing[{0.01,0.05}],PlotRange->{{0,20},{0,70000}}],
Plot[Evaluate[yst[t]+ysu[t]/.Solution[t,1],{t,0,250}],PlotStyle->Thickness[0.01],PlotRange->{{0,20},{0,70000}}]
]

Numbers of HIV- (thin curve), untreated HIV+ (short dashed curve), treated HIV+ (long dashed curve), total HIV+ (thick curve) if risky behavior increases by 50%:

Show[
Plot[Evaluate[x[t]/.Solution[t,0.5],{t,0,250}],PlotRange->{Automatic,{0,70000}}],
Plot[Evaluate[ysu[t]/.Solution[t,0.5],{t,0,250}],PlotStyle->Dashing[{0.01}],PlotRange->{Automatic,{0,70000}}],
Plot[Evaluate[yst[t]/.Solution[t,0.5],{t,0,250}],PlotStyle->Dashing[{0.01,0.05}],PlotRange->{Automatic,{0,70000}}],
Plot[Evaluate[yst[t]+ysu[t]/.Solution[t,0.5],{t,0,250}],PlotStyle->Thickness[0.01],PlotRange->{Automatic,{0,70000}}]
]

Numbers of HIV- (thin curve), untreated HIV+ (short dashed curve), treated HIV+ (long dashed curve), total HIV+ (thick curve) if risky behavior doubles:

Show[
Plot[Evaluate[x[t]/.Solution[t,1],{t,0,250}],PlotRange->{{0,20},{0,70000}}],
Plot[Evaluate[ysu[t]/.Solution[t,1],{t,0,250}],PlotStyle->Dashing[{0.01}],PlotRange->{{0,20},{0,70000}}],
Plot[Evaluate[yst[t]/.Solution[t,1],{t,0,250}],PlotStyle->Dashing[{0.01,0.05}],PlotRange->{{0,20},{0,70000}}],
Plot[Evaluate[yst[t]+ysu[t]/.Solution[t,1],{t,0,250}],PlotStyle->Thickness[0.01],PlotRange->{{0,20},{0,70000}}]
]

Without treatment the number of AIDS deaths would be:

tablewithout=Table[Sum[vsu*ysu0/.VarSet[0],{t,0,j}],{j,0,10}]//N

tableaxes=Table[j,{j,0,10}]

With treatment but with a (1+i) factor change in the amount of risky sexual behavior, the cumulative number of AIDS deaths from the beginning of treatment to year j would be:

tablewith[i_]:=tablewith[i]=Table[Sum[Evaluate[vsu*ysu[t]+vst*yst[t]/.Solution[t,i]/.VarSet[i]],{t,0,j}][[1]],{j,0,10}]

tablewith[0.5]

Percent decrease in the number of AIDS deaths:

Clear[AIDSaverted]
	AIDSaverted[i_]:=AIDSaverted[i]=Transpose[{tableaxes,100*(tablewithout-tablewith[i])/tablewithout}]

Show[ListPlot[AIDSaverted[0.],Joined->True],ListPlot[AIDSaverted[0.25],Joined->True],ListPlot[AIDSaverted[0.5],Joined->True],ListPlot[AIDSaverted[0.75],Joined->True],ListPlot[AIDSaverted[1],Joined->True],PlotRange->{{-0.2,10.5},{0,50}}]

If risky behavior increases more:

Show[ListPlot[AIDSaverted[0.],Joined->True],ListPlot[AIDSaverted[0.5],Joined->True],ListPlot[AIDSaverted[1],Joined->True],ListPlot[AIDSaverted[2],Joined->True],ListPlot[AIDSaverted[3],Joined->True],PlotRange->{{-0.2,10.5},{-40,100}}]

Mathematica’s manipulate function allows us to vary the change in risky behavior (i) dynamically if desired, to see how this alters the percentage of cases that are averted.  Note that when this becomes negative, the change in risky behavior cancels the benefits of treatment in terms of the number of AIDS deaths.
[Click the + sign by the toggle bar to see the exact values of i used):

Manipulate[ListPlot[AIDSaverted[i],Joined->True,PlotRange->{{-0.2,10.5},{-50,50}}],{{i,0},0,3}]

Without treatment the cumulative number of HIV cases would be:

tablewithout2=Table[Sum[x0*(c*(bsu*ysu0+bst*yst0)/(x0+ysu0+yst0))/.VarSet[0],{t,0,j}],{j,0,10}]//N

With treatment but with a (1+i) factor change in the amount of risky sexual behavior, the number of HIV cases would be:

tablewith2[i_]:=tablewith2[i]=Table[Sum[Evaluate[x[t]*(c*(bsu*ysu[t]+bst*yst[t])/(x[t]+ysu[t]+yst[t]))/.Solution[t,i]/.VarSet[i]],{t,0,j}][[1]],{j,0,10}]

tablewith2[0.5]

Percent decrease in the cumulative number of new HIV cases:

NEWaverted[i_]:=NEWaverted[i]=Transpose[{tableaxes,100*(tablewithout2-tablewith2[i])/tablewithout2}]

Show[ListPlot[NEWaverted[0.],Joined->True],ListPlot[NEWaverted[0.25],Joined->True],ListPlot[NEWaverted[0.5],Joined->True],ListPlot[NEWaverted[0.75],Joined->True],ListPlot[NEWaverted[1],Joined->True],PlotRange->{{-0.2,10.5},{-100,100}}]

Mathematica’s manipulate function allows us to vary the change in risky behavior (i) dynamically if desired, to see how this alters the percentage of cases that are averted.  Note that when this becomes negative, the change in risky behavior cancels the benefits of treatment in terms of the cumulative number of new HIV cases.
[Click the + sign by the toggle bar to see the exact values of i used):

Manipulate[ListPlot[NEWaverted[i],Joined->True,PlotRange->{{-0.2,10.5},{-100,100}}],{{i,0},0,3}]

Figure 1.7 - Spread of HIV within the Hlabisa community of women (based on Williams et al. (2001) Statistics in Medicine 20:2003-2016)

Clear["Global`*"]

From Williams et al (2001) Table 1:

total={7,32,86,206,276,219,241,205,193,203,152,148,117,130,123,90,119,62,99,60,58,71,55,47,65,27,29,13,8,3,4,8,1,3,1,2};

infected={0,2,11,31,70,59,79,84,81,82,65,69,39,46,39,28,30,20,21,8,14,17,12,9,16,7,3,2,2,0,0,2,0,0,0,0};

prevalence=infected/total

To help get rid of noise in the prevalence data, Williams et al smooth it using piecewise polynomials to the log-transformed data, resulting in the following smoothed estimates for prevalence:

smoothed={0.016,0.047,0.102,0.173,0.244,0.304,0.348,0.377,0.392,0.397,0.395,0.387,0.376,0.362,0.346,0.329,0.311,0.293,0.275,0.256,0.238,0.220,0.203,0.187,0.171,0.156,0.142,0.128,0.116,0.105,0.094,0.085,0.076,0.068,0.061,0.054};

Here are lists of the relevant ages:

age=Table[i,{i,14,49}]

age2=Table[i,{i,15,49}]

To get a table with {x,y} coordinates of {age,prevalence}, we put together the two lists {{age},{prevalance}} and then transpose this table (making rows into columns and vice versa):

Transpose[{age,prevalence}]

Show[ListPlot[Transpose[{age,prevalence}],PlotStyle->PointSize[0.02]],ListPlot[Transpose[{age,smoothed}],Joined->True],PlotRange->{{10,50},{0,0.5}},AxesOrigin->{10,0}]

Spread rate from 1997 to 1998:

spread=1.45;

"Various studies indicate that there is very little mortality over the first two years after which survivorship starts to drop with a median life-expectancy after infection of about seven years corresponding to an average mortality of
about 10 per cent per year [12–14]."

Noting that the median survival is seven years, we can ask what the survival rate must be per year such that 50% of individuals survive for 7 years:

Solve[survive^7==0.5,survive]

The only real solution is:

survival=0.9057236642639066

The risk is estimated by equation (S10.2.8) in the supplementary material to Chatper 10:

incidence=Table[(smoothed[[i]](1-(1-survival)*smoothed[[i-1]]/spread)-smoothed[[i-1]]/spread*survival)/(1-smoothed[[i-1]]/spread),{i,2,Length[smoothed]}]

Show[ListPlot[Transpose[{age,prevalence}],PlotStyle->PointSize[0.02]],ListPlot[Transpose[{age,smoothed}],Joined->True],ListPlot[Transpose[{age2,incidence}],Joined->True,PlotStyle->Dashing[{0.01,0.01}]],PlotRange->{{10,50},{0,0.5}},AxesOrigin->{10,0}]

Next, we estimate the relative amount of error if we were to use the shape of the prevalence curve (the data) to infer the shape of the incidence curve (giving the unknown risk).  To do this, we first normalize the above curves, so that the areas under them are equal (in this case, the sums over each age class are equal):

sum1=Sum[incidence[[i]],{i,1,Length[incidence]}]

sum2=Sum[smoothed[[i+1]],{i,1,Length[incidence]}]

relative=Table[(incidence[[i]]/sum1)/(smoothed[[i+1]]/sum2),{i,1,Length[incidence]}]

Show[
ListPlot[Transpose[{age2,relative}]],
Plot[1,{x,10,49}],
PlotRange->{{10,50},{0,2}},
AxesOrigin->{10,0}
]