chap18


Modeling and Simulation in Python

Chapter 18

Copyright 2017 Allen Downey

License: Creative Commons Attribution 4.0 International


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 *

Code from the previous chapter

Read the data.


In [2]:
data = pd.read_csv('data/glucose_insulin.csv', index_col='time');

Interpolate the insulin data.


In [3]:
I = interpolate(data.insulin)

The glucose minimal model

I'll cheat by starting with parameters that fit the data roughly; then we'll see how to improve them.


In [4]:
params = Params(G0 = 290,
                k1 = 0.03,
                k2 = 0.02,
                k3 = 1e-05)

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)

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)

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

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')

Numerical solution

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)

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

results is a TimeFrame with one row for each time step and one column for each state variable:


In [17]:
results2

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

In [20]:
max(abs(percent_diff))

Exercises

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 goes here

In [22]:
# Solution goes here

In [23]:
# Solution goes here

In [24]:
# Solution goes here

Under the hood

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 [ ]: