Integrating ordinary differential equations


Contents

Calling odeint
Specifying parameters
An example of a genetic switch

Calling odeint

Within SciPy, odeint is used to integrate systems of (first-order) ordinary differential equations. In IPython, we must first import the command:

from scipy.integrate import odeint

We can then use, for example,

y= odeint(model, initial_values, t)

to simulate a model of k variables or chemical species. The model is specifed by the function model, which returns a list where each item in the list is the time-derivative of one of the variables. This function must be a function of y and of t in that order:

def model(y, t):

For example, to simulate the inactivation of a protein from its activate state Aa to its inactive state Ai, we would use:

import numpy as np

def model(y,t):
  dydt= np.empty(len(y))
  dydt[0]= -k*y[0]  # rate equation for Aa
  dydt[1]= k*y[0]   # rate equation for Ai
  return dydt

The initial values of the variables to be simulated are specified in the list initial_values (and would be the initial numbers of Aa and Ai in our example), and the time points for which values of the variables should be returned are given in the array t. For example,

t= np.arange(0, 1000, 10)

The SciPy function odeint returns a matrix which has k columns, one for each of the k variables, and a row for each time point specified in t.

Specifying parameters

If you wish to define the model function to depend on variables other than y and t, then these variables should be passed from odeint in a tuple using the args= option. Typically, these variables would be parameters for the model, such as k in the inactivation example above. If the model function has, for example, two additional parameters, such as

def model(y, t, a, b):

where a and b are the additional parameters, then we call odeint with

parameters= (0.1, 0.001)
y= odeint(model, init, t, args= parameters)

or

y= odeint(model, init, t, args= (0.1, 0.001))

more succinctly. Specifying just one parameter requires

y= odeint(model, init, t, args= (0.1,))

with a trailing comma to denote a tuple. Note that these parameters must be specified after t in the definition of the model function.

An example of a genetic switch

In this example, we will draw the nullclines and direction arrows, find the steady-states, and give a sample trajectory for a two-variable model of a genetic switch. We will consider positive autoregulation and the variables are protein, $P$, and mRNA, $M$.

The differential equations are

$${dP}/{dt} = M - d_P P$$

and

$${dM}/{dt} = v + u P^4/{K^4 + P^4} - d_M M$$

with basal transcription, $v$, and positive feedback of strength determined by $u$.

The script to generate the figure is:


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import fsolve

# define differential equations
(the function is defined in this way to enable the easy calculation of the direction arrows)
def genswitch(y, t):
     p, m = y
     dpdt= m - dp*p
     dmdt= v + u*p**n/(K**n + p**n) - dm*m
     return np.array([dpdt, dmdt])

# define steady-state equation for P
def steadystateeqn(p):
     return v + u*p**n/(K**n + p**n) - dm*dp*p


# specify parameters
n= 4
v= 0.01
dm= 0.008
dp= 0.002
K= 2.0e3
u= 0.06

# specify range of variables
npts= 25
pmin, pmax= 0, 5.0e3
mmin, mmax= 0, 10.0
p= np.linspace(pmin, pmax, npts)
m= np.linspace(mmin, mmax, npts)

plt.figure()
# plot nullclines
plt.plot(p, v/dm + u*p**n/(K**n + p**n)/dm, 'b')
plt.plot(p, dp*p, 'b')

# find steady-state solutions using an initial guess
sol1= fsolve(steadystateeqn, K/100)[0]
sol2= fsolve(steadystateeqn, K)[0]
sol3= fsolve(steadystateeqn, 5*K)[0]
plt.plot(sol1, dp*sol1, 'go')
plt.plot(sol2, dp*sol2, 'ro')
plt.plot(sol3, dp*sol3, 'go')

# find and plot direction arrows (dP, dM)
P, M= np.meshgrid(p, m)
dP, dM= genswitch([P, M], 0.0)
# calculate magnitude of total rate of change
norm= np.sqrt(dP**2 + dM**2)
# colour of direction arrow is magnitude of total rate of change
plt.streamplot(P, M, dP, dM, color= norm, cmap= plt.cm.cool, density= 2)

# label and show figure
plt.xlim([pmin, pmax])
plt.ylim([mmin, mmax])
plt.xlabel('P')
plt.ylabel('M')
plt.show()