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]:
Initialize the parameters
In [4]:
G0 = 290
k1 = 0.03
k2 = 0.02
k3 = 1e-05
Out[4]:
To estimate basal levels, we'll use the concentrations at t=0.
In [5]:
Gb = data.glucose[0]
Ib = data.insulin[0]
Out[5]:
Create the initial condtions.
In [6]:
init = State(G=G0, X=0)
Out[6]:
Make the System object.
In [7]:
t_0 = get_first_label(data)
t_end = get_last_label(data)
Out[7]:
In [8]:
system = System(G0=G0, k1=k1, k2=k2, k3=k3,
init=init, Gb=Gb, Ib=Ib, I=I,
t_0=t_0, t_end=t_end, dt=2)
Out[8]:
In [9]:
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)
In [10]:
def run_simulation(system, update_func):
"""Runs a simulation of the system.
system: System object
update_func: function that updates state
returns: TimeFrame
"""
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
In [11]:
results = run_simulation(system, update_func);
In the previous chapter, we approximated the differential equations with difference equations, and solved them using run_simulation.
In this chapter, we solve the differential equation numerically using run_ode_solver, which is a wrapper for the SciPy ODE solver.
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 [12]:
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 [13]:
slope_func(init, 0, system)
Out[13]:
Here's how we run the ODE solver.
In [14]:
results2, details = run_ode_solver(system, slope_func, t_eval=data.index);
details is a ModSimSeries object with information about how the solver worked.
In [15]:
details
Out[15]:
results is a TimeFrame with one row for each time step and one column for each state variable:
In [16]:
results2
Out[16]:
Plotting the results from run_simulation and run_ode_solver, we can see that they are not very different.
In [17]:
plot(results.G, '-')
plot(results2.G, '-')
plot(data.glucose, 'bo')
Out[17]:
The differences in G are less than 1%.
In [18]:
diff = results.G - results2.G
percent_diff = diff / results2.G * 100
percent_diff.dropna()
Out[18]:
Now let's find the parameters that yield the best fit for the data.
We'll use these values as an initial estimate and iteratively improve them.
In [19]:
params = Params(G0 = 290,
k1 = 0.03,
k2 = 0.02,
k3 = 1e-05)
Out[19]:
make_system takes the parameters and actual data and returns a System object.
In [20]:
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
"""
# params might be a Params object or an array,
# so we have to unpack it like this
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(G0=G0, k1=k1, k2=k2, k3=k3,
init=init, Gb=Gb, Ib=Ib, I=I,
t_0=t_0, t_end=t_end, dt=2)
In [21]:
system = make_system(params, data)
Out[21]:
error_func takes the parameters and actual data, makes a System object, and runs odeint, then compares the results to the data. It returns an array of errors.
In [22]:
system = make_system(params, data)
results, details = run_ode_solver(system, slope_func, t_eval=data.index)
details
Out[22]:
In [23]:
def error_func(params, data):
"""Computes an array of errors to be minimized.
params: sequence of parameters
data: DataFrame of values to be matched
returns: array of errors
"""
print(params)
# make a System with the given parameters
system = make_system(params, data)
# solve the ODE
results, details = run_ode_solver(system, slope_func, t_eval=data.index)
# compute the difference between the model
# results and actual data
errors = results.G - data.glucose
return errors
When we call error_func, we provide a sequence of parameters as a single object.
Here's how that works:
In [24]:
error_func(params, data)
Out[24]:
leastsq is a wrapper for scipy.optimize.leastsq
Here's how we call it.
In [25]:
best_params, fit_details = leastsq(error_func, params, data)
The first return value is a Params object with the best parameters:
In [26]:
best_params
Out[26]:
The second return value is a ModSimSeries object with information about the results.
In [27]:
fit_details
Out[27]:
Now that we have best_params, we can use it to make a System object and run it.
In [28]:
system = make_system(best_params, data)
results, details = run_ode_solver(system, slope_func, t_eval=data.index)
details.message
Out[28]:
Here are the results, along with the data. The first few points of the model don't fit the data, but we don't expect them to.
In [29]:
plot(results.G, label='simulation')
plot(data.glucose, 'bo', label='glucose data')
decorate(xlabel='Time (min)',
ylabel='Concentration (mg/dL)')
savefig('figs/chap08-fig04.pdf')
In [30]:
def indices(params):
"""Compute glucose effectiveness and insulin sensitivity.
params: sequence of G0, k1, k2, k3
data: DataFrame with `glucose` and `insulin`
returns: State object containing S_G and S_I
"""
G0, k1, k2, k3 = params
return State(S_G=k1, S_I=k3/k2)
Here are the results.
In [31]:
indices(best_params)
Out[31]:
In [32]:
source_code(run_ode_solver)
In [33]:
source_code(leastsq)
Exercise: Since we don't expect the first few points to agree, it's probably better not to make them part of the optimization process. We can ignore them by leaving them out of the Series returned by error_func. Modify the last line of error_func to return errors.loc[8:], which includes only the elements of the Series from t=8 and up.
Does that improve the quality of the fit? Does it change the best parameters by much?
Note: You can read more about this use of loc in the Pandas documentation.
Exercise: How sensitive are the results to the starting guess for the parameters? If you try different values for the starting guess, do we get the same values for the best parameters?
Related reading: You might be interested in this article about people making a DIY artificial pancreas.
In [ ]: