In [1]:
# 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 *
In [2]:
data = pd.read_csv('data/glucose_insulin.csv', index_col='time');
Interpolate the insulin data.
In [3]:
I = interpolate(data.insulin)
Out[3]:
In [4]:
params = Params(G0 = 290,
k1 = 0.03,
k2 = 0.02,
k3 = 1e-05)
Out[4]:
Here's a version of make_system that takes the parameters and data:
In [5]:
def make_system(params, data):
"""Makes a System object with the given parameters.
params: sequence of G0, k1, k2, k3
data: DataFrame with `glucose` and `insulin`
returns: System object
"""
G0, k1, k2, k3 = params
Gb = data.glucose[0]
Ib = data.insulin[0]
I = interpolate(data.insulin)
t_0 = get_first_label(data)
t_end = get_last_label(data)
init = State(G=G0, X=0)
return System(params,
init=init, Gb=Gb, Ib=Ib, I=I,
t_0=t_0, t_end=t_end, dt=2)
In [6]:
system = make_system(params, data)
Out[6]:
And here's the update function.
In [7]:
def update_func(state, t, system):
"""Updates the glucose minimal model.
state: State object
t: time in min
system: System object
returns: State object
"""
G, X = state
k1, k2, k3 = system.k1, system.k2, system.k3
I, Ib, Gb = system.I, system.Ib, system.Gb
dt = system.dt
dGdt = -k1 * (G - Gb) - X*G
dXdt = k3 * (I(t) - Ib) - k2 * X
G += dGdt * dt
X += dXdt * dt
return State(G=G, X=X)
Before running the simulation, it is always a good idea to test the update function using the initial conditions. In this case we can veryify that the results are at least qualitatively correct.
In [8]:
update_func(system.init, system.t_0, system)
Out[8]:
Now run_simulation is pretty much the same as it always is.
In [9]:
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, dt)
for t in ts:
frame.row[t+dt] = update_func(frame.row[t], t, system)
return frame
And here's how we run it.
In [10]:
results = run_simulation(system, update_func);
The results are in a TimeFrame object with one column per state variable.
In [11]:
results
Out[11]:
The following plot shows the results of the simulation along with the actual glucose data.
In [12]:
subplot(2, 1, 1)
plot(results.G, 'b-', label='simulation')
plot(data.glucose, 'bo', label='glucose data')
decorate(ylabel='Concentration (mg/dL)')
subplot(2, 1, 2)
plot(results.X, 'C1', label='remote insulin')
decorate(xlabel='Time (min)',
ylabel='Concentration (arbitrary units)')
savefig('figs/chap18-fig01.pdf')
Now let's solve the differential equation numerically using run_ode_solver, which is an implementation of Ralston's method.
Instead of an update function, we provide a slope function that evaluates the right-hand side of the differential equations.
We don't have to do the update part; the solver does it for us.
In [13]:
def slope_func(state, t, system):
"""Computes derivatives of the glucose minimal model.
state: State object
t: time in min
system: System object
returns: derivatives of G and X
"""
G, X = state
k1, k2, k3 = system.k1, system.k2, system.k3
I, Ib, Gb = system.I, system.Ib, system.Gb
dGdt = -k1 * (G - Gb) - X*G
dXdt = k3 * (I(t) - Ib) - k2 * X
return dGdt, dXdt
We can test the slope function with the initial conditions.
In [14]:
slope_func(system.init, 0, system)
Out[14]:
Here's how we run the ODE solver.
In [15]:
results2, details = run_ode_solver(system, slope_func)
details is a ModSimSeries object with information about how the solver worked.
In [16]:
details
Out[16]:
results is a TimeFrame with one row for each time step and one column for each state variable:
In [17]:
results2
Out[17]:
Plotting the results from run_simulation and run_ode_solver, we can see that they are not very different.
In [18]:
plot(results.G, 'C0', label='run_simulation')
plot(results2.G, 'C2--', label='run_ode_solver')
decorate(xlabel='Time (min)', ylabel='Concentration (mg/dL)')
savefig('figs/chap18-fig02.pdf')
The differences in G are less than 2%.
In [19]:
diff = results.G - results2.G
percent_diff = diff / results2.G * 100
percent_diff
Out[19]:
In [20]:
max(abs(percent_diff))
Out[20]:
Exercise: Our solution to the differential equations is only approximate because we used a finite step size, dt=2 minutes.
If we make the step size smaller, we expect the solution to be more accurate. Run the simulation with dt=1 and compare the results. What is the largest relative error between the two solutions?
In [21]:
# Solution
system3 = System(system, dt=1)
results3, details = run_ode_solver(system3, slope_func)
details
Out[21]:
In [22]:
# Solution
plot(results2.G, 'C2--', label='run_ode_solver (dt=2)')
plot(results3.G, 'C3:', label='run_ode_solver (dt=1)')
decorate(xlabel='Time (m)', ylabel='mg/dL')
In [23]:
# Solution
diff = (results2.G - results3.G).dropna()
percent_diff = diff / results2.G * 100
Out[23]:
In [24]:
# Solution
max(abs(percent_diff))
Out[24]:
Here's the source code for run_ode_solver if you'd like to know how it works.
Notice that run_ode_solver is another name for run_ralston, which implements Ralston's method.
In [25]:
source_code(run_ode_solver)
Related reading: You might be interested in this article about people making a DIY artificial pancreas.
In [ ]: