<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 The Insulin Minimal Model
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License
Step1:  The following cells download and read the data.
Step2:  In Chapter 17 I present the glucose minimal model; in Chapter 18 we implemented it using run_simulation and run_solve_ivp.
In this case study, we'll implement the other half of the Minimal Model, which describes the concentration of insulin.
In the insulin minimal model, the concentration of insulin, $I$, is governed by this differential equation
Step3:  Exercise
Step4:  Exercise
Step6:  Exercise
Step7:  Exercise
Step8:  Exercise
Step9:  According to Pacini and Bergman, here are the normal ranges for these quantities.
<ASSISTANT_TASK:>
Python Code:
# install Pint if necessary
try:
    import pint
except ImportError:
    !pip install pint
# download modsim.py if necessary
from os.path import basename, exists
def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)
    
download('https://raw.githubusercontent.com/AllenDowney/' +
         'ModSim/main/modsim.py')
# import functions from modsim
from modsim import *
Explanation: The Insulin Minimal Model
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International
End of explanation
download('https://raw.githubusercontent.com/AllenDowney/' +
         'ModSim/main/data/glucose_insulin.csv')
data = pd.read_csv('glucose_insulin.csv', index_col='time');
Explanation: The following cells download and read the data.
End of explanation
I0 = 360   
k = 0.25
gamma = 0.004
G_T = 80
params = I0, k, gamma, G_T
# Solution
def make_system(params, data):
    I0, k, gamma, G_T = params
    
    init = State(I=I0)
    
    t_0 = data.index[0]
    t_end = data.index[-1]
    G=interpolate(data.glucose)
    
    system = System(I0=I0, k=k, gamma=gamma, G_T=G_T, G=G,
                    init=init, t_0=t_0, t_end=t_end, dt=1)
    return system
# Solution
system = make_system(params, data)
Explanation: In Chapter 17 I present the glucose minimal model; in Chapter 18 we implemented it using run_simulation and run_solve_ivp.
In this case study, we'll implement the other half of the Minimal Model, which describes the concentration of insulin.
In the insulin minimal model, the concentration of insulin, $I$, is governed by this differential equation:
$ \frac{dI}{dt} = -k I(t) + \gamma (G(t) - G_T) t $
where $G(t)$ is the concentration of glucose at time $t$, and
$k$, $\gamma$, and $G_T$ are positive-valued parameters:
$k$ controls the rate of insulin disappearance.
$G_T$ determines the glucose threshold: when $G(t)$ exceeds this threshold, it causes insulin to appear; when $G(t)$ is below this threshold, it causes insulin to disappear.
$\gamma$ controls how quickly insulin appears or disappears when the concentration of glucose is elevated or depressed.
Notice that this equation depends on time, $t$, since the initial injection.  It has the effect of increasing glucose sensitivity over time.  If you are familiar with control systems, the effect of this term is similar to the integral term in a PID controller.
In addition to the three parameters in the equation, we will also consider the initial concentration of insulin, $I_0$, to be a free parameter; that is, we will choose the value of $I_0$, and the other parameters, that best fit the data.
Exercise:  Write a version of make_system that takes the parameters of the 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_solve_ivp.
Use it to make a System object with the following parameters:
End of explanation
# Solution
def slope_func(t, state, system):
    I = state
    k, gamma = system.k, system.gamma
    G, G_T = system.G, system.G_T
        
    dIdt = -k * I + gamma * (G(t) - G_T) * t 
    return dIdt
# Solution
slope_func(system.t_0, system.init, system)
Explanation: Exercise: Write a slope function that takes a time stamp, a State object, and a System object, and returns the derivative of I with respect to time.  Test your function with the initial conditions from system.
End of explanation
# Solution
results, details = run_solve_ivp(system, slope_func,
                                 t_eval=data.index)
details.message
# Solution
results.tail()
# Solution
data.insulin.plot(style='o', color='C2',
                  alpha=0.5, label='insulin data')
results.I.plot(style='-', color='C2', label='simulation')
decorate(xlabel='Time (min)',
         ylabel='Concentration ($\mu$U/mL)')
Explanation: Exercise: Run run_solve_ivp with your System object and slope function, and plot the results, along with the measured insulin levels.  Use the keyword argument t_eval=data.index so the results are evaluated as the same time stamps as the data.
End of explanation
# 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, details = run_solve_ivp(system, slope_func,
                                     t_eval=data.index)
    # compute the difference between the model
    # results and actual data
    errors = results.I - data.insulin
    return errors.iloc[2:]
# Solution
error_func(params, data)
Explanation: 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, call run_solve_ivp, 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 with the glucose model, you might want to drop the first 2-3 elements from the sequence of errors.
End of explanation
# Solution
best_params, details = leastsq(error_func, params, data)
# Solution
details.mesg
# Solution
system2 = make_system(best_params, data)
# Solution
results2, details = run_solve_ivp(system2, slope_func, 
                                  t_eval=data.index)
details.message
# Solution
data.insulin.plot(style='o', color='C2',
                  alpha=0.5, label='insulin data')
results2.I.plot(style='-', color='C2', label='simulation')
decorate(xlabel='Time (min)',
         ylabel='Concentration ($\mu$U/mL)')
Explanation: Exercise: Use 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.
End of explanation
# Solution
I0, k, gamma, G_T = best_params
# 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 = 272
Gb = data.glucose[0]
G0, Gb
# Solution
phi_1 = (I_max - Ib) / k / (G0 - Gb)
phi_1
# Solution
phi_2 = gamma * 1e4
phi_2
Explanation: 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 $
For $G_0$, use the best estimate from the glucose model, 272 mg/dL.  For $G_b$ and $I_b$, use the initial measurements from the data.
For $I_{max}$ is the maximum measurement of insulin concentration.
End of explanation
phi_1_interval = 2, 4
phi_1_interval
phi_2_interval = 20, 35
phi_2_interval
Explanation: According to Pacini and Bergman, here are the normal ranges for these quantities.
End of explanation
<END_TASK>