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>
../_images/Tutorials_Timecourse_10_1.png
../_images/Tutorials_Timecourse_10_2.png
../_images/Tutorials_Timecourse_10_3.png
../_images/Tutorials_Timecourse_10_4.png

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>
../_images/Tutorials_Timecourse_12_1.png

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>
../_images/Tutorials_Timecourse_14_1.png
../_images/Tutorials_Timecourse_14_2.png

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>
../_images/Tutorials_Timecourse_16_1.png

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>
../_images/Tutorials_Timecourse_19_1.png

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>
../_images/Tutorials_Timecourse_25_1.png

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