PyCoTools

PyCoTools is a python package that was developed as an alternative interface into COPASI, simulation software for modelling biochemical systems. The PyCoTools paper can be found here and describes in detail the intentions and functionality of PyCoTools. There are some important differences between the PyCoTools version that is described in the publication and the current version. The first and most important is that PyCoTools is now a python 3 only package. If using Python 2.7 you should create a virtual Python 3.6 environment using conda or virtualenv. My preference is conda. The other major difference is the interface into COPASI’s parameter estimation task which has a new interface and has been enhanced to support all features of COPASI’s parameter estimation task.

Installation

Use:

$ pip install pycotools3

Remember to source activate your python 3.6 environment if you need to.

To install from source:

$ git clone https://github.com/CiaranWelsh/pycotools3.git
$ cd pycotools3
$ python setup.py install

The procedure is the same in linux, mac and windows.

Documentation

This is a guide to PyCoTools version >2.0.1.

Getting started

As PyCoTools only provides an alternative interface into some of COPASI’s tasks, if you are unfamiliar with COPASI <http://copasi.org/>_ then it is a good idea to become acquainted, prior to proceeding. As much as possible, arguments to PyCoTools functions follow the corresponding option in the COPASI user interface.

In addition to COPASI, PyCoTools depends on tellurium which is a Python package for modelling biological systems. While tellurium and COPASI have some of the same features, generally they are complementary and productivity is enhanced by using both together.

More specifically, tellurium uses antimony strings to define a model which is then converted into SBML. PyCoTools provides the model.BuildAntimony class which is a wrapper around this tellurium feature, which creates a Copasi model and parses it into a PyCoTools model.Model.

Since antimony is described elsewhere we will focus here on using antimony to build a copasi model.

Build a model with antimony

import site, os
site.addsitedir('D:\pycotools3')
from pycotools3 import model

working_directory = os.path.abspath('')
copasi_filename = os.path.join(working_directory, 'NegativeFeedbackModel.cps')
antimony_string =  '''
    model negative_feedback()
        // define compartments
        compartment cell = 1.0
        //define species
        var A in cell
        var B in cell
        //define some global parameter for use in reactions
        vAProd = 0.1
        kADeg  = 0.2
        kBProd = 0.3
        kBDeg  = 0.4
        //define initial conditions
        A      = 0
        B      = 0
        //define reactions
        AProd:  => A; cell*vAProd
        ADeg: A =>  ; cell*kADeg*A*B
        BProd:  => B; cell*kBProd*A
        BDeg: B =>  ; cell*kBDeg*B
    end
    '''
with model.BuildAntimony(copasi_filename) as loader:
    negative_feedback = loader.load(antimony_string)
print(negative_feedback)
assert os.path.isfile(copasi_filename)

Model(name=negative_feedback, time_unit=s, volume_unit=l, quantity_unit=mol)

Create an antmiony string from an existing model

The Copasi user interface is an excellant way of constructing a model and it is easy to convert this model into an antimony string that can be pasted into a document.

print(negative_feedback.to_antimony())

// Created by libAntimony v2.9.4
function Constant_flux__irreversible(v)
  v;
end

function Function_for_ADeg(A, B, kADeg)
  kADeg*A*B;
end

function Function_for_BProd(A, kBProd)
  kBProd*A;
end


model *negative_feedback()

  // Compartments and Species:
  compartment cell;
  species A in cell, B in cell;

  // Reactions:
  AProd:  => A; cell*Constant_flux__irreversible(vAProd);
  ADeg: A => ; cell*Function_for_ADeg(A, B, kADeg);
  BProd:  => B; cell*Function_for_BProd(A, kBProd);
  BDeg: B => ; cell*kBDeg*B;

  // Species initializations:
  A = 0;
  B = 0;

  // Compartment initializations:
  cell = 1;

  // Variable initializations:
  vAProd = 0.1;
  kADeg = 0.2;
  kBProd = 0.3;
  kBDeg = 0.4;

  // Other declarations:
  const cell, vAProd, kADeg, kBProd, kBDeg;
end

One paradigm of model development is to use antimony to ‘hard code’ perminent changes to the model and the Copasi user interface for experimental changes. The Model.open() method is useful for this paradigm as it opens the model with whatever configurations have been defined.

## negative_feedback.open()

cmd D:\pycotools3\pycotools3\COPASI\windows\CopasiUI.exe "D:\pycotools3\docs\source\NegativeFeedbackModel.cps"

Simulate a time course

Since we have used an antimony string, we can simulate this model with either tellurium or Copasi. Simulating with tellurium uses a library called roadrunner which is described in detail elsewhere. To run a simulation with Copasi we need to configure the time course task, make the task executable (i.e. tick the check box in the top right of the time course task) and run the simulation with CopasiSE. This is all taken care of by the tasks.TimeCourse class.

from pycotools3 import tasks
time_course = tasks.TimeCourse(negative_feedback, end=100, step_size=1, intervals=100)
time_course

<pycotools3.tasks.TimeCourse at 0x1c9ac327a20>

However a more convenient interface is provided by the model.simulate method, which is a wrapper around tasks.TimeCourse which additionally parses the resulting data from file and returns a pandas.DataFrame

from pycotools3 import tasks
sim_data = negative_feedback.simulate(0, 100, 1) ##start, end, by
sim_data.head()

	A	B
Time
0	0.000000	0.000000
1	0.099932	0.013181
2	0.199023	0.046643
3	0.295526	0.093275
4	0.387233	0.147810

The results are saved in a file defined by the report_name option, which defaults to timecourse.txt in the same directory as the copasi model.

Visualise a time course

PyCoTools also provides facilities for visualising simulation output. To plot a timecourse, pass the task.TimeCourse object to the viz.PlotTimeCourse object.

from pycotools3 import viz
viz.PlotTimeCourse(time_course, savefig=True)

<pycotools3.viz.PlotTimeCourse at 0x1c9ac3279e8>

More information about running time courses with PyCoTools and Copasi can be found in the time course tutorial

Tutorials

In this section of the documentation I provide detailed explainations of how PyCoTools works, with examples. The tutorials are split into sections which are linked to below.

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

import os
import site
site.addsitedir('D:\pycotools3')
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

Model(name=michaelis_menten, time_unit=s, volume_unit=l, quantity_unit=mol)

Deterministic Time Course

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()

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.

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()

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.

data = mm.simulate(start=0, stop=100, by=0.1)
data.head()

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

viz.PlotTimeCourse(TC)

It is also possible to plot these on the same axis by specifying separate=False

viz.PlotTimeCourse(TC, separate=False)

or to choose the y variables,

viz.PlotTimeCourse(TC, y=['E', 'S'], separate=False)
viz.PlotTimeCourse(TC, y=['ES', 'P'], separate=False)

Plot in Phase Space

Choose the x variable to plot phase space. Same arguments apply as above.

viz.PlotTimeCourse(TC, x='E', y='ES', separate=True)

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.

viz.PlotTimeCourse(TC, y=['S', 'P'], separate=False, savefig=True, filename='MyTimeCourse.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

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)

Run a Time Course Using Direct Method

data = mm.simulate(0, 100, 1, method='direct')
data.head(n=10)

	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

import matplotlib
import seaborn
seaborn.set_context('notebook')
data.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1fa3c225d68>

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

Insert Parameters

Parameters can be inserted automatically into a Copasi model from python code using PyCoTools

Build a demonistration model

While antimony or the COPASI user interface are the preferred ways to build a model, PyCoTools does have a mechanism for constructing COPASI models. For variation and demonstration, this method is used here.

import os
import site
site.addsitedir('D:\pycotools3')
from pycotools3 import model, tasks, viz
## Choose a working directory for model
working_directory = os.path.abspath('')
copasi_file = os.path.join(working_directory, 'MichaelisMenten.cps')

if os.path.isfile(copasi_file):
    os.remove(copasi_file)


kf = 0.01
kb = 0.1
kcat = 0.05
with model.Build(copasi_file) as m:
    m.name = 'Michaelis-Menten'
    m.add('compartment', name='Cell')

    m.add('metabolite', name='P', concentration=0)
    m.add('metabolite', name='S', concentration=30)
    m.add('metabolite', name='E', concentration=10)
    m.add('metabolite', name='ES', concentration=0)

    m.add('reaction', name='S bind E', expression='S + E -> ES', rate_law='kf*S*E',
          parameter_values={'kf': kf})

    m.add('reaction', name='S unbind E', expression='ES -> S + E', rate_law='kb*ES',
         parameter_values={'kb': kb})

    m.add('reaction', name='ES produce P', expression='ES -> P + E', rate_law='kcat*ES',
          parameter_values={'kcat': kcat})

mm = model.Model(copasi_file)
mm

Model(name=Michaelis-Menten, time_unit=s, volume_unit=ml, quantity_unit=mmol)

Insert Parameters from Python Dictionary

params = {'E': 100,
          'P': 150}

## Insert into model
I = model.InsertParameters(mm, parameter_dict=params)
##format the parameters for displaying nicely
I.parameters.index = ['Parameter Value']
I.parameters.transpose()

	Parameter Value
E	100
P	150

Alternatively use inplace=True argument (analogous to the pandas library) to modify the object inplace, rather than needing to assign

model.InsertParameters(mm, parameter_dict=params, inplace=True)

<pycotools3.model.InsertParameters at 0x1b41f67a9e8>

Insert Parameters from Pandas DataFrame

import pandas
params = {'(S bind E).kf': 50,
          '(S unbind E).kb': 96}
df = pandas.DataFrame(params, index=[0])
df

	(S bind E).kf	(S unbind E).kb
0	50	96

model.InsertParameters(mm, df=df, inplace=True)

<pycotools3.model.InsertParameters at 0x1b41f67a048>

Insert Parameters from Parameter Estimation Output

First we’ll get some parameter estimation data by fitting a model to simulated data.

fname = os.path.join(os.path.abspath(''), 'timecourse.txt')
data = mm.simulate(0, 50, 1, report_name=fname)
assert os.path.isfile(fname)

with  tasks.ParameterEstimation.Context(copasi_file, fname, context='s', parameters='a') as context:
    context.randomize_start_values = True
    context.lower_bound = 0.01
    context.upper_bound = 100
    context.run_mode = True
    config = context.get_config()
PE = tasks.ParameterEstimation(config)

mm

Model(name=Michaelis-Menten, time_unit=s, volume_unit=ml, quantity_unit=mmol)

Now we can insert the estimated parameters using:

##index=0 for best parameter set (i.e. lowest RSS)
model.InsertParameters(mm, parameter_path=PE.results_directory, index=0, inplace=True)

---------------------------------------------------------------------------
InputError                                Traceback (most recent call last)
<ipython-input-24-d908d0566dc2> in <module>()
      1 ##index=0 for best parameter set (i.e. lowest RSS)
----> 2 model.InsertParameters(mm, parameter_path=PE.results_directory, index=0, inplace=True)

D:\pycotools3\pycotools3\model.py in __init__(self, model, parameter_dict, df, parameter_path, index, quantity_type, inplace)
   4674         self._do_checks()
   4675
-> 4676         self.model = self.insert()
   4677         if self.inplace:
   4678             self.model.save()

D:\pycotools3\pycotools3\model.py in insert(self)
   4867
   4868         """
-> 4869         self.model = self.insert_locals()
   4870         self.model = self.insert_compartments()
   4871         self.model = self.insert_global_quantities()

D:\pycotools3\pycotools3\model.py in insert_locals(self)
   4774         # print self.parameters
   4775
-> 4776         locals = [j for i in self.model.reactions for j in i.parameters if
   4777                   j.global_name in list(self.parameters.keys())]
   4778         if locals == []:

D:\pycotools3\pycotools3\model.py in <listcomp>(.0)
   4775
   4776         locals = [j for i in self.model.reactions for j in i.parameters if
-> 4777                   j.global_name in list(self.parameters.keys())]
   4778         if locals == []:
   4779             return self.model

D:\pycotools3\pycotools3\cached_property.py in __get__(self, obj, cls)
     38         if obj is None:
     39             return self
---> 40         value = obj.__dict__[self.func.__name__] = self.func(obj)
     41         return value
     42

D:\pycotools3\pycotools3\model.py in parameters(self)
   4760
   4761         if self.parameter_path != None:
-> 4762             P = viz.Parse(self.parameter_path, copasi_file=self.model.copasi_file)
   4763             if isinstance(self.index, int):
   4764                 return pandas.DataFrame(P.data.iloc[self.index]).transpose()

D:\pycotools3\pycotools3\viz.py in __init__(self, cls_instance, log10, copasi_file, alpha, rss_value, num_data_points)
    519             raise errors.InputError('{} not in {}'.format(
    520                 self.cls_instance,
--> 521                 accepted_types)
    522             )
    523

InputError: {'MichaelisMenten': 'D:\\pycotools3\\docs\\source\\Tutorials\\Problem1\\Fit1\\MichaelisMenten\\ParameterEstimationData'} not in [<class 'pycotools3.tasks.TimeCourse'>, <class 'pycotools3.tasks.Scan'>, <class 'pycotools3.tasks.ParameterEstimation'>, <class 'str'>, <class 'pycotools3.viz.Parse'>, <class 'pycotools3.tasks.ProfileLikelihood'>, <class 'pandas.core.frame.DataFrame'>, <class 'pycotools3.tasks.ChaserParameterEstimations'>]

Insert Parameters using the model.Model().insert_parameters method

The same means of inserting parameters can be used from the model object itself

mm.insert_parameters(parameter_dict=params, inplace=True)

Change parameters using model.Model().set

Individual parameters can also be changed using the set method. For example, we could set the metabolite with name S concentration or particle numbers to 55

mm.set('metabolite', 'S', 55, 'name', 'concentration')

## or

mm.set('metabolite', 'S', 55, 'name', 'particle_numbers')

Parameter Scan

Copasi supports three types of scan, a regular parameter scan, a repeat scan and sampling from a parametric distributions.

We first build a model to work with throughout the tutorial.

import os
import site
site.addsitedir('D:\pycotools3')
from pycotools3 import model, tasks, viz
import pandas

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

Model(name=michaelis_menten, time_unit=s, volume_unit=l, quantity_unit=mol)

S = tasks.Scan(
    mm, scan_type='scan', subtask='time_course', report_type='time_course',
    report_name = 'ParameterScanOfTimeCourse.txt', variable='S',
    minimum=1, maximum=20, number_of_steps=8, run=True,
)

## Now check parameter scan data exists
os.path.isfile(S.report_name)

True

Two Way Parameter Scan

By default, scan tasks are removed before setting up a new scan. To set up dual scans, set clear_scans to False in a second call to Scan so that the first is not removed prior to adding the second.

## Clear scans for setting up first scan
tasks.Scan(
    mm, scan_type='scan', subtask='time_course', report_type='time_course',
    variable='E', minimum=1, maximum=20, number_of_steps=8, run=False, clear_scan=True,
)

## do not clear tasks when setting up the second
S = tasks.Scan(
    mm, scan_type='scan', subtask='time_course', report_type='time_course',
    report_name = 'TwoWayParameterScanOfTimeCourse.csv', variable='S',
    minimum=1, maximum=20, number_of_steps=8, run=True, clear_scan=False,
)

## check the output exists
os.path.isfile(S.report_name)

True

An arbitrary number of scans can be setup this way. Further, its possible to chain together scans with repeat or random distribution scans.

Repeat Scan Items

Repeat scans are very useful for running multiple parameter estimations and for running stochastic time courses.

## Assume Parameter Estimation task already configured
tasks.Scan(
    mm, scan_type='repeat', subtask='parameter_estimation', report_type='parameter_estimation',
    number_of_steps=6, run=False, ##set run to True to run via CopasiSE
)


## Assume model runs stochastically and time course settings are already configured
tasks.Scan(
    mm, scan_type='repeat', subtask='time_course', report_type='time_course',
    number_of_steps=100, run=False,  ##set run to True to run via CopasiSE
)

<pycotools3.tasks.Scan at 0x246af8b7be0>

Parameter Estimation Tutorial

import os, glob
import site
site.addsitedir(r'/home/ncw135/Documents/pycotools3')
site.addsitedir('D:\pycotools3')
from pycotools3 import viz, model, misc, tasks, models
from io import StringIO
import pandas
%matplotlib inline

Build a Model

working_directory = os.path.abspath('')

copasi_file = os.path.join(working_directory, 'negative_feedback.cps')

with model.BuildAntimony(copasi_file) as loader:
    mod = loader.load(
        """
        model negative_feedback
            compartment cell = 1.0
            var A in cell
            var B in cell

            vAProd = 0.1
            kADeg = 0.2
            kBProd = 0.3
            kBDeg = 0.4
            A = 0
            B = 0

            AProd: => A; cell*vAProd
            ADeg: A =>; cell*kADeg*A*B
            BProd: => B; cell*kBProd*A
            BDeg: B => ; cell*kBDeg*B
        end
        """
    )

## open model in copasi
#mod.open()
mod

Model(name=negative_feedback, time_unit=s, volume_unit=l, quantity_unit=mol)

Collect some experimental data

Organise your experimental data into delimited text files

experimental_data = StringIO(
    """
Time,A,B
 0, 0.000000, 0.000000
 1, 0.099932, 0.013181
 2, 0.199023, 0.046643
 3, 0.295526, 0.093275
 4, 0.387233, 0.147810
 5, 0.471935, 0.206160
 6, 0.547789, 0.265083
 7, 0.613554, 0.322023
 8, 0.668702, 0.375056
 9, 0.713393, 0.422852
10, 0.748359, 0.464639
    """.strip()
)

df = pandas.read_csv(experimental_data, index_col=0)

fname = os.path.join(os.path.abspath(''), 'experimental_data.csv')
df.to_csv(fname)

assert os.path.isfile(fname)

The Config Object

The interface to COPASI’s parameter estimation using pycotools3 revolves around the ParameterEstimation.Config object. ParameterEstimation.Config is a dictionary-like object which allows the user to define their parameter estimation problem. All features of COPASI’s parameter estimations task are supported, including configuration of validation experiments, affected experiments, affected validation experiments and constraints as well additional features such as the configuration of multiple models simultaneously via the affected_models keyword.

The ParameterEstimation.Config object expects at the bare minimum some information about the models being configured, some experimental data, some fit items and a working directory. The remaining options are automatically filled in with defaults.

config = tasks.ParameterEstimation.Config(
    models=dict(
        negative_feedback=dict(
            copasi_file=copasi_file
        )
    ),
    datasets=dict(
        experiments=dict(
            first_dataset=dict(
                filename=fname,
                separator=','
            )
        )
    ),
    items=dict(
        fit_items=dict(
            A={},
            B={},
        )
    ),
    settings=dict(
        working_directory=working_directory
    )
)
config

datasets:
    experiments:
        first_dataset:
            affected_models: negative_feedback
            filename: D:\pycotools3\docs\source\Tutorials\experimental_data.csv
            mappings:
                A:
                    model_object: A
                    object_type: Metabolite
                    role: dependent
                B:
                    model_object: B
                    object_type: Metabolite
                    role: dependent
                Time:
                    model_object: Time
                    role: time
            normalize_weights_per_experiment: true
            separator: ','
    validations: {}
items:
    fit_items:
        A:
            affected_experiments:
            - first_dataset
            affected_models:
            - negative_feedback
            affected_validation_experiments: []
            lower_bound: 1.0e-06
            start_value: model_value
            upper_bound: 1000000
        B:
            affected_experiments:
            - first_dataset
            affected_models:
            - negative_feedback
            affected_validation_experiments: []
            lower_bound: 1.0e-06
            start_value: model_value
            upper_bound: 1000000
models:
    negative_feedback:
        copasi_file: D:\pycotools3\docs\source\Tutorials\negative_feedback.cps
        model: Model(name=negative_feedback, time_unit=s, volume_unit=l, quantity_unit=mol)
settings:
    calculate_statistics: false
    config_filename: config.yml
    cooling_factor: 0.85
    copy_number: 1
    create_parameter_sets: false
    fit: 1
    iteration_limit: 50
    lower_bound: 1.0e-06
    max_active: 3
    method: genetic_algorithm
    number_of_generations: 200
    number_of_iterations: 100000
    overwrite_config_file: false
    pe_number: 1
    pf: 0.475
    population_size: 50
    prefix: null
    problem: 1
    quantity_type: concentration
    random_number_generator: 1
    randomize_start_values: false
    report_name: PEData.txt
    results_directory: ParameterEstimationData
    rho: 0.2
    run_mode: false
    save: false
    scale: 10
    seed: 0
    start_temperature: 1
    start_value: 0.1
    std_deviation: 1.0e-06
    swarm_size: 50
    tolerance: 1.0e-05
    update_model: false
    upper_bound: 1000000
    use_config_start_values: false
    validation_threshold: 5
    validation_weight: 1
    weight_method: mean_squared
    working_directory: D:\pycotools3\docs\source\Tutorials

The COPASI user will be familiar with most of these settings, though there are also a few additional options.

Once built, a ParameterEstimation.Config object can be passed to ParameterEstimation object.

PE = tasks.ParameterEstimation(config)

By default, the run_mode setting is set to False. To run the parameter estimation in background processes using CopasiSE, set run_mode to True or parallel.

config.settings.run_mode = True
PE = tasks.ParameterEstimation(config)
viz.Parse(PE)['negative_feedback']
# config

	A	B	RSS
0	0.000001	0.000001	7.955450e-12

Running multiple parameter estimations

With pycotools, parameter estimations are run via the scan task interface so that we have the option of running the same problem pe_number times. Additionally, pycotools provides a way of copying a model copy_number times so that the final number of parameter estimations that get executed is \(pe\_number\)cdot`copy_number``.

config.settings.copy_number = 4
config.settings.pe_number = 2
config.settings.run_mode = True
PE = tasks.ParameterEstimation(config)

And sure enough we have ran the problem 8 times.

viz.Parse(PE)['negative_feedback']

	A	B	RSS
0	0.000001	0.000001	7.955450e-12
1	0.000001	0.000001	7.955450e-12
2	0.000001	0.000001	7.955450e-12
3	0.000001	0.000001	7.955450e-12
4	0.000001	0.000001	7.955450e-12
5	0.000001	0.000001	7.955450e-12
6	0.000001	0.000001	7.955450e-12
7	0.000001	0.000001	7.955450e-12

A shortcut for configuring the ParameterEstimation.Config object

Manually configuring the ParameterEstimation.Config object can take some time as it is bulky, but necessarily so in order to enable users to configure any type of parameter estimation. The ParameterEstimation.Config class should be used directly when a lower level interface into COPASI configurations are required. For instance, if you want to configure different boundaries for each parameter, choose which parameters are affected by which experiment, mix timecourse and steasy state experiments, define independent variables, add constraints or choose which models are affected by which experiments, you can use the ParameterEstimation.Config class directly.

However, if you want a more standard configuration such as all parameters estimated between the same boundaries, all experiments affecting all parameters and models etc.. then you can use the ParameterEstimation.Context class to build the ParameterEstimation.Config class for you. The ParameterEstimation.Context class has a context argument that defaults to 's' for simple. While not yet implemented, everntually, alternative options for context will be provided to support other common patters, such as cross_validation or chaser_estimations (global followed by local algorithm). Note that an option is no longer required for model_selection since it is innately incorporated via the affected_models argument.

To use the ParameterEstimation.Context object

with tasks.ParameterEstimation.Context(mod, fname, context='s', parameters='g') as context:
    context.set('method', 'genetic_algorithm')
    context.set('population_size', 10)
    context.set('copy_number', 4)
    context.set('pe_number', 2)
    context.set('run_mode', True)
    context.set('separator', ',')
    config = context.get_config()

pe = tasks.ParameterEstimation(config)

The parameters keyword provides an easy interface for parameter selection. Here are the available options:

* `g` specifies that all global variables are to be estimated
* `l` specifies that all local parameters are to be estimated
* `m` specifies that all metabolites are to be estimated
* `c` specifies that all compartment volumes are to be estimated
* `a` specifies that all of the above will be estimated

These options can also be combined. For example, parameters='cgm' means that compartment volumes, global quantities and metabolite concentrations (or particle numbers) will be estimated.

viz.Parse(pe)['negative_feedback']

	RSS	kADeg	kBDeg	kBProd	vAProd
0	8.851340e-13	0.2	0.4	0.3	0.1
1	8.851340e-13	0.2	0.4	0.3	0.1
2	8.851340e-13	0.2	0.4	0.3	0.1
3	8.851340e-13	0.2	0.4	0.3	0.1
4	8.851340e-13	0.2	0.4	0.3	0.1
5	8.851340e-13	0.2	0.4	0.3	0.1
6	8.851340e-13	0.2	0.4	0.3	0.1
7	8.851340e-13	0.2	0.4	0.3	0.1

Examples

Simple Parameter Estimation

This is an example of how to configure a simple parameter estimation using pycotools. We first create a toy model for demonstration, then simulate some experimental data from it and fit it back to the model, using pycotools for configuration.

import os, glob
import pandas, numpy
import matplotlib.pyplot as plt
import seaborn
from pycotools3 import model, tasks, viz
seaborn.set_context(context='talk')

## Choose a directory for our model and analysis
working_directory = os.path.dirname(__file__)

## In this model, A gets reversibly converted to B but the backwards reaction is additionally regulated by C.
## B is reversibly converted into C.
antimony_string = """
model simple_parameter_estimation()
    compartment Cell = 1;

    A in Cell;
    B in Cell;
    C in Cell;

    // reactions
    R1: A => B ; Cell * k1 * A;
    R2: B => A ; Cell * k2 * B * C;
    R3: B => C ; Cell * k3 * B;
    R4: C => B ; Cell * k4 * C;

    // initial concentrations
    A = 100;
    B = 1;
    C = 1;

    // reaction parameters
    k1 = 0.1;
    k2 = 0.1;
    k3 = 0.1;
    k4 = 0.1;
end
"""

copasi_file = os.path.join(working_directory, 'example_model.cps')

## build model
with model.BuildAntimony(copasi_file) as builder:
    mod = builder.load(antimony_string)

assert isinstance(mod, model.Model)

## simulate some data, returns a pandas.DataFrame
data = mod.simulate(0, 20, 1)

## write data to file
experiment_filename = os.path.join(working_directory, 'experiment_data.txt')
data.to_csv(experiment_filename)

## We now have a model and some experimental data and can
## configure a parameter estimation

Parameter estimation configuration in pycotools3 revolves around the tasks.ParameterEstimation.Config object which is the input to the parameter estimation task. The object necessarily takes a lot of manual configuration to ensure it is flexible enough for any parameter estimation configuration. However, the ParameterEstimation.Context class is a tool for simplifying the construction of a Config object.

with tasks.ParameterEstimation.Context(mod, experiment_filename, context='s', parameters='g') as context:
    context.set('separator', ',')
    context.set('run_mode', True)
    context.set('randomize_start_values', True)
    context.set('method', 'genetic_algorithm')
    context.set('population_size', 100)
    context.set('lower_bound', 1e-1)
    context.set('upper_bound', 1e1)

    config = context.get_config()

pe = tasks.ParameterEstimation(config)

data = viz.Parse(pe).data
print(data)

Parameter estimation with multiple models

This is an example of how to configure a parameter estimation for multiple COPASI models using pycotools. We first create two similar but different toy models for demonstration, then simulate some experimental data from one of them and and fit it back to both models.

import os, glob
import pandas, numpy
import matplotlib.pyplot as plt
import seaborn
from pycotools3 import model, tasks, viz

seaborn.set_context(context='talk')

## Choose a directory for our model and analysis
working_directory = os.path.dirname(__file__)

model1_string = """
model model1()

    R1:   => A ; k1*S;
    R2: A =>   ; k2*A;
    R3:   => B ; k3*A;
    R4: B =>   ; k4*B*C; //feedback term
    R5:   => C ; k5*B;
    R6: C =>   ; k6*C;

    S = 1;
    k1 = 0.1;
    k2 = 0.1;
    k3 = 0.1;
    k4 = 0.1;
    k5 = 0.1;
    k6 = 0.1;
end
"""

model2_string = """
model model2()
    R1:   => A ; k1*S;
    R2: A =>   ; k2*A*C; //feedback term
    R3:   => B ; k3*A;
    R4: B =>   ; k4*B;
    R5:   => C ; k5*B;
    R6: C =>   ; k6*C;

    S = 1;
    k1 = 0.1;
    k2 = 0.1;
    k3 = 0.1;
    k4 = 0.1;
    k5 = 0.1;
    k6 = 0.1;
end
"""
copasi_file1 = os.path.join(working_directory, 'model1.cps')
copasi_file2 = os.path.join(working_directory, 'model2.cps')

antimony_strings = [model1_string, model2_string]
copasi_files = [copasi_file1, copasi_file2]

model_list = []
for i in range(len(copasi_files)):
    with model.BuildAntimony(copasi_files[i]) as builder:
        model_list.append(builder.load(antimony_strings[i]))

## simulate some data, returns a pandas.DataFrame
data = model_list[0].simulate(0, 20, 1)

## write data to file
experiment_filename = os.path.join(working_directory, 'data_from_model1.txt')
data.to_csv(experiment_filename)

with tasks.ParameterEstimation.Context(model_list, experiment_filename, context='s', parameters='g') as context:
    context.set('separator', ',')
    context.set('run_mode', True)
    context.set('randomize_start_values', True)
    context.set('method', 'genetic_algorithm')
    context.set('population_size', 25)
    context.set('lower_bound', 1e-1)
    context.set('upper_bound', 1e1)

    config = context.get_config()

pe = tasks.ParameterEstimation(config)

data = viz.Parse(pe).data

print(data)

API documentation

Here you will find detailed information about every module class and method in the pycotools3 package.

The model module

The pycotools3 model is of central importance in pycotools.

Model	Construct a pycotools3 model from a copasi file
ImportSBML	Import from sbml file
InsertParameters	Parse parameters into a copasi model
BuildAntimony	Build a copasi model using antimony
Build	Build a copasi model.

The tasks module

TimeCourse	Simulate a time course
ParameterEstimation.Config	A class for holding a parameter estimation configuration
ParameterEstimation	Interface to COPASI’s parameter estimation task
ParameterEstimation.Context	High level interface to create a ParameterEstimation.Config object.
Sensitivities	Interface to COPASI sensitivity task
Scan	Interface to COPASI scan task
Reports	Creates reports in copasi output specification section.

The viz module

The viz module exists to make visualising simulation output quick and easy for common patterns, such as plotting time courses or comparing parameter estimation output to experimental data. However it should be emphasised that the matplotlib and seaborn libraries are always close to hand in Python.

The viz module is currently in a state of rebuilding and so I only describe here the features which currently work.

Parse	General class for parsing copasi output into Python.
PlotTimeCourse	Plot time course data
Boxplots	Plot a boxplot for multi parameter estimation data.

Support

Users can post a question on stack-overflow using the pycotools tag. I get email notifications for these questions and will respond.

People

PyCoTools has been developed by Ciaran Welsh in Daryl Shanley’s lab at Newcastle University.

Caveats

• Non-ascii characters are minimally supported and can break PyCoTools
• Do not use unusual characters or naming systems (i.e. A reaction name called “A -> B” will break pycotools)
• In COPASI we can have (say) a global quantity and a metaboltie with the same name because they are different entities. This is not supported in Pycotools and you must use unique names for every model component

Citing PyCoTools

If you made use of PyCoTools, please cite this article using:

• Welsh, C.M., Fullard, N., Proctor, C.J., Martinez-Guimera, A., Isfort, R.J., Bascom, C.C., Tasseff, R., Przyborski, S.A. and Shanley, D.P., 2018. PyCoTools: a Python toolbox for COPASI. Bioinformatics, 34(21), pp.3702-3710.

And also please remember to cite COPASI:

• Hoops, S., Sahle, S., Gauges, R., Lee, C., Pahle, J., Simus, N., Singhal, M., Xu, L., Mendes, P. and Kummer, U., 2006. COPASI—a complex pathway simulator. Bioinformatics, 22(24), pp.3067-3074.

and tellurium:

• Medley, J.K., Choi, K., König, M., Smith, L., Gu, S., Hellerstein, J., Sealfon, S.C. and Sauro, H.M., 2018. Tellurium notebooks—An environment for reproducible dynamical modeling in systems biology. PLoS computational biology, 14(6), p.e1006220.