A simplified python interface to COPASI.

Hello, I've been trying to load and simulate the curated section of biomodels using basico. It seems there is a segfault when trying to simulate model 516 after 513 has already been loaded/simulated. Here's the code:

from basico import *
import os
import glob
import time
try:
    from roadrunner.tests.download_biomodels import download_biomodels
except ImportError:
    raise ImportError("Do a `pip install libroadrunner==2.2.0`")

biomodels_folder = os.path.join(os.path.dirname(__file__), "biomodels")

if not os.path.isdir(biomodels_folder):
    biomodels_folder = download_biomodels(os.path.dirname(__file__))

if not os.path.isdir(biomodels_folder):
    raise NotADirectoryError(biomodels_folder)

biomodels_files = glob.glob(os.path.join(biomodels_folder, "*.xml"))

print(len(biomodels_files))

start = time.time()
for model_file in [biomodels_files[512], biomodels_files[515]]:
    try:
        print(model_file)
        model = load_model(model_file)
        result = run_time_course(duration=100)
    except Exception:
        print(f"Model {model_file} failed")

end = time.time() - start

output:

1036
D:\roadrunner\roadrunner\docs\source\parallel\roadrunner_map\biomodels\BIOMD0000000513.xml
D:\roadrunner\roadrunner\docs\source\parallel\roadrunner_map\biomodels\BIOMD0000000516.xml

Process finished with exit code -1073740791 (0xC0000409)

Hello,

When doing parameter estimation, one loads data with the add_experiment() function which takes a parameter "file_name" to chose a specific storage location for the preprocessed file.

If i want my model and experiment files to be in a subfolder compared to the current directory i am running my script in, i would typically do:

from basico import *
import pandas as pd

model_io.load_model("copasi\\example.cps")
data = pd.read_csv("copasi\\example.txt", header=0, sep='\t')
add_experiment("example", data, file_name="copasi\\exp_example.tsv")

result = run_parameter_estimation()   # raises file not found (from C code ?)

The experiment saving works fine, but when run_parameter_estimation() looks for the file, it looks at a wrong path:

ERROR:root:Error while initializing the simulation: >ERROR 2022-01-14T10:00:09<
  CFitting (8): Failure reading file 'C:/Users/lbecquey/Desktop/minimal_example/copasi/copasi/exp_example.tsv'.

Notice the duplicated name of the subfolder in the path.

Note that using an absolute path everywhere does not change the problem.

Hello, I have a model which I can simulate without problem using the direct method:

from basico import *
load_model('mymodel.txt');
run_time_course(duration=1000.0, method='directMethod')

Next, I want to set stepsize=1000.0, to ensure that only the final datapoint is saved:

run_time_course(duration=1000.0, stepsize=1000.0, method='directMethod')

however, this gives an error:

ERROR:root:Error while running the simulation: >EXCEPTION 2022-02-22T11:21:58<
  CTrajectoryMethod (12): Internal step limit exceeded.

and an altered output (final timepoint is around t=775), see image:

This seems unexpected to me, as I only change when we save the output, I am uncertain why this causes an error.

Thie file is attached (originally had the .XML extension, but changed to .txt so that GitHub would accept it): mymodel.txt

Hi,

Following is the description of task (sensitivities) which has been extracted from a copasi model file. I want to retrieve the information enclosed in <problem> section.

<Task key="Task_16" name="Sensitivities" type="sensitivities" scheduled="false" updateModel="false"> <Report reference="Report_3" target="" append="1" confirmOverwrite="0"/> <Problem> <Parameter name="SubtaskType" type="unsignedInteger" value="1"/> <ParameterGroup name="TargetFunctions"> <Parameter name="SingleObject" type="cn" value="CN=Root,Model=MAPK cascade,Vector=Compartments[compartment],Vector=Metabolites[PP-MAPK],Reference=Concentration"/> <Parameter name="ObjectListType" type="unsignedInteger" value="1"/> </ParameterGroup> <ParameterGroup name="ListOfVariables"> <ParameterGroup name="Variables"> <Parameter name="SingleObject" type="cn" value="CN=Root,Model=MAPK cascade,Vector=Reactions[binding of MAPKKK activator],ParameterGroup=Parameters,Parameter=k1,Reference=Value"/> <Parameter name="ObjectListType" type="unsignedInteger" value="1"/> </ParameterGroup> <ParameterGroup name="Variables"> <Parameter name="SingleObject" type="cn" value=""/> <Parameter name="ObjectListType" type="unsignedInteger" value="0"/> </ParameterGroup> </ParameterGroup> </Problem> <Method name="Sensitivities Method" type="SensitivitiesMethod"> <Parameter name="Delta factor" type="unsignedFloat" value="0.001"/> <Parameter name="Delta minimum" type="unsignedFloat" value="9.9999999999999998e-13"/> </Method> </Task>

I used the function get_task_settings('Sensitivities', basic_only=False) with and without setting "basic_only" flag to False, but it gives the same following results:

{'scheduled': False, 'update_model': False, 'problem': {'SubtaskType': 1}, 'method': {'Delta factor': 0.001, 'Delta minimum': 1e-12, 'name': 'Sensitivities Method'}, 'report': {'filename': '', 'report_definition': 'Sensitivities', 'append': True, 'confirm_overwrite': True}}

Here you can see, the above function is returning only the name and value attributes of the sub-section (<Parameter >) within the <Problem> section.

Can you please suggest how can I retrieve all the required information within <Problem> section?

Thanks

Hi all, great library. I want to have timecourse results with SBML identifiers and be able to access either amount/concentration of species.

Currently there is the function

def __build_result_from_ts(time_series, use_concentrations=True):
    # type: (COPASI.CTimeSeries) -> pandas.DataFrame
    col_count = time_series.getNumVariables()
    row_count = time_series.getRecordedSteps()

    column_names = []
    column_keys = []
    for i in range(col_count):
        column_keys.append(time_series.getKey(i))
        column_names.append(time_series.getTitle(i))
        # FIXME: SBML ids for comparison
        column_names.append(time_series.getSBMLId(i))

    concentrations = np.empty([row_count, col_count])
    for i in range(row_count):
        for j in range(col_count):
            if use_concentrations:
                concentrations[i, j] = time_series.getConcentrationData(i, j)
            else:
                concentrations[i, j] = time_series.getData(i, j)

    df = pd.DataFrame(data=concentrations, columns=column_names)
    df = df.set_index('Time')

    return df

I tried to get the SBML identifiers but need a CDataModel for that in time_serices.getSBMLId(i). How would I do this? Is there some documention for the python-copasi API? Are there some examples and use cases?

How would I get amounts/concentrations for species?

Best Matthias

Based on the example in #17 I tried to set a custom selection to include the things I want to see in the results via

tc = run_time_course_with_output(['Time', '[S1]', '[S2]', 'Values[k1]'], 0, 5, 5)

But this results in the following error due to the Values[k1] in the selection.

--- c1 ---
ERROR:root:Error while initializing the simulation: >WARNING 2022-04-05T12:44:44<
  CCopasiTask (7): Problems compiling output.
Empty DataFrame
Columns: [Time, [S1], [S2], Values[k1]]
Index: []

Hi,

I want to retrieve complete information of methods used in any specific task. For example, copasi model with optimization task has the following method details:

<Method name="Genetic Algorithm" type="GeneticAlgorithm"> <Parameter name="Log Verbosity" type="unsignedInteger" value="0"/> <Parameter name="Number of Generations" type="unsignedInteger" value="200"/> <Parameter name="Population Size" type="unsignedInteger" value="20"/> <Parameter name="Random Number Generator" type="unsignedInteger" value="1"/> <Parameter name="Seed" type="unsignedInteger" value="0"/> <Parameter name="Mutation Variance" type="float" value="0.10000000000000001"/> <Parameter name="Stop after # Stalled Generations" type="unsignedInteger" value="0"/> </Method>

I want to retrieve/set all the information including the "type" of the method. I tried to use get_opt_settings() but i returns on the following information related to method: 'method': {'Number of Generations': 200, 'Population Size': 20, 'name': 'Genetic Algorithm'}

I tried this code:

add_species('B');
set_species('B', exact=True,  initial_concentration=0.4, status='fixed');

and the result is that species 'B' ends up with initial concentration of 1 (default value, set by the add_species), not 0.4 as should have been changed by set_species ...

I had to write:

add_species('B', initial_concentration=0.4);
set_species('B', exact=True, status='fixed');

to have 'B' actually get the initial concentration of 0.4.

Hi @fbergmann ,

Square brackets are missing in the columns names of sim_data returned by the function get_simulation_results for time course data.

e.g.

sim_data>
Index(['ADP', 'fructose 6-phosphate', 'NAD', 'glyceraldehyde 3-phosphate',
       'pyruvate{cytosol}', 'dihydroxyacetone phosphate', '2-phosphoglycerate',
       dtype='object')

exp_data>
Index(['Time', '[d-glucose{"external compartment (blood)"}]_0', '[ATP]',
       '[dihydroxyacetone phosphate]', '[fructose 6-phosphate]',
       dtype='object')

For steady-state data both sim_data and exp_data have brackets for the concentration of species in the column names.

Could you please have a look?

Hello basiCO team,

when I'm running the function 'plot_per_dependent_variable()' I get the following error:

'AttributeError: 'CDataObject' object has no attribute 'setDblValue'' with the following traceback:

(Environment: Anaconda 3; Python 3.8.3 (default, Jul 2 2020, 17:30:36) [MSC v.1916 64 bit (AMD64)]; IPython 7.16.1)

AttributeError                            Traceback (most recent call last)
<ipython-input-12-0f3cf678761d> in <module>
----> 1 plot_per_dependent_variable();

C:\anaconda3\lib\site-packages\basico\task_parameterestimation.py in plot_per_dependent_variable(**kwargs)
   1176         return result
   1177
-> 1178     exp_data, sim_data = get_simulation_results(**kwargs)
   1179
   1180     dependent_variables = {}

C:\anaconda3\lib\site-packages\basico\task_parameterestimation.py in get_simulation_results(values_only, **kwargs)
   1027             model.updateInitialValues(change_set)
   1028
-> 1029         _update_fit_parameters_from(dm, solution, exp_name)
   1030
   1031         if experiment.getExperimentType() == COPASI.CTaskEnum.Task_steadyState:

C:\anaconda3\lib\site-packages\basico\task_parameterestimation.py in _update_fit_parameters_from(dm, solution, exp_name)
   1078                 obj.getObjectParent().setInitialConcentration(value)
   1079             elif type(obj.getObjectParent()) is COPASI.CCopasiParameter:
-> 1080                 obj.setDblValue(value)
   1081                 model.updateInitialValues(obj)
   1082             else:

C:\anaconda3\lib\site-packages\COPASI\__init__.py in <lambda>(self, name)
   3235     for _s in [CObjectInterface, CUndoObjectInterface]:
   3236         __swig_getmethods__.update(getattr(_s, '__swig_getmethods__', {}))
-> 3237     __getattr__ = lambda self, name: _swig_getattr(self, CDataObject, name)
   3238     __repr__ = _swig_repr
   3239     Container = _COPASI.CDataObject_Container

C:\anaconda3\lib\site-packages\COPASI\__init__.py in _swig_getattr(self, class_type, name)
     83     if method:
     84         return method(self)
---> 85     raise AttributeError("'%s' object has no attribute '%s'" % (class_type.__name__, name))
     86
     87

I'm using the example model 'LM' and I'm following the steps of the parameter estimation as listed in the documentation here: https://basico.readthedocs.io/en/latest/notebooks/Parameter_Estimation_Example.html Input lines 1 - 11 all work without any issue.

When changing a function on a reaction, it should automatically add modifiers to the reaction scheme.

• currently, when a reaction r1 with scheme A -> B is assigned a kinetic using modifiers this fails

In the gui, this is allowed, however, this results in an invalid function, since no mapping is specified at that point. I would rather want to avoid creating invalid states, as it would be hard to track down later.

So i propose to silently add modifiers to the reaction scheme, iff a mapping for modifier to species is defined. that means:

• set_reaction('r1', function='function with modifier') would still fail, however
• set_reaction('r1', function='function with modifier', mapping={'formulaTokenUsedAsModifier', 'speciesName' }) would work

@pmendes would that work for you?

Hello @fbergmann ,

I tried to make a few changes in the function def plot_per_dependent_variable to create plots for the steady-state data and I would like to share the modified function (please see below) with you.

Thanks so much for the new release, this helps me a lot in obtaining the steady-state data in the required format.

def plot_per_dependent_variable_ss(**kwargs):
    """
    This function creates a figure for each dependent variable, with traces for all experiments.

    :param kwargs:

    - | `model`: to specify the data model to be used (if not specified
      | the one from :func:`.get_current_model` will be taken)

    :return: array of tuples (fig, ax) for each figure created
    """
    dm = model_io.get_model_from_dict_or_default(kwargs)

    task = dm.getTask(TASK_PARAMETER_ESTIMATION)
    assert (isinstance(task, COPASI.CFitTask))

    problem = task.getProblem()
    assert (isinstance(problem, COPASI.CFitProblem))

    experiments = problem.getExperimentSet()
    assert (isinstance(experiments, COPASI.CExperimentSet))

    result = []
    num_experiments = experiments.getExperimentCount()
    if num_experiments == 0:
        return result

    data = get_simulation_results(**kwargs)
    exp_data = data[0]  # experimental data
    sim_data = data[1]  # simulation data

    dependent_variables = {}

    for i in range(num_experiments):
        experiment = experiments.getExperiment(i)
        mapping = get_experiment_mapping(experiment)

        # set independent values for that experiment
        independent = mapping[mapping.type == 'independent']
        xname = independent.iloc[0].mapping

        # set dependent values for that experiment
        dependent = mapping[mapping.type == 'dependent']
        num_dependent = dependent.shape[0]
        for j in range(num_dependent):
            name = dependent.iloc[j].mapping
            if name not in dependent_variables:
                dependent_variables[name] = []
            dependent_variables[name].append(i)

    for dependent in dependent_variables:
        fig, ax = plt.subplots()
        cycler = plt.cycler("color", plt.cm.tab20c.colors)()
        ax.set_title(dependent)
        experiment_indices = dependent_variables[dependent]

        for i in experiment_indices:
            experiment = experiments.getExperiment(i)
            exp_name = experiment.getObjectName()
            nextval = next(cycler)['color']
            name = dependent
            if name not in sim_data[i].columns:
                name = name[1:-1]

            sim_data[i][xname] = exp_data[i][xname]
            # print(sim_data[i])

            sim_data[i].reset_index().plot(x=xname, y=name,
                                           label="{0} Fit".format(exp_name), ax=ax, color=nextval)
            exp_data[i].plot.scatter(x=xname, y=dependent, ax=ax, color=nextval,
                                     label='{0} Measured'.format(exp_name))
        # plt.savefig('ss.png')
        result.append((fig, ax))

    return result