Four-variable HIV model.sagews

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}$

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,
)
vector_field(V, R, L, E) = system
initial_state = (4E-7, 200, 0, 0)

Error in lines 12-12 Traceback (most recent call last): File "/cocalc/lib/python3.9/site-packages/smc_sagews/sage_server.py", line 1230, in execute exec( File "", line 1, in <module> File "/ext/sage/9.4/local/lib/python3.9/site-packages/sage/calculus/desolvers.py", line 1713, in desolve_odeint all_vars = set().union(*[de.variables() for de in des]) File "/ext/sage/9.4/local/lib/python3.9/site-packages/sage/calculus/desolvers.py", line 1713, in <listcomp> all_vars = set().union(*[de.variables() for de in des]) File "sage/structure/element.pyx", line 493, in sage.structure.element.Element.__getattr__ (build/cythonized/sage/structure/element.c:4708) return self.getattr_from_category(name) File "sage/structure/element.pyx", line 506, in sage.structure.element.Element.getattr_from_category (build/cythonized/sage/structure/element.c:4820) return getattr_from_other_class(self, cls, name) File "sage/cpython/getattr.pyx", line 367, in sage.cpython.getattr.getattr_from_other_class (build/cythonized/sage/cpython/getattr.c:2551) raise AttributeError(dummy_error_message) AttributeError: 'FreeModule_ambient_field_with_category.element_class' object has no attribute 'variables'

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.

%auto
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, default=0.00027, label="$\\beta$"),
                    alpha=slider(0.00036, 0.036, 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,
    )
    vector_field(V, R, L, E) = system

    # 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(vector_field, 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)"))

Interact: please open in CoCalc