code/chap18.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 18\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code from the previous chapter\n",
    "\n",
    "Read the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('data/glucose_insulin.csv', index_col='time');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolate the insulin data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "I = interpolate(data.insulin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "G0 = 290\n",
    "k1 = 0.03\n",
    "k2 = 0.02\n",
    "k3 = 1e-05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate basal levels, we'll use the concentrations at `t=0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Gb = data.glucose[0]\n",
    "Ib = data.insulin[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the initial condtions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "init = State(G=G0, X=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_0 = get_first_label(data)\n",
    "t_end = get_last_label(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "system = System(init=init, \n",
    "                k1=k1, k2=k2, k3=k3,\n",
    "                I=I, Gb=Gb, Ib=Ib,\n",
    "                t_0=t_0, t_end=t_end, dt=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Updates the glucose minimal model.\n",
    "    \n",
    "    state: State object\n",
    "    t: time in min\n",
    "    system: System object\n",
    "    \n",
    "    returns: State object\n",
    "    \"\"\"\n",
    "    G, X = state\n",
    "    unpack(system)\n",
    "        \n",
    "    dGdt = -k1 * (G - Gb) - X*G\n",
    "    dXdt = k3 * (I(t) - Ib) - k2 * X\n",
    "    \n",
    "    G += dGdt * dt\n",
    "    X += dXdt * dt\n",
    "\n",
    "    return State(G=G, X=X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "        \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    unpack(system)\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t_0] = init\n",
    "    ts = linrange(t_0, t_end, dt)\n",
    "    \n",
    "    for t in ts:\n",
    "        frame.row[t+dt] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time results = run_simulation(system, update_func);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical solution\n",
    "\n",
    "In the previous chapter, we approximated the differential equations with difference equations, and solved them using `run_simulation`.\n",
    "\n",
    "In this chapter, we solve the differential equation numerically using `run_ode_solver`, which is a wrapper for the SciPy ODE solver.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes derivatives of the glucose minimal model.\n",
    "    \n",
    "    state: State object\n",
    "    t: time in min\n",
    "    system: System object\n",
    "    \n",
    "    returns: derivatives of G and X\n",
    "    \"\"\"\n",
    "    G, X = state\n",
    "    unpack(system)\n",
    "    \n",
    "    dGdt = -k1 * (G - Gb) - X*G\n",
    "    dXdt = k3 * (I(t) - Ib) - k2 * X\n",
    "    \n",
    "    return dGdt, dXdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_func(init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we run the ODE solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time results2, details = run_ode_solver(system, slope_func, t_eval=data.index);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`details` is a `ModSimSeries` object with information about how the solver worked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`results` is a `TimeFrame` with one row for each time step and one column for each state variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "results2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the results from `run_simulation` and `run_ode_solver`, we can see that they are not very different."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(results.G, 'g-')\n",
    "plot(results2.G, 'b-')\n",
    "plot(data.glucose, 'bo')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The differences in `G` are less than 1%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = results.G - results2.G\n",
    "percent_diff = diff / results2.G * 100\n",
    "percent_diff.dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's find the parameters that yield the best fit for the data.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use these values as an initial estimate and iteratively improve them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = Params(G0 = 290,\n",
    "                k1 = 0.03,\n",
    "                k2 = 0.02,\n",
    "                k3 = 1e-05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`make_system` takes the parameters and actual data and returns a `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_system(params, data):\n",
    "    \"\"\"Makes a System object with the given parameters.\n",
    "    \n",
    "    params: sequence of G0, k1, k2, k3\n",
    "    data: DataFrame with `glucose` and `insulin`\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    G0, k1, k2, k3 = params\n",
    "    \n",
    "    Gb = data.glucose[0]\n",
    "    Ib = data.insulin[0]\n",
    "    \n",
    "    t_0 = get_first_label(data)\n",
    "    t_end = get_last_label(data)\n",
    "\n",
    "    init = State(G=G0, X=0)\n",
    "    \n",
    "    return System(G0=G0, k1=k1, k2=k2, k3=k3,\n",
    "                  init=init, Gb=Gb, Ib=Ib,\n",
    "                  t_0=t_0, t_end=t_end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def error_func(params, data):\n",
    "    \"\"\"Computes an array of errors to be minimized.\n",
    "    \n",
    "    params: sequence of parameters\n",
    "    data: DataFrame of values to be matched\n",
    "    \n",
    "    returns: array of errors\n",
    "    \"\"\"\n",
    "    print(params)\n",
    "    \n",
    "    # make a System with the given parameters\n",
    "    system = make_system(params, data)\n",
    "    \n",
    "    # solve the ODE\n",
    "    results, details = run_ode_solver(system, slope_func, t_eval=data.index)\n",
    "    \n",
    "    # compute the difference between the model\n",
    "    # results and actual data\n",
    "    errors = results.G - data.glucose\n",
    "    return errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we call `error_func`, we provide a sequence of parameters as a single object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how that works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_func(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`fit_leastsq` is a wrapper for `scipy.optimize.leastsq`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we call it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params, fit_details = fit_leastsq(error_func, params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first return value is a `Params` object with the best parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second return value is a `ModSimSeries` object with information about the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have `best_params`, we can use it to make a `System` object and run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "system = make_system(best_params, data)\n",
    "results, details = run_ode_solver(system, slope_func, t_eval=data.index)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(results.G, label='simulation')\n",
    "plot(data.glucose, 'bo', label='glucose data')\n",
    "\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='Concentration (mg/dL)')\n",
    "\n",
    "savefig('figs/chap08-fig04.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting parameters\n",
    "\n",
    "Based on the parameters of the model, we can estimate glucose effectiveness and insulin sensitivity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def indices(params):\n",
    "    \"\"\"Compute glucose effectiveness and insulin sensitivity.\n",
    "    \n",
    "    params: sequence of G0, k1, k2, k3\n",
    "    data: DataFrame with `glucose` and `insulin`\n",
    "    \n",
    "    returns: State object containing S_G and S_I\n",
    "    \"\"\"\n",
    "    G0, k1, k2, k3 = params\n",
    "    return State(S_G=k1, S_I=k3/k2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "indices(best_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Under the hood\n",
    "\n",
    "Here's the source code for `run_ode_solver` and `fit_leastsq`, if you'd like to know how they work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "%psource run_ode_solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%psource fit_leastsq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "**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.\n",
    "\n",
    "Does that improve the quality of the fit?  Does it change the best parameters by much?\n",
    "\n",
    "Note: You can read more about this use of `loc` [in the Pandas documentation](https://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-integer)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**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?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Related reading:** You might be interested in this article about [people making a DIY artificial pancreas](https://www.bloomberg.com/news/features/2018-08-08/the-250-biohack-that-s-revolutionizing-life-with-diabetes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}