<jupyter_start><jupyter_text>Modeling and Simulation in PythonInsulin minimal modelCopyright 2017 <NAME>License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)<jupyter_code># 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 *<jupyter_output><empty_output><jupyter_text>DataWe have data from Pacini and Bergman (1986), "MINMOD: a computer program to calculate insulin sensitivity and pancreatic responsivity from the frequently sampled intravenous glucose tolerance test", *Computer Methods and Programs in Biomedicine*, 23: 113-122..<jupyter_code>data = pd.read_csv('data/glucose_insulin.csv', index_col='time');<jupyter_output><empty_output><jupyter_text>The return value from `interpolate` is a function.  The insulin minimal modelIn addition to the glucose minimal mode, Pacini and Bergman present an insulin minimal model, in which the concentration of insulin, $I$, is governed by this differential equation:$ \frac{dI}{dt} = -k I(t) + \gamma (G(t) - G_T) t $ **Exercise:**  Write a version of `make_system` that takes the parameters of this model, `I0`, `k`, `gamma`, and `G_T` as parameters, along with a `DataFrame` containing the measurements, and returns a `System` object suitable for use with `run_simulation` or `run_odeint`.Use it to make a `System` object with the following parameters:<jupyter_code>I0 = 360
k = 0.25
gamma = 0.004
G_T = 80
# Solution

def make_system(I0, k, gamma, G_T, data):
    init = State(I=I0)
    system = System(init=init, 
                    k=k, gamma=gamma, G_T=G_T,
                    G=interpolate(data.glucose),
                    ts=data.index)
    return system
# Solution

system = make_system(I0, k, gamma, G_T, data)<jupyter_output><empty_output><jupyter_text>**Exercise:** Write a slope function that takes state, t, system as parameters and returns the derivative of `I` with respect to time.  Test your function with the initial condition $I(0)=360$.<jupyter_code># Solution

def slope_func(state, t, system):
    I = state
    unpack(system)
    
    dIdt = - k * I + gamma * (G(t) - G_T) * t
    
    return dIdt
# Solution

init = State(I=I0)
slope_func(init, 0, system)<jupyter_output><empty_output><jupyter_text>**Exercise:** Run `run_odeint` with your `System` object and slope function, and plot the results, along with the measured insulin levels.<jupyter_code># Solution

results = run_odeint(system, slope_func)
# Solution

plot(results, 'g-', label='simulation')
plot(data.insulin, 'go', label='insulin data')

decorate(xlabel='Time (min)',
         ylabel='Concentration ($\mu$U/mL)')<jupyter_output><empty_output><jupyter_text>**Exercise:**  Write an error function that takes a sequence of parameters as an argument, along with the `DataFrame` containing the measurements.  It should make a `System` object with the given parameters, run it, and compute the difference between the results of the simulation and the measured values.  Test your error function by calling it with the parameters from the previous exercise.Hint: As we did in a previous exercise, you might want to drop the errors for times prior to `t=8`.<jupyter_code># Solution

def error_func(params, data):
    """Computes an array of errors to be minimized.
    
    params: sequence of parameters
    actual: array 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 = run_odeint(system, slope_func)
    
    # compute the difference between the model
    # results and actual data
    errors = results.I - data.insulin
    return errors.loc[8:]
# Solution

params = [I0, k, gamma, G_T]
params
# Solution

error_func(params, data)<jupyter_output>[360, 0.25, 0.004, 80]<jupyter_text>**Exercise:** Use `fit_leastsq` to find the parameters that best fit the data.  Make a `System` object with those parameters, run it, and plot the results along with the measurements.<jupyter_code># Solution

best_params = fit_leastsq(error_func, params, data)
# Solution

system = make_system(*best_params, data)
# Solution

results = run_odeint(system, slope_func)
# Solution

plot(results, 'g-', label='simulation')
plot(data.insulin, 'go', label='insulin data')

decorate(xlabel='Time (min)',
         ylabel='Concentration ($\mu$U/mL)')<jupyter_output><empty_output><jupyter_text>**Exercise:** Using the best parameters, estimate the sensitivity to glucose of the first and second phase pancreatic responsivity:$ \phi_1 = \frac{I_{max} - I_b}{k (G_0 - G_b)} $$ \phi_2 = \gamma \times 10^4 $<jupyter_code># Solution

I_max = data.insulin.max()
Ib = data.insulin[0]
I_max, Ib
# Solution

# The value of G0 is the best estimate from the glucose model
G0 = 289
Gb = data.glucose[0]
G0, Gb
# Solution

phi_1 = (I_max - Ib) / k / (G0 - Gb)
phi_1
# Solution

phi_2 = gamma * 1e4
phi_2<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.102321083172147, 'fraction_numerical': 0.01528046421663443, 'mean_word_length': 4.596320346320346, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '35265', 'license': 'MIT', 'n_tokens_mistral': 1715, 'n_tokens_neox': 1600, 'n_words': 565, 'path': 'code/soln/insulin_soln.ipynb', 'repo_name': 'rsmxingu/ModSimPy', 'size': '5176'}