sys-bio / tellurium

Python Environment for Modeling and Simulating Biological Systems
http://tellurium.analogmachine.org/
Apache License 2.0
109 stars 36 forks source link

boundary and floating species #593

Open ahmedibatta opened 1 month ago

ahmedibatta commented 1 month ago

Hi, I have a question related to boundary and floating species. I have large model , I read the equations from a text file (which contains ODE equations with species (') and assignment rules(:=)) and an Excel sheet that has parameters and initial values of species. However, some species do not have ODE equations. After I read both files and built the model, I found that the model considers the species without ODEs as floating, and others as boundary species. The simulation of all is good and reasonable.

But after I omitted the species that don't have related ODE equations and added them as parameters, the model treated all species as floating but the simulation was not reasonable at all and toatlly different than previous one. Why is this behavior happening?

part of the whole model antimony format :

Full model in Antimony format:
// Created by libAntimony v2.14.0
// Compartments and Species:
compartment lung, blood;
species $I in lung, $M1 in lung, $pDC in lung, $AT2 in lung, $V in lung, $IFNb in lung;
// Assignment Rules:
virus_endocytosis := k_int*AT2*V*(1 - IFNb);
deg_I := b_I*I;
// Rate Rules:
I' = virus_endocytosis - deg_I - kill_CTL_I - k_V_innate*I*(M1 + pDC) - damage_ROS_AT2*I;
// Species initializations:
I = 0;
M1 = 1869.4;
// Compartment initializations:
lung = 1;
blood = 1;

// Variable initializations:
k_V_innate = 1e-11;
f_int = 1;
/ Other declarations:
var virus_endocytosis, deg_I,
const lung, blood, k_V_innate
matthiaskoenig commented 1 month ago

The most likely reason is that species can be either amounts or concentrations, depending on how they are defined. The hasOnlySubstrateUnits flag determines whether they are defined as amounts (substance_unit) or concentrations (substrate_unit/compartment_unit). All reactions are in substance_unit/time_unit. If you replace your concentrations with paramters, your equations won't match anymore, but you have to replace the use of S1 with p1/V1 in all equations, where p1 is the parameter that replaces S1 and V1 is the compartment S1 is in.

An easy way to make things correct is to just use units for all your objects, which will check all your equations for unit consistency.

TD/LR: You cannot just replace species with compartments without keeping track of what amounts and concentrations are and adjusting the ODEs accordingly. Use units when refactoring SBML models.

ahmedibatta commented 1 month ago

Thanks, @matthiaskoenig for your reply. I am currently building the model by initially ignoring the units. I’ve read the parameters and species with their values without adding units to any objects in the model. Previously, I tried a similar model, but without the intermediate rules. The equation was quite long, but it worked fine with almost the same code for reading the species and parameters. Now, I’ve separated the assignment rules from the ODE equation to investigate some intemediate steps in the reactions separately.

I do not think the problem in the code for read the equations itself,

def read_equations(file_path):
    with open(file_path, 'r') as file:
        lines = file.readlines()

    modified_lines = []
    for line in lines:
        line = line.strip()
        if '=' in line:
            variable, equation = line.split('=', 1)
            variable = variable.strip()
            equation = equation.strip()

            if not variable.endswith("'"):
                modified_lines.append(f"{variable} := {equation}") # to define them as assignment rules
            else:
                modified_lines.append(f"{variable} = {equation}")
        else:
            modified_lines.append(line)

    return '\n'.join(modified_lines)

# Then at the end : 
param_string = '\n'.join([f"{param} = {value};" for param, value in parameters.items()])
species_string = '\n'.join(species_string)
full_equations = f"{equations}\n\n{compartment_string}{species_string}\n\n{param_string}"
model = te.loada(full_equations)
luciansmith commented 1 month ago

If your compartment sizes aren't '1', and especially if your compartments change size in time, the difference in your two models is almost certainly what Matthias indicated: that when the symbol was a species, it was interpreted as a concentration, and when the symbol was a parameter, it was interpreted as an amount.

However, I can't be certain without an example. If you take your example you gave us define the unset parameters to '1', and perform a default simulation and plot, we get:

image

And if I change it to the following:

// Compartments and Species:
compartment lung, blood;
// Assignment Rules:
virus_endocytosis := k_int*AT2*V*(1 - IFNb);
deg_I := b_I*I;
// Rate Rules:
I' = virus_endocytosis - deg_I - kill_CTL_I - k_V_innate*I*(M1 + pDC) - damage_ROS_AT2*I;
// Species initializations:
I = 0;
M1 = 1869.4;
// Compartment initializations:
lung = 1;
blood = 1;

// Variable initializations:
k_V_innate = 1e-11;
f_int = 1;
// Other declarations:
var virus_endocytosis, deg_I
const lung, blood, k_V_innate

k_int = 1
b_I = 1
kill_CTL_I = 1
damage_ROS_AT2 = 1
pDC = 1
AT2 = 1
V  = 1
IFNb = 1

I get exactly the same plot:

image

However! If I set 'lung' to '2', I get two different plots, as expected:

image

image

(note the x axis.) This is because 'I' means two different things in those two different models.