Four-variable HIV model.ipynb

Kernel: SageMath 9.4

Simulation of the 4-variable HIV model, from Further Exercises 12 and 13 in Section 1.4 of Modeling Life.

Equations:

$$\begin{cases} V' = 100E - 2V \\ R' = 0.272 - 0.00136R - \underline{\qquad} \\ L' = \underline{\qquad} - 0.00136L - \underline{\qquad} \\ E' = \underline{\qquad} + \underline{\qquad} - 0.33E \end{cases}$$

$V = \text{number of viruses}$

$R = \text{number of uninfected T-cells}$

$L = \text{number of latently infected T-cells}$

$E = \text{number of actively infected T-cells}$

import numpy

# Define our state variables:
state_vars = list(var("V,R,L,E"))

# Define the vector field for our system of differential equations:
system = (
    100*E - 2*V,
    0.272 - 0.00136*R - 0.00027*R*V,
    0.1*0.00027*R*V - 0.00136*L - 0.036*L,
    0.9*0.00027*R*V + 0.036*L -0.33*E,
)

# Initial state: 0.0000004 viruses, 200 uninfected T-cells, 0 infected T-cells
initial_state = (4E-7, 200, 0, 0)

# Set up to run the simulation up to t = 400 (days)
t_range = srange(0, 400, 0.01)

# Run the simulation
solution = desolve_odeint(system, initial_state, t_range, state_vars)
solution = numpy.insert(solution, 0, t_range, axis=1)

# Graph the results (time series)
V_plot = list_plot(solution[::100,(0,1)] * (1, 0.1), plotjoined=True,
        color="blue", legend_label="Viruses ($V$, $\\frac{1}{10}$ scale)")
R_plot = list_plot(solution[::100,(0,2)], plotjoined=True,
        color="red", legend_label="Uninfected cells ($R$)")
L_plot = list_plot(solution[::100,(0,3)], plotjoined=True,
        color="gold", legend_label="Latently infected cells ($L$)")
E_plot = list_plot(solution[::100,(0,4)], plotjoined=True,
        color="purple", legend_label="Actively infected cells ($E$)")

# Display the graph
show(V_plot + R_plot + L_plot + E_plot, ymax=550,
        axes_labels=("$t$", "Populations\n(per mm$^3$ of blood)"))

If you want to play around with this simulation, you can copy and paste the code below into a worksheet of your own, and run it. It will give you an interactive that allows you to change the initial values of the three variables, or to run the simulation for a longer period than 400 days.

import numpy

@interact
def HIV_interactive(initV=slider(1E-7, 1E-6, default=4E-7, label="Initial $V$"),
                    initR=slider(100, 1000, default=200, label="Initial $R$"),
                    initL=slider(0, 100, default=0, label="Initial $L$"),
                    initE=slider(0, 100, default=0, label="Initial $E$"),
                    beta=slider(0.1*0.00027, 10*0.00027, 0.00027, default=0.00027, label="$\\beta$"),
                    alpha=slider(0.00036, 0.036, 0.00036, default=0.0036, label="$\\alpha$"),
                    tmax=slider(10, 3650, default=400, label="Simulation length")):
    # Define our state variables:
    state_vars = list(var("V,R,L,E"))

    # Define the vector field for our system of differential equations:
    system = (
        100*E - 2*V,
        0.272 - 0.00136*R - beta*R*V,
        0.1*beta*R*V - 0.00136*L - alpha*L,
        0.9*beta*R*V + alpha*L -0.33*E,
    )

    # Set the initial state
    initial_state = (initV, initR, initL, initE)

    # Set up to run the simulation up to t = 400 (days)
    t_range = srange(0, tmax, 0.01)

    # Run the simulation
    solution = desolve_odeint(system, initial_state, t_range, state_vars)
    solution = numpy.insert(solution, 0, t_range, axis=1)

    # Graph the results (time series)
    V_plot = list_plot(solution[::100,(0,1)] * (1, 0.1), plotjoined=True,
            color="blue", legend_label="Viruses ($V$, $\\frac{1}{10}$ scale)")
    R_plot = list_plot(solution[::100,(0,2)], plotjoined=True,
            color="red", legend_label="Uninfected cells ($R$)")
    L_plot = list_plot(solution[::100,(0,3)], plotjoined=True,
            color="gold", legend_label="Latently infected cells ($L$)")
    E_plot = list_plot(solution[::100,(0,4)], plotjoined=True,
            color="purple", legend_label="Actively infected cells ($E$)")

    # Display the graph
    show(V_plot + R_plot + L_plot + E_plot, ymax=550,
            axes_labels=("$t$", "Populations\n(per mm$^3$ of blood)"))

from sage.calculus.desolvers import desolve_odeint
@interact
def f(t1=text_control("Initial amounts"), 
      initV=slider(10^-7, 10^-6, step_size=10^-7, default=4*10^-7, label="Virus"), 
      initR=slider(100, 1000, default=200, label="Uninfected T-cells"), 
      initL=slider(0,100, label="Latently infected T-cells"), 
      initE=slider(0,100, label="Actively infected T-cells"), 
      tblank=text_control(""), t2=text_control("Parameters"), 
      beta=slider(27/10000, 10*27/10000, step_size=27/100000, default=27/10000, label='beta'), 
      alpha=slider(36/10000, 36/100, step_size=36/10000, default=36/10000, label='alpha'), 
      tmax=slider(10,3650, default=400, label="Simulation length"), 
      showV=checkbox(label="Plot virus levels")):
    
    #Integration
    var("V,R,L,E")
    eqns = [100*E-2*V, 0.272-0.00136*R-beta*R*V, 0.1*beta*R*V - 0.00136*L - alpha*L, 0.9*beta*R*V + alpha*L -0.33*E]
    inits = [initV, initR, initL, initE]
    t=srange(0,tmax,0.025)
    sol=desolve_odeint(eqns, inits, t, [V,R,L,E])

    #Plotting
    #Pick every hundredth element (for speed)
    dect = t[0::100]
    decSol0 = sol[:,0][0::100]
    decSol1 = sol[:,1][0::100]
    decSol2 = sol[:,2][0::100]
    decSol3 = sol[:,3][0::100]

    p1=list_plot(list(zip(dect,decSol0/10)), legend_label="Virus (1/10th)", plotjoined=True)
    p2=list_plot(list(zip(dect,decSol1)), color="red", legend_label="Uninfected T-cells", plotjoined=True)
    p3=list_plot(list(zip(dect,decSol2)), color="purple", legend_label="Latently infected T-cells", plotjoined=True)
    p4=list_plot(list(zip(dect,decSol3)), color="orange", legend_label="Actively infected T-cells", plotjoined=True)
    if showV==True:
        show(p1+p2+p3+p4, show_legend=True)
    else:
        show(p2+p3+p4, show_legend=True)