# 2
#
#
# See the exercise "Repressilator.pdf" from Repressilator.html
# in http://www.physics.cornell.edu/~myers/teaching/ComputationalMethods/ComputerExercises/
#
#
import RepressilatorAnimation
RepressilatorAnimation.ClearScene()
animator = RepressilatorAnimation.RepressilatorAnimator()
#
import scipy, scipy.integrate
import random
#
class Chemical:
"""Chemical is a class describing chemicals and/or chemical complexes.
Chemical.amount: the amount of the Chemical in the cell.
Chemical.reactions: a list of reactions whose rates are changed if
Chemical.amount is altered."""
def __init__(self, amount):
# Chemical.amount represents the amount of this Chemical in the cell
self.amount = amount
# Chemical.reactions is a list of reactions whose rates are changed if
# Chemical.amount is altered.
self.reactions = []
#
class DegradationReaction:
"""DegradationReaction describes the removal of one molecule of
the specified substrate, with specified rate_constant.
Overall rate for this reaction to occur is
substrate.amount * rate_constant."""
def __init__(self, substrate, rate_constant):
self.stoichiometry = {substrate: -1}
self.substrate = substrate
self.rate_constant = rate_constant
substrate.reactions.append(self)
def GetRate(self):
return self.substrate.amount * self.rate_constant
#
class CatalyzedSynthesisReaction:
"""CatalyzedSynthesisReaction describes the synthesis of product in
the presence of a catalyst: catalyst -> catalyst + product,
with specified rate_constant.
Overall rate for this reaction to occur is
catalyst.amount * rate_constant."""
def __init__(self, catalyst, product, rate_constant):
self.stoichiometry = {product: 1}
self.catalyst = catalyst
self.rate_constant = rate_constant
product.reactions.append(self)
def GetRate(self):
return self.catalyst.amount * self.rate_constant
#
class HeterodimerBindingReaction:
"""HeterodimerBindingReaction describes the binding of two distinct
types of chemicals, A and B, to form a product dimer: A + B -> dimer,
with specified rate_constant.
Overall rate for this reaction to occur is
A.amount * B.amount * rate_constant."""
def __init__(self, A, B, dimer, rate_constant):
self.stoichiometry = {A: -1, B: -1, dimer: 1}
self.A = A
self.B = B
self.rate_constant = rate_constant
A.reactions.append(self)
B.reactions.append(self)
dimer.reactions.append(self)
def GetRate(self):
return self.A.amount * self.B.amount * self.rate_constant
#
class HeterodimerUnbindingReaction:
"""HeterodimerBindingReaction describes the unbinding of a
heterodimer into two distinct types of chemicals, A and B:
dimer -> A + B, with specified rate_constant.
Overall rate for this reaction to occur is
dimer.amount * rate_constant."""
def __init__(self, dimer, A, B, rate_constant):
self.stoichiometry = {A: 1, B: 1, dimer: -1}
self.dimer = dimer
self.rate_constant = rate_constant
A.reactions.append(self)
B.reactions.append(self)
dimer.reactions.append(self)
def GetRate(self):
return self.dimer.amount * self.rate_constant
#
# ---------------------------------------------------------------------
#
class Repressilator:
"""
Repressilator is a base class for the three-gene repressilator
system introduced by Elowitz and Leibler, ref. All reactions
in the stochastic version of the Repressilator are implemented, with
rates as given in the paper.
The set of Chemical states (all 3 mRNAs, all 3 proteins, and all 9
promoter states) have a defined ordering, as stored in the
dictionary chemIndex, which maps Chemicals to integer indices in
the ordering. This ordering is necessary to communicate with utilities
(such as scipy.integrate.odeint) that require chemical amounts to
be stored as a 1D array. The GetStateVector method returns the
Chemical amounts as a scipy array in the appropriate order, and
the SetStateVector method sets the Chemical amounts based on a
supplied array.
"""
def __init__(self,
mRNA_degradation_rate=scipy.log(2.)/120.,
protein_degradation_rate=scipy.log(2.)/600.,
translation_rate=0.167,
unocc_transcription_rate=0.5,
occ_transcription_rate=5.0e-04,
P_binding=1.0,
P1_unbinding=224.0,
P2_unbinding=9.0):
self.t = 0.
self.chemIndex = {}
self.reactions = []
self.mRNA_degradation_rate = mRNA_degradation_rate
self.protein_degradation_rate = protein_degradation_rate
self.translation_rate = translation_rate
self.unocc_transcription_rate = unocc_transcription_rate
self.occ_transcription_rate = occ_transcription_rate
self.P_binding = P_binding
self.P1_unbinding = P1_unbinding
self.P2_unbinding = P2_unbinding
self.lacI = Chemical(30.0)
self.tetR = Chemical(20.0)
self.cI = Chemical(10.0)
self.LacI = Chemical(300.0)
self.TetR = Chemical(200.0)
self.CI = Chemical(100.0)
self.PlacI = Chemical(1.0)
self.PtetR = Chemical(1.0)
self.PcI = Chemical(1.0)
self.PlacI_CI = Chemical(0.0)
self.PtetR_LacI = Chemical(0.0)
self.PcI_TetR = Chemical(0.0)
self.PlacI_CI_CI = Chemical(0.0)
self.PtetR_LacI_LacI = Chemical(0.0)
self.PcI_TetR_TetR = Chemical(0.0)
self.mRNAs = [self.lacI, self.tetR, self.cI]
self.proteins = [self.LacI, self.TetR, self.CI]
self.P0 = [self.PlacI, self.PtetR, self.PcI]
self.P1 = [self.PlacI_CI, self.PtetR_LacI, self.PcI_TetR]
self.P2 = [self.PlacI_CI_CI, self.PtetR_LacI_LacI, self.PcI_TetR_TetR]
for chem in self.mRNAs:
self.AddChemical(chem)
for chem in self.proteins:
self.AddChemical(chem)
for chem in self.P0:
self.AddChemical(chem)
for chem in self.P1:
self.AddChemical(chem)
for chem in self.P2:
self.AddChemical(chem)
for i in range(3):
self.AddReaction(\
DegradationReaction(self.mRNAs[i], mRNA_degradation_rate))
self.AddReaction(\
DegradationReaction(self.proteins[i],protein_degradation_rate))
self.AddReaction(\
CatalyzedSynthesisReaction(self.mRNAs[i], self.proteins[i], \
translation_rate))
self.AddReaction(\
CatalyzedSynthesisReaction(self.P0[i],
self.mRNAs[i],
unocc_transcription_rate))
self.AddReaction(\
CatalyzedSynthesisReaction(self.P1[i],
self.mRNAs[i],
occ_transcription_rate))
self.AddReaction(\
CatalyzedSynthesisReaction(self.P2[i],
self.mRNAs[i],
occ_transcription_rate))
j = (i+1)%3
self.AddReaction(\
HeterodimerBindingReaction(self.P0[i],
self.proteins[j],
self.P1[i],
P_binding))
self.AddReaction(\
HeterodimerBindingReaction(self.P1[i],
self.proteins[j],
self.P2[i],
P_binding))
self.AddReaction(\
HeterodimerUnbindingReaction(self.P1[i],
self.P0[i],
self.proteins[j],
P1_unbinding))
self.AddReaction(\
HeterodimerUnbindingReaction(self.P2[i],
self.P1[i],
self.proteins[j],
P2_unbinding))
self.rates = scipy.zeros(len(self.reactions), float)
for rIndex, r in enumerate(self.reactions):
self.rates[rIndex] = r.GetRate()
def AddChemical(self,chemical):
self.chemIndex[chemical] = len(self.chemIndex)
def GetChemicalIndex(self, chemical):
return self.chemIndex[chemical]
def AddReaction(self,reaction):
self.reactions.append(reaction)
def GetStateVector(self):
c = scipy.zeros(len(self.chemIndex), float)
for chem, index in list(self.chemIndex.items()):
c[index] = chem.amount
return c
def SetFromStateVector(self, c):
for chem, index in list(self.chemIndex.items()):
chem.amount = c[index]
#
class StochasticRepressilator (Repressilator):
"""
StochasticRepressilator is a stochastic implementation of
the three-gene repressilator system, with time evolution implemented
using Gillespie's Direct Method, as described in detail in Step.
"""
def ComputeReactionRates(self):
"""ComputeReactionRates computes the current rate for every
reaction defined in the network, and stores the rates in self.rates."""
for index, r in enumerate(self.reactions):
self.rates[index] = r.GetRate()
def Step(self, dtmax):
"""Step(self, dtmax) implements Gillespie's Direct Simulation Method,
executing at most one reaction and returning the time increment
required for that reaction to take place. If no reaction is executed,
the specified maximal time increment dtmax is returned.
(1) all reaction rates are computed
(2) a total rate for all reactions is found
(3) a random time is selected, to be drawn from an exponential
distribution with mean value given by the inverse of the total
rate, e.g., ran_dtime = -scipy.log(1.-random.random())/total_rate
(4) if the random time is greater than the time interval under
consideration (dtmax), then no reaction is executed and dtmax
is returned
(5) otherwise, a reaction is chosen at random with relative
probabilities given by the relative reaction rates;
this is done by
(5a) uniformly drawing a random rate from the interval from
[0., total rate)
(5b) identifying which reaction rate interval corresponds to
the randomly drawn rate, e.g.,
|<-------------------total rate---------------------->|
|<----r0----->|<-r1->|<--r2-->|<-----r3----->|<--r4-->|
| X |
Randomly drawn rate X lands in interval r3
(6) the chosen reaction is executed
(7) the time at which the reaction is executed is returned
"""
self.ComputeReactionRates()
total_rate = sum(self.rates)
# get exponentially distributed time
ran_time = -scipy.log(1.-random.random())/total_rate
if ran_time > dtmax:
return dtmax
# get uniformly drawn rate in interval defined by total_rate
ran_rate = total_rate*random.random()
# find interval corresponding to random rate
reac_index = len(self.rates) - sum(scipy.cumsum(self.rates) > ran_rate)
reaction = self.reactions[reac_index]
# execute specified reaction
for chem, dchem in list(reaction.stoichiometry.items()):
chem.amount += dchem
# return time at which reaction takes place
return ran_time
def Run(self, T, delta_t=0.0):
"""Run(self, T, delta_t) runs the StochasticRepressilator for
a specified time interval T, returning the trajectory at
specified time intervals delta_t (or, if delta_t == 0.0, after
every reaction)
"""
tfinal = self.t + T
if delta_t == 0.:
ts = [self.t]
trajectory = [self.GetStateVector()]
while self.t < tfinal:
dt = self.Step(tfinal-self.t)
self.t += dt
ts.append(self.t)
trajectory.append(self.GetStateVector())
return scipy.array(ts), scipy.array(trajectory)
else:
eps=1.0e-06
ts = scipy.arange(0., T+eps, delta_t)
trajectory = scipy.zeros((len(ts), len(self.chemIndex)),
float)
trajectory[0] = self.GetStateVector()
tindex = 0
while self.t < tfinal:
dt = self.Step(ts[tindex+1]-self.t)
self.t += dt
if self.t >= ts[tindex+1]:
tindex += 1
for chem, cindex in list(self.chemIndex.items()):
trajectory[tindex][cindex] = chem.amount
return ts, trajectory
#
class DeterministicRepressilator (Repressilator):
"""
DeterministicRepressilator is a deterministic implementation of
the three-gene repressilator system, with time evolution implemented
by summing up all reactions as appropriate to form a differential
equation describing the time rate of change of all chemical
constituents in the model.
"""
def dcdt(self, c, t):
"""dcdt(self, c, t) returns the instantaneous time rate of
change of the DeterministicRepressilator system, given chemical
concentration state vector c and current time t, for use in
integration by scipy.integrate.odeint.
dcdt loops through all reactions defined in the Repressilator
system, computes the rates of those reactions, and increments
those elements in a dc_dt array that are affected by the
reaction under consideration.
the fully assembled dc_dt array is returned by this method.
"""
pass
def Run(self, tmax, dt):
"""Run(self, tmax, dt) integrates the DeterministicRepressilator
for a time tmax, returning the trajectory at time steps as
specified by dt, by calling scipy.integrate.odeint with the
self.dcdt method describing the time derivative of the system
Run should return the time array on which the trajectory is computed,
along with the trajectory corresponding to those time points.
"""
pass
#
def RunStochasticRepressilator(T=100., dt=1., plots=False, plotPromoter=False):
"""RunStochasticRepressilator(tmax, dt, plots=False, plotPromoter=False)
creates and runs a StochasticRepressilator for the specified time
interval T, returning the trajectory in time increments dt,
optionally using pylab to make plots of mRNA,
protein, and promoter amounts along the trajectory.
"""
sr = StochasticRepressilator()
sts, straj = sr.Run(T, dt)
curvetypes = ['r-', 'g-', 'b-']
if plots:
import pylab
pylab.figure(1)
for i in range(3):
pylab.plot(sts, straj[:,sr.chemIndex[sr.mRNAs[i]]], curvetypes[i])
pylab.figure(2)
for i in range(3):
pylab.plot(sts, straj[:,sr.chemIndex[sr.proteins[i]]], curvetypes[i])
if plotPromoter:
pylab.figure(3)
for i in range(3):
promoter_state = (0.99+0.01*i)*\
(straj[:,sr.chemIndex[sr.P1[i]]]\
+2.*straj[:,sr.chemIndex[sr.P2[i]]])
pylab.plot(sts, promoter_state, curvetypes[i])
pylab.show()
return sr, sts, straj
#
def RunDeterministicRepressilator(T=100., dt=1., plots=False,
plotPromoter=False):
"""RunDeterministicRepressilator(tmax, dt, plots=False, plotPromoter=False)
creates and runs a DeterministicRepressilator for the specified time
interval T, returning the trajectory in time increments dt,
optionally using pylab to make plots of mRNA, protein,
and promoter amounts along the trajectory.
"""
pass
#
def CompareDetSto(tmax=100., dt=1., RNAFactor=1., TelegraphFactor=1.,
plots=True, plotPromoter=False):
"""
CompareDetSto compares the Deterministic and Stochastic variants of
the Repressilator by running each, and plotting the trajectories
alongside each other. It is useful to return both of the Repressilator
instances as well as both of the trajectories created.
CompareDetSto(tmax, dt) runs each Repressilator for tmax seconds,
returning the trajectory for each in increments of dt.
The optional arguments RNAFactor and TelegraphFactor are used to explore
the effects of shot noise and telegraph noise on the dynamics of
the StochasticRepressilator; in particular,
- initial RNA concentrations and the transcription rate are multiplied by
RNAFactor, and the translation rate is divided by RNAFactor, thereby
increasing RNA concentrations without otherwise affecting the continuum
equations
- promoter binding and unbinding rates are multiplied by TelegraphFactor
The argument plots is a boolean variable affecting whether or not
the RNA and protein concentrations for each Repressilator are plotted.
The argument plotPromoter is a boolean variable affecting whether or not
the various promoter states are plotted. Because the promoter states
switch discretely (in the StochasticRepressilator) between 0, 1, and 2,
it is useful to separate them from each other so what the state of
each can be ascertained. For a trajectory straj, such a separation
might look like:
if plotPromoter:
pylab.figure(3)
pylab.plot(sts, 0.99*(straj[:,9]+2.*straj[:,12]), 'r-')
pylab.plot(sts, 1.0*(straj[:,10]+2.*straj[:,13]), 'b-')
pylab.plot(sts, 1.01*(straj[:,11]+2.*straj[:,14]), 'g-')
"""
pass
#
# ---------------------------------------------------------------------
#
def Animate(trajectory):
"""Animate(trajectory) will animate a Repressilator trajectory using
objects from VPython/visual to represent protein, mRNA and promoter
states. A trajectory should be a scipy array consisting of a number
of timesteps, where each time step is a 15-element state vector for
the Repressilator, in the defined ordering.
NOTE: VPython/visual uses threads to separate control of the
graphics from the Python interpreter. The VPython thread seems
to conflict with threads initiated by pylab, and can lead to
crashing of the Python interpreter. Therefore it is generally not a
good idea to do both VPython animations and pylab plots within
the same session.
Since trajectories can take a significant amount of time to generate,
it is sometimes useful to save a generated trajectory to a file
(especially given the propensity of thread-induced crashes as described
above). Trajectories are simply arrays, so they can be written to
and read from disk using the ArrayIO module in the Scientific.IO package:
from Scientific.IO import ArrayIO
# write trajectory traj to file traj.dat
ArrayIO.writeArray(traj, 'traj.dat')
# read file traj.dat and store trajectory as newtraj
newtraj = ArrayIO.readArray('traj.dat')
"""
#import RepressilatorAnimation
#RepressilatorAnimation.ClearScene()
#animator = RepressilatorAnimation.RepressilatorAnimator()
animator.AnimateTrajectory(trajectory)
#
# Copyright (C) Cornell University
# All rights reserved.
# Apache License, Version 2.0