Model.ipynb

Kernel: Python 3

# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *

def make_system(gamma,mu,beta,tau,roe,alpha,pi,sigma,delta):
    """Make a system object for the SIR model.
    
    beta: contact rate in days
    gamma: recovery rate in days
    
    returns: System object
    """
    init = State(R=200, L=0, E=0, V=4e-7)

    t_0 = 0
    t_end = 120

    return System(init=init, t_0=t_0, t_end=t_end,
                  beta=beta, gamma=gamma,mu=mu,tau=tau,roe=roe,alpha=alpha,sigma=sigma, pi=pi, delta = delta, dt=.0541)

system = make_system(1.36,0.00136,.00027,0.2,0.1,0.036,100,2,.33)

	values
init	R 2.000000e+02 L 0.000000e+00 E 0.000...
t_0	0
t_end	120
beta	0.00027
gamma	1.36
mu	0.00136
tau	0.2
roe	0.1
alpha	0.036
sigma	2
pi	100
delta	0.33
dt	0.0541

def update_func(state, t, system):
    R, L, E, V = state
    diffR = R
    diffL = L
    diffE = E
    diffV = V
    diffR = system.gamma*system.tau-system.mu*R-system.beta*R*V
a    diffE = (1-system.roe)*system.beta* R * V + system.alpha * L - system.delta * E
    diffV = system.pi * E - system.sigma * V
    R += diffR * system.dt
    L += diffL * system.dt
    E += diffE * system.dt
    V += diffV * system.dt
    return State(R=R, L=L, E=E, V=V)

def run_simulation(system, update_func):
    """Runs a simulation of the system.
        
    system: System object
    update_func: function that updates state
    
    returns: TimeFrame
    """
    init = system.init
    t_0, t_end, dt = system.t_0, system.t_end, system.dt
    
    frame = TimeFrame(columns=init.index)
    frame.row[t_0] = init
    ts = linrange(t_0, t_end, system.dt)
    
    for t in ts:
        frame.row[t+dt] = update_func(frame.row[t], t, system)
    
    return frame

results = run_simulation(system, update_func)

	R	L	E	V
0.0000	200.000000	0.000000e+00	0.000000e+00	4.000000e-07
0.0541	200.000000	1.168560e-10	1.051704e-09	3.567200e-07
0.1082	200.000000	2.208320e-10	1.971065e-09	3.238126e-07
0.1623	200.000000	3.149843e-10	2.787693e-09	2.994396e-07
0.2164	200.000000	4.018259e-10	3.525843e-09	2.821216e-07
...	...	...	...	...
119.7774	18.008168	6.612353e-01	2.387806e-01	1.200246e+01
119.8315	18.018401	6.602145e-01	2.386470e-01	1.199560e+01
119.8856	18.028633	6.591958e-01	2.385137e-01	1.198875e+01
119.9397	18.038865	6.581792e-01	2.383808e-01	1.198193e+01
119.9938	18.049096	6.571646e-01	2.382483e-01	1.197512e+01

2219 rows × 4 columns

plot(results.V)
decorate(title = 'Free Virions',
         xlabel = 'Days',
         ylabel = 'Virions')

plot(results.R)
decorate(title = 'Uninfected CD4',
         xlabel = 'Days',
         ylabel = 'CD4 Lymphosytes')

plot(results.E)
decorate(title = 'Infected CD4',
         xlabel = 'Days',
         ylabel = 'Infected CD4',
         yscale = 'log')

plot(results.L)

decorate(title = 'Free Virions',
         xlabel = 'Days',
         ylabel = 'Latently Infected CD4')