Merging Kinetic Models

[Andrew Glanville]

Hi,

I'm attempting to merge the following -

BIOMD0000000009_url.xml

BIOMD0000000426_url.xml

BIOMD0000000093_url.xml

(forgive me I'm a novice [obviously] and this is more a computer exercise than a biological one) I spoke with Dr Lobo and he advised these are kinetic models / reason why mergem doesn't work etc. I tried adding back rate laws (bits) after mergem but didn't work.

Evidently COPASI can merge the three using the UI by simply Add to model. My question= how do I "add to model" using the bindings? Presumably there is a way of doing this programatically?

Many thanks in advance,

Andy

[Frank Bergmann]

We have not yet simplified the API for this, so it is not yet in basico or corc. However, I've just tried it and it does work. As long as you are careful, and merge always objects of the same type (i.e compartments with compartments, species with species and so on) The following uses  a mix of basico and the normal COPASI python API.

1. first download all the biomodels and save them locally as copasi files: 

from basico import *
import COPASI

dms = {}
for m in ['BIOMD0000000009', 'BIOMD0000000426', 'BIOMD0000000093']:
    dms[m] = load_biomodel(m)
    save_model(m + '.cps', model=dms[m])

2. choose one of those models to start out with, or choose a new model to add all of the other models to:

nm = new_model(name='Combined')
for m in ['BIOMD0000000009', 'BIOMD0000000426', 'BIOMD0000000093']:
    nm.addModel(m + '.cps')

3. now you have all the elements of the three models together. For example you see that there are 4 compartments now when running: 

get_compartments(model=nm)

4. what you need to do from then on is to go through all items that you identified as being the same and merge them together. So lets merge the compartments 'compartment_[merge]' and 'compartment[merge]': 

model = nm.getModel()
expa = COPASI.CModelExpansion(model)
emap = COPASI.CModelExpansion_ElementsMap()

c1 = model.getCompartment('compartment[merge]')
c2 = model.getCompartment('compartment_[merge]')

emap.add(c1, c2)
expa.replaceInModel(emap, True)

5. now look again, and you will find only 3 compartments in the model. 

Let us know if you have further questions. 

Frank

[Andrew Glanville]

Many thanks Frank.

I'll give it a try!

[Andrew Glanville]

Yes, it works!

[Andrew Glanville]

Hi again,

Another question= Taking the resultant merged model and using basico to run a timecourse gives a result the same as the last model added.

Presumably looks only at the last compartment?

Is there an easy way of getting a more merged looking output? (get the compartments to interact somehow) please see image and code-

from basico import *
from datetime import datetime
import sys
#import pandas as pd
import matplotlib.pyplot as plt
#%matplotlib inline

# Examine how many functions are non-reversible or reversible
#get_functions(reversible=False)
#print ("Non-Reversible Functions")
#print ("-------------------------")
#print (get_functions(reversible=False))
#print ("")
#print ("Reversible Functions")
#print ("-------------------------")
#get_functions(reversible=True)
#print (get_functions(reversible=True))


load_model('combined.cps');
#load_model('BIOMD0000000009_url.xml'); # MAPK model
#load_model('BIOMD0000000426_url.xml'); # PI3K_AKT_mTOR model
#load_model('BIOMD0000000093_url.xml'); # JAK/STAT model

# Here's the time course - Some functions are reversible hence kinetic function Reversible Hill
data = run_time_course(duration=10, method='Reversible Hill', use_numbers=True, use_seed=False);

# All data is sent to pandas for the chart
ax = data.plot(lw=2, colormap='jet', marker='.', markersize=1, title='Three Merged Models combined.cps');
#ax = data.plot(lw=2, colormap='jet', marker='.', markersize=1, title='BIOMD0000000093_url.xml');
ax.set_xlabel("Time");
ax.set_ylabel("Concentrations");

print('there are {0} reports in the model'.format(len(get_reports())))
get_reports(ignore_automatic=True)
get_reports(task=T.STEADY_STATE)
assign_report('Time-Course', task=T.TIME_COURSE, filename='out.txt', append=True)

set_task_settings(T.TIME_COURSE,
   settings = {'report':
              {'report_definition': 'X Time-Course'}})

get_task_settings(T.TIME_COURSE)['report']['report_definition']

# record the plot by making a *.png file with the date and time as name
grafik = datetime.now().strftime("%Y_%m_%d-%I_%M_%S_%p")
print (grafik)
plt.savefig(grafik + '.png')

[Andrew Glanville]

Sorry, here's some more output

[Frank Bergmann]

If you simply merge the models, i would expect the results to remain quite similar until you add additional reactions. So I'm not quite sure what you mean. Copasi will simulate the system of ODEs for all compartments. 

Just one word for the run_time_course function. You call it with: 

    data = run_time_course(duration=10, method='Reversible Hill', use_numbers=True, use_seed=False); 

Here the method refers to the implemented integration method / stochastic method with which to simulate the model. Valid arguments are: 

   * deterministic / lsoda: the LSODA implementation

   * stochastic: the Gibson & Bruck Gillespie implementation

   * directMethod: Gillespie Direct Method

   * others: hybrid, hybridode45, hybridlsoda, adaptivesa, tauleap, radau5, sde

You can also pass along the model argument when running the time course, if you think that the wrong model might be simulated. 

best
Frank

[Andrew Glanville]

Thanks Frank,

Yes, apologies for my noviceness, I thought what I was seeing was limited to one compartment but I started altering some values as attached.

Thanks again for your help.

Andy

--
You received this message because you are subscribed to a topic in the Google Groups "COPASI User Forum" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/copasi-user-forum/MMvEtciQI0o/unsubscribe.
To unsubscribe from this group and all its topics, send an email to copasi-user-fo...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/copasi-user-forum/49a0a5d2-d89b-4509-bdd8-ca82d6b37711n%40googlegroups.com.

[Andrew Glanville]

Another question. I'm "assuming" if I poke one of my component models with a species stick (Pi in this instance) and the reaction is different (in some way for the given species) in the combined model I have found crosstalk?

See image - Or is this more likely to be some sort of unit missmatch i.e. 0.8 / 8 - Or is this something else?