adding ODE exercise

docs/requirements.txt

@@ -4,6 +4,7 @@ robpol86-sphinxcontrib-googleanalytics
 matplotlib
 numpy
 optlang
+scipy
 sphinx
 sphinx-rtd-theme
 sqlalchemy

docs/tutorials/cell_modeling/simulation/index.rst

@@ -28,7 +28,15 @@ There are also several generalizations of Boolean and logical-valued models incl
 
 Ordinary differential equations (ODEs)
 --------------------------------------
-Ordinary differential equations (ODEs) is one of the most commonly used approaches for modeling dynamical systems. ODE models are based on microscopic analyses of how the concentration of each species in the cell changes over time in response to the concentrations of other species. Because ODE assume that cells are well-mixed and that they behave deterministically, ODE models are most appropriate for small systems that involve large concentrations and high fluxes. ODE models can be simulated by numerically integrating the differential equations.  
+Ordinary differential equations (ODEs) is one of the most commonly used approaches for modeling dynamical systems. ODE models are based on microscopic analyses of how the concentration of each species in the cell changes over time in response to the concentrations of other species. Because ODE assume that cells are well-mixed and that they behave deterministically, ODE models are most appropriate for small systems that involve large concentrations and high fluxes. 
+
+ODE models can be simulated by numerically integrating the differential equations. The most basic ODE integration method is Euler's method. Euler's method is a time stepped algorithm in which the next state is computed by adding the current state and the multiplication of the current differentials with the timestep.
+
+.. math::
+
+    y(t+\Delta t) = y(t) + \Delta t \frac{dy}{dt}
+
+Euler's method estimates :math:`y(t+\Delta t)` using a first-order approximation. ODE models can be simulated more accurately using higher order estimates, or Taylor series, of :math:`y(t+\Delta t)`. One of the most popular algorithms which implements this approach is the Runge-Kutta 4th order method. Yet more advanced integration methods select the time step adaptively. Some of the most sophisticated ODE integration packages include ODEPACK (lsoda) and Sundials (vode). These packages can be used in Python via scipy's integrate module.
 
 
 Stochastic simulation
@@ -98,7 +106,7 @@ dFBA enables dynamic simulations by (1) assuming that cells quickly reach pseudo
         Solve for the maximum growth rate and optimal fluxes
         Update the biomass concentration based on the predicted growth rate
         Update the environmental conditions based on the predicted exchange fluxes
-        
+
 
 Multi-algorithmic simulation
 ----------------------------
@@ -179,7 +187,8 @@ For Ubuntu, the following commands can be used to install all of the software re
     sudo pip install \
         matplotlib \
         numpy \
-        optlang 
+        optlang \
+        scipy
 
 
 Boolean simulation
@@ -316,13 +325,46 @@ Finally, compare the simulation results from the different update schemes. How d
 
 ODE simulation
 ^^^^^^^^^^^^^^
+In this exercise we will use lsoda to simulate the Goldbeter 1991 cell cycle model of cyclin, cdc2, and cyclin protease (`doi: 10.1073/pnas.88.20.9107 <https://doi.org/10.1073/pnas.88.20.9107>`_, `BIOMD0000000003 <http://www.ebi.ac.uk/biomodels-main/BIOMD0000000003>`_).
+
+.. image:: ode-model.png
+
+First, open the `BioModels entry <http://www.ebi.ac.uk/biomodels-main/BIOMD0000000003>`_ for the model in your web browser. Identify the reactions, their rate laws, the parameter values, and the initial conditions of the model. Note, that the model uses two assignment rules for :math:`V1` and :math:`V3` which are not displayed on the BioModels page. These assignment rules must be identified from the SBML version of the model which can be exported from the BioModels page.
+
+Second, write a function to calculate the time derivative of the cyclin, cd2, and protease concentrations/activities by completing the code fragment below::
+
+    def d_conc_d_t(concs, time):
+        """ Calculate differentials for Goldbeter 1991 cell cycle model 
+        (`BIOMD0000000003 <http://www.ebi.ac.uk/biomodels-main/BIOMD0000000003>`_)
+
+        Args:
+            time (obj:`float`): time
+            concs (:obj:`numpy.ndarray`): array of current concentrations
+
+        Returns:
+            :obj:`numpy.ndarray
+        """
+        ...
+
+Next, create a vector to represent the initial conditions by completing this code fragment::
+
+    init_concs = ...
+
+Next, use ``scipy.integrate.odeint`` which implements the lsoda algorithm to simulate the model by completing this code fragment::
+
+    time_max = 100
+    time_step = 0.1
+    time_hist = numpy.linspace(0., time_max, time_max / time_step + 1)
+    conc_hist = scipy.integrate.odeint(...)
+
+Finally, use ``matplotlib`` to plot the simulation results. You should see results similar to those below.
 
-* Write a Runge-kutta 4th order ODE integrator
+.. image:: ode-results.png
 
 
 Stochastic simulation
 ^^^^^^^^^^^^^^^^^^^^^
-In this exercise we will simulate the stochastic synthesis and degradation of a single RNA using the Gillespie algorith.
+In this exercise we will simulate the stochastic synthesis and degradation of a single RNA using the Gillespie algorithm.
 
 .. math::
 

intro_to_wc_modeling/cell_modeling/simulation/ode.py

@@ -1,11 +1,101 @@
 """ ODE simulation tutorial
 
 :Author: Jonathan Karr <jonrkarr@gmail.com>
-:Date: 2017-06-22
+:Date: 2017-06-23
 :Copyright: 2017, Karr Lab
 :License: MIT
 """
 
+import matplotlib
+matplotlib.use('Agg')
 import matplotlib.pyplot
 import numpy
+import scipy.integrate
 import os
+
+def d_conc_d_t(concs, time):
+    """ Calculate differentials for Goldbeter 1991 cell cycle model 
+    (`BIOMD0000000003 <http://www.ebi.ac.uk/biomodels-main/BIOMD0000000003>`_)
+
+    Args:
+        time (obj:`float`): time
+        concs (:obj:`numpy.ndarray`): array of current concentrations
+
+    Returns:
+        :obj:`numpy.ndarray`
+    """
+    cell = 1.0
+    vi = 0.025
+    kd = 0.01
+    vd = 0.25
+    Kd = 0.02
+    K1 = 0.005
+    V2 = 1.5
+    K2 = 0.005
+    K3 = 0.005
+    V4 = 0.5
+    K4 = 0.005
+    VM1 = 3.
+    VM3 = 1.
+    Kc = 0.5
+
+    C = concs[0] # cyclin
+    M = concs[1] # cdc2
+    X = concs[2] # cyclin protease
+
+    V1 = C * VM1 / (C + Kc)
+    V3 = M * VM3
+
+    r_cyclin_creation = cell * vi
+    r_default_cyclin_degradation = C * cell * kd
+    r_cdc2_triggered_cyclin_degradation = C * cell * vd * X / (C + Kd)
+    r_activation_of_cdc2 = cell * (1 - M) * V1 / (K1 - M + 1)
+    r_deactivation_of_cdc2 = cell * M * V2 / (K2 + M)
+    r_activation_of_cyclin_protease = cell * V3 * (1 - X) / (K3 - X + 1)
+    r_deactivation_of_cyclin_protease = cell * V4 * X / (K4 + X)
+
+    d_cyclin_dt = \
+        + r_cyclin_creation \
+        - r_default_cyclin_degradation \
+        - r_cdc2_triggered_cyclin_degradation
+    d_cdc2_dt = \
+        + r_activation_of_cdc2 \
+        - r_deactivation_of_cdc2
+    d_cyclin_protease_dt = \
+        + r_activation_of_cyclin_protease \
+        - r_deactivation_of_cyclin_protease
+
+    return numpy.array([
+        d_cyclin_dt, 
+        d_cdc2_dt,
+        d_cyclin_protease_dt
+        ])
+
+# initial conditions
+init_concs = numpy.array([0.01, 0.01, 0.01])
+
+# integrate model
+time_max = 100
+time_step = 0.1
+time_hist = numpy.linspace(0., time_max, time_max / time_step + 1)
+conc_hist = scipy.integrate.odeint(d_conc_d_t, init_concs, time_hist)
+
+# plot results
+line_cyclin, = matplotlib.pyplot.plot(time_hist, conc_hist[:, 0], 'b-', label='Cyclin')
+line_cdc2, = matplotlib.pyplot.plot(time_hist, conc_hist[:, 1], 'r-', label='Cdc2')
+line_cyclin_protease, = matplotlib.pyplot.plot(time_hist, conc_hist[:, 2], 'g-', label='Protease')
+matplotlib.pyplot.legend([line_cyclin, line_cdc2, line_cyclin_protease], ['Cyclin', 'Cdc2', 'Protease'])
+matplotlib.pyplot.xlim(0, time_max)
+matplotlib.pyplot.xlabel('Time (min)')
+matplotlib.pyplot.ylabel('Cyclin concentration and fraction of\nactive Cdc2 and cyclin protease')
+matplotlib.pyplot.gca().spines['top'].set_visible(False)
+matplotlib.pyplot.gca().spines['right'].set_visible(False)
+
+# display figure
+# matplotlib.pyplot.show()
+
+# save figure
+filename = os.path.join(os.path.dirname(__file__), '..', '..', '..', 'docs', 'tutorials',
+                        'cell_modeling', 'simulation', 'ode-results.png')
+matplotlib.pyplot.savefig(filename, transparent=True, bbox_inches='tight')
+matplotlib.pyplot.close()

requirements.txt

@@ -1,10 +1,11 @@
 cement                                                   # framework for building command line programs
 ipython                                                  # interactive interpret
 matplotlib                                               # plotting
-numpy                                                    # numerical computations
+numpy                                                    # numerical computing
 optlang                                                  # numerical optimization
 # pgmpy                                                  # graphical models package
 # scikit-learn                                           # data science package
+scipy                                                    # scientific computing
 sqlalchemy                                               # relational databases
 git+https://github.com/KarrLab/wc_lang.git#egg=wc_lang   # wc-modeling language
 git+https://github.com/KarrLab/wc_utils.git#egg=wc_utils # wc-modeling utilities