Running Time-Courses¶
Copasi enables users to simulate their model with a range of different solvers.
Create a model¶
Here we do our imports and create the model we use for the tutorial
[1]:
import os
import site
from pycotools3 import model, tasks, viz
working_directory = r'/home/ncw135/Documents/pycotools3/docs/source/Tutorials/timecourse_tutorial'
if not os.path.isdir(working_directory):
os.makedirs(working_directory)
copasi_file = os.path.join(working_directory, 'michaelis_menten.cps')
if os.path.isfile(copasi_file):
os.remove(copasi_file)
antimony_string = """
model michaelis_menten()
compartment cell = 1.0
var E in cell
var S in cell
var ES in cell
var P in cell
kf = 0.1
kb = 1
kcat = 0.3
E = 75
S = 1000
SBindE: S + E => ES; kf*S*E
ESUnbind: ES => S + E; kb*ES
ProdForm: ES => P + E; kcat*ES
end
"""
with model.BuildAntimony(copasi_file) as loader:
mm = loader.load(antimony_string)
mm
The BuildAntimony context manager is deprecated and will be removed in future versions. Please use model.loada instead.
[1]:
Model(name=michaelis_menten, time_unit=s, volume_unit=l, quantity_unit=mol)
Deterministic Time Course¶
[2]:
TC = tasks.TimeCourse(
mm, report_name='mm_simulation.txt',
end=1000, intervals=50, step_size=20
)
## check its worked
os.path.isfile(TC.report_name)
import pandas
df = pandas.read_csv(TC.report_name, sep='\t')
df.head()
[2]:
Time | [E] | [S] | [ES] | [P] | Values[kf] | Values[kb] | Values[kcat] | |
---|---|---|---|---|---|---|---|---|
0 | 0 | 75.00000 | 1000.000000 | 1.000000 | 1.000 | 0.1 | 1 | 0.3 |
1 | 20 | 2.00306 | 479.797000 | 73.996900 | 448.206 | 0.1 | 1 | 0.3 |
2 | 40 | 12.80010 | 62.104800 | 63.199900 | 876.695 | 0.1 | 1 | 0.3 |
3 | 60 | 75.13160 | 0.119830 | 0.868371 | 1001.010 | 0.1 | 1 | 0.3 |
4 | 80 | 75.99560 | 0.000604 | 0.004434 | 1001.990 | 0.1 | 1 | 0.3 |
When running a time course, you should ensure that the number of intervals times the step size equals the end time, i.e.:
- $$intervals \cdot step\_size = end$$
The default behaviour is to output all model variables as they can easily be filtered later in the Python environment. However, the metabolites
, global_quantities
and local_parameters
arguments exist to filter the variables that are simulated prior to running the time course.
[3]:
TC=tasks.TimeCourse(
mm,
report_name='mm_timecourse.txt',
end=100,
intervals=50,
step_size=2,
global_quantities = ['kf'], ##recall that antimony puts all parameters as global quantities
)
##check that we only kf as a global variables
pandas.read_csv(TC.report_name,sep='\t').head()
[3]:
Time | [E] | [S] | [ES] | [P] | Values[kf] | |
---|---|---|---|---|---|---|
0 | 0 | 75.00000 | 1000.000 | 1.0000 | 1.0000 | 0.1 |
1 | 2 | 1.10437 | 881.378 | 74.8956 | 45.7263 | 0.1 |
2 | 4 | 1.16265 | 836.516 | 74.8373 | 90.6465 | 0.1 |
3 | 6 | 1.22737 | 791.698 | 74.7726 | 135.5300 | 0.1 |
4 | 8 | 1.29962 | 746.928 | 74.7004 | 180.3720 | 0.1 |
An alternative and more convenient interface into the tasks.TimeCourse
class is using the model.Model.simulate
method. This is simply a wrapper and is used like so.
[4]:
data = mm.simulate(start=0, stop=100, by=0.1)
data.head()
[4]:
E | ES | P | S | |
---|---|---|---|---|
Time | ||||
0.0 | 75.00000 | 1.0000 | 1.00000 | 1000.000 |
0.1 | 1.05984 | 74.9402 | 3.02124 | 924.039 |
0.2 | 1.05665 | 74.9433 | 5.26956 | 921.787 |
0.3 | 1.05920 | 74.9408 | 7.51783 | 919.541 |
0.4 | 1.06175 | 74.9382 | 9.76601 | 917.296 |
This mechanism of running a time course has the advantage that 1) pycotools parses the data back into python in the form of a pandas.DataFrame
and 2) the column names are automatically pruned to remove the copasi reference information.
Visualization¶
[5]:
viz.PlotTimeCourse(TC)
[5]:
<pycotools3.viz.PlotTimeCourse at 0x16c80294358>
It is also possible to plot these on the same axis by specifying separate=False
[6]:
viz.PlotTimeCourse(TC, separate=False)
[6]:
<pycotools3.viz.PlotTimeCourse at 0x16ca06f91d0>
or to choose the y variables,
[7]:
viz.PlotTimeCourse(TC, y=['E', 'S'], separate=False)
viz.PlotTimeCourse(TC, y=['ES', 'P'], separate=False)
[7]:
<pycotools3.viz.PlotTimeCourse at 0x16ca0760748>
Plot in Phase Space¶
Choose the x variable to plot phase space. Same arguments apply as above.
[8]:
viz.PlotTimeCourse(TC, x='E', y='ES', separate=True)
[8]:
<pycotools3.viz.PlotTimeCourse at 0x16ca07c1b38>
Save to file¶
Use the savefig=True
option to save the figure to file and give an argument to the filename option to choose the filename.
[9]:
viz.PlotTimeCourse(TC, y=['S', 'P'], separate=False, savefig=True, filename='MyTimeCourse.png')
[9]:
<pycotools3.viz.PlotTimeCourse at 0x16ca0847cf8>
Alternative Solvers¶
Valid arguments for the method
argument of TimeCourse
are:
- deterministic
- direct
- gibson_bruck
- tau_leap
- adaptive_tau_leap
- hybrid_runge_kutta
- hybrid_lsoda
Copasi also includes a hybrid_rk45
solver but this is not yet supported by Pycotools. To use an alternative solver, pass the name of the solver to the method
argument.
Stochastic MM¶
For demonstrating simulation of stochastic time courses we build another michaelis-menten type reaction schema. We need to do this so we can set unit substance = item
, or in other words, change the model to particle numbers - otherwise there are too may molecules in the system to simulate a stochastic model
[10]:
copasi_file = os.path.join(working_directory, 'michaelis_menten_stochastic.cps')
antimony_string = """
model michaelis_menten()
compartment cell = 1.0;
var E in cell;
var S in cell;
var ES in cell;
var P in cell;
kf = 0.1;
kb = 1;
kcat = 0.3;
E = 75;
S = 1000;
SBindE: S + E => ES; kf*S*E;
ESUnbind: ES => S + E; kb*ES;
ProdForm: ES => P + E; kcat*ES;
unit substance = item;
end
"""
with model.BuildAntimony(copasi_file) as loader:
mm = loader.load(antimony_string)
The BuildAntimony context manager is deprecated and will be removed in future versions. Please use model.loada instead.
Run a Time Course Using Direct Method¶
[11]:
data = mm.simulate(0, 100, 1, method='direct')
data.head(n=10)
[11]:
E | ES | P | S | |
---|---|---|---|---|
Time | ||||
0 | 75 | 1 | 1 | 1000 |
1 | 0 | 76 | 18 | 908 |
2 | 2 | 74 | 36 | 892 |
3 | 0 | 76 | 57 | 869 |
4 | 1 | 75 | 81 | 846 |
5 | 0 | 76 | 100 | 826 |
6 | 0 | 76 | 122 | 804 |
7 | 1 | 75 | 143 | 784 |
8 | 1 | 75 | 167 | 760 |
9 | 1 | 75 | 200 | 727 |
Plot stochastic time course¶
Note that we can also use the pandas
, matplotlib
and seaborn
libraries for plotting
[12]:
import matplotlib
import seaborn
seaborn.set_context('notebook')
data.plot()
[12]:
<matplotlib.axes._subplots.AxesSubplot at 0x16ca0810c88>
Notice how similar the stochastic simulation is to the deterministic. As the number of molecules being simulated increases, the stochastic simulation converges to the deterministic solution. Another way to put it, stochastic effects are greatest when simulating small molecule numbers