# Primary Infection by HIV
Phillips (1996) Science 271:497-499

"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 Γ tauτ,
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 == π E[t] - σ V[t]

S;

Hint: Type "Esc" then spell out a greek letter like "mu" then "Esc" to format a greek letter.

## Equations

Type the set of equations:
Hint: Use command-1 to insert a comment.

eqns: [
    dRdt :  Γ*τ - μ*R- β*R*V,
    dLdt : p*β*R*V - μ*L - α*L,
    dEdt :  (1-p)*β*R*V + α*L - δ*E,
    dVdt : π*E - σ*V];

Variables:

vars : [R, L, E, V];

Maxima knows many of the rules of algebra and calculus and allows us to focus on the insights that math provides.  For
example, we could calculate the derivative of an equation with respect to a variable

diff(eqns[1],V);

## Calculating an equilibrium

As we'll learn later in class, an equilibrium is a point where the system does not change over time. 

Maxima can be used to solve systems of equations like these for when nothing changes:

equil:solve([ 0= Γ*τ - μ*R- β*R*V,
    0=p*β*R*V - μ*L - α*L,
    0=  (1-p)*β*R*V + α*L - δ*E,
    0= π*E - σ*V],vars);

The above shows two equilibria - how do they differ?

Let's look at how the equilibrium density of viruses changes with increasing infectivity β:

part(equil,2,4,2);

toplot(β):=subst([Γ = 1.36, μ =1.36*10^(-3),τ = 0.2, p = 0.1, α= 3.6*10^(-2), 
        σ= 2, δ= 0.33, π = 100], part(equil,2,4,2))$

wxplot2d([toplot(β)], [β, 0.00001, 0.00030],[ylabel,"Equilibrium virus"]);

## Example plot

Here we susbtitute in some parameter values into the equations:

example:subst([Γ = 1.36, μ =1.36*10^(-3),τ = 0.2,β = 0.00027, p = 0.1, α= 3.6*10^(-2), 
        σ= 2, δ= 0.33, π = 100], eqns);

Here we give an initial starting point for the set of variables (in the order defined by vars):

plotstart : [200, 0, 0,  4*10^(-7)];

Then we plug the equations, the variables, and the starting point into numerical solver of differential equations
(rkf45):

load(rkf45);

 sol:rkf45(example,vars,plotstart,[t,0,100],
             report=true,absolute_tolerance=1e-8)$

Here we pull out the list from the solution for each variable ("first" is time, then each of the variables in our list,
starting with susceptibles):

time:makelist(part(sol,t,1),t,1,length(sol))$
susceptibles: makelist(part(sol,t,2),t,1,length(sol))$
active: makelist(part(sol,t,3),t,1,length(sol))$
latent: makelist(part(sol,t,4),t,1,length(sol))$
virus: makelist(part(sol,t,5),t,1,length(sol))$

Plotting the curves, coloured according to the flow diagram (blue for susceptibles, pink for exposed, red for
presymptomatic, brown for infected, and orange for asymptomatic):

wxdraw2d(point_type=-1,points_joined=true,color=blue,points(time,susceptibles),
    color=pink,points(time,active),color=red,points(time,latent),
    color=brown,points(time,virus),
    xlabel="time",ylabel="type");

Alternatively, we could plot this on a log-scale along the y-axis (as did Philips):

wxdraw2d(logy=true,point_type=-1,points_joined=true,color=blue,points(time,susceptibles),
    color=pink,points(time,active),color=red,points(time,latent),
    color=brown,points(time,virus),
    xlabel="time",ylabel="type");