import numpy
state_vars = list(var("V, R, E"))
system = (
100*E - 2*V,
0.272 - 0.00136*R - 0.00027*R*V,
0.00027*R*V-0.33*E,
)
@interact
def HIV_interactive(initV=slider(10^-7, 10^-6, 10^-7, default=4E-7, label="Initial $V$"),
initR=slider(100, 1000, 10, default=200, label="Initial $R$"),
initE=slider(0, 100, 1, default=0, label="Initial $E$"),
tmax=slider(10, 3660, 10, default=400, label="Simulation length (days)")):
initial_state = (initV, initR, initE)
t_range = srange(0, tmax, 0.025)
solution = desolve_odeint(system, initial_state, t_range, state_vars)
solution = numpy.insert(solution, 0, t_range, axis=1)
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$)")
E_plot = list_plot(solution[::100,(0,3)], plotjoined=True,
color="purple", legend_label="Infected cells ($E$)")
show(V_plot + R_plot + E_plot, ymax=550,
axes_labels=("$t$", "Populations\n(per mm$^3$ of blood)"))