Parameter Estimation

import os, glob
import site
from pycotools3 import viz, model, misc, tasks
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')
ant = """
        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
        """
mod = model.loada(ant, copasi_file)

## 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
    context: s
    cooling_factor: 0.85
    copy_number: 1
    create_parameter_sets: false
    cross_validation_depth: 1
    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
    pl_lower_bound: 1000
    pl_upper_bound: 1000
    population_size: 50
    prefix: null
    problem: Problem1
    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.955430e-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)

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

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.