Merge pull request #109 from sebapersson/Extend_SBML

Extend sbml support

Project.toml

@@ -1,7 +1,8 @@
 name = "PEtab"
 uuid = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
 authors = ["Viktor Hasselgren", "Sebastian Persson", "Rafael Arutjunjan", "Damiano Ognissanti"]
-version = "2.0.1"
+version = "2.0.2"
+
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -52,10 +53,10 @@ PEtabCatalystExtension = ["Catalyst"]
 PEtabFidesExtension = ["PyCall"]
 PEtabIpoptExtension = ["Ipopt"]
 PEtabOptimExtension = ["Optim"]
+PEtabPlotsExtension = ["Plots"]
 PEtabSciMLSensitivityExtension = ["SciMLSensitivity", "Zygote"]
 PEtabSelectExtension = ["PyCall"]
 PEtabZygoteExtension = ["SciMLSensitivity", "Zygote"]
-PEtabPlotsExtension = ["Plots"]
 
 [compat]
 CSV = "0.10"

src/SBML/Common.jl

@@ -51,10 +51,17 @@ function rewrite_piecewise_to_ifelse(rule_formula, variable, model_dict, base_fu
 
         if condition[1:2] == "lt" || condition[1:2] == "gt" || condition[1:2] == "eq" || condition[1:3] == "neq" || condition[1:3] == "geq" || condition[1:3] == "leq" 
             eq_syntax_dict[variable_change] = simple_piecewise_to_ifelse(condition, variable_change, value_active, value_inactive, model_dict, base_functions)
+
         elseif condition[1:3] == "and" || condition[1:2] == "if" || condition[1:2] == "or" || condition[1:3] == "xor" || condition[1:3] == "not"
             eq_syntax_dict[variable_change] = complex_piecewise_to_ifelse(condition, variable, value_active, value_inactive, model_dict, base_functions)
+
+        # Recursion to handle nested piecewise in condition             
+        elseif length(condition) ≥ 9 && condition[1:9] == "piecewise"
+            condition = rewrite_piecewise_to_ifelse(condition, "foo", model_dict, base_functions, model_SBML, ret_formula=true)
+            condition = rewrite_derivatives(condition, model_dict, base_functions, model_SBML)
+            eq_syntax_dict[variable_change] = simple_piecewise_to_ifelse(condition, variable_change, value_active, value_inactive, model_dict, base_functions)
         else
-            @error "Somehow we cannot process the piecewise expression"
+            @error "Somehow we cannot process the piecewise expression, condition = $condition"
         end
     end
 
@@ -124,6 +131,7 @@ end
 
 function simple_piecewise_to_ifelse(condition, variable, value_active, value_inactive, dicts, base_functions)
 
+    nested_condition::Bool = false
     if "leq" == condition[1:3]
         stripped_condition = condition[5:end-1]
         inEqUse = " <= "
@@ -146,13 +154,18 @@ function simple_piecewise_to_ifelse(condition, variable, value_active, value_ina
         return "true"
     elseif "false" == condition[1:5]
         return "false"
+    elseif "ifelse" == condition[1:6]
+        nested_condition = true
     else
         @error "Cannot recognize form of inequality, condition = $condition"
     end
 
-    parts = split_between(stripped_condition, ',')
-    # Trigger of event
-    expression = "ifelse(" * parts[1] * inEqUse * parts[2] * ", " * value_active * ", " * value_inactive * ")"
+    if nested_condition == false
+        parts = split_between(stripped_condition, ',')
+        expression = "ifelse(" * parts[1] * inEqUse * parts[2] * ", " * value_active * ", " * value_inactive * ")"
+    else
+        expression = "ifelse(" * condition * ", " * value_active * ", " * value_inactive * ")"
+    end
 
     return expression
 end

src/SBML/Process_rules.jl

@@ -18,6 +18,8 @@ function process_assignment_rule!(model_dict::Dict, rule_formula::String, rule_v
         return 
     end
 
+    rule_formula = replace_reactionid_with_math(rule_formula, model_SBML)
+
     #=
         If the rule does not involve a piecewise expression simply encode it as a function which downsteram
         is integrated into the equations when inserting "functions" into the reactions 
@@ -69,6 +71,8 @@ function process_rate_rule!(model_dict::Dict, rule_formula::String, rule_variabl
         end
     end
 
+    rule_formula = replace_reactionid_with_math(rule_formula, model_SBML)
+
     # Add rate rule as part of model derivatives and remove from parameters dict if rate rule variable
     # is a parameter
     if rule_variable ∈ keys(model_dict["parameters"])
@@ -153,7 +157,7 @@ function _time_dependent_ifelse_to_bool(formula_with_ifelse::String, model_dict:
         time_left = check_for_time(string(lhsRule))
         rewrite_ifelse = true
         if time_left == false && time_left == false
-            println("Have ifelse statements which does not contain time. Hence we do not rewrite as event, but rather keep it as an ifelse.")
+            @info "Have ifelse statements which does not contain time. Hence we do not rewrite as event, but rather keep it as an ifelse." maxlog=1
             rewrite_ifelse = false
             continue
         elseif time_left == true

src/SBML/SBML_to_ModellingToolkit.jl

@@ -48,6 +48,10 @@ function as_trigger(trigger_formula, model_dict, model_SBML)
         parts[1] = replace_whole_word(parts[1], "time", "t")
     end
 
+    # Account for potential reaction-math kinetics making out the trigger 
+    parts[1] = replace_reactionid_with_math(parts[1], model_SBML)
+    parts[2] = replace_reactionid_with_math(parts[2], model_SBML)
+
     # States in ODE-system are typically in substance units, but formulas in 
     # concentratio. Thus, each state is divided with its corresponding 
     # compartment 
@@ -104,7 +108,10 @@ function process_initial_assignment(model_SBML, model_dict::Dict, base_functions
     initally_assigned_parameter = Dict{String, String}()
     for (assignId, initialAssignment) in model_SBML.initial_assignments
 
+        # The formula might refer to reaciton, in this case insert reaction dynamics 
         _formula = SBML_math_to_str(initialAssignment)
+        _formula = replace_reactionid_with_math(_formula, model_SBML)
+
         formula = rewrite_derivatives(_formula, model_dict, base_functions, model_SBML)
         # Initial time i zero 
         formula = replace_whole_word(formula, "t", "0.0")
@@ -250,6 +257,13 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
         if model_dict["isBoundaryCondition"][state_id] == true
            model_dict["derivatives"][state_id] *= "0.0"
         end
+
+        # In case the conc. is given in initial conc, but the state should be in amounts this 
+        # must be acounted for with initial values
+        if model_dict["stateGivenInAmounts"][state_id][1] == false && model_dict["hasOnlySubstanceUnits"][state_id] == true
+            model_dict["stateGivenInAmounts"][state_id] = (true, state.compartment)
+            model_dict["states"][state_id] = string(state.initial_concentration) * " * " * state.compartment
+        end
     end
 
     # Extract model parameters and their default values. In case a parameter is non-constant 
@@ -306,6 +320,8 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
             event_assign_to[i] = event_assignment.variable
             event_formulas[i] = replace_function_with_formula(SBML_math_to_str(event_assignment.math), model_dict["modelFunctions"])
             event_formulas[i] = replace_whole_word(event_formulas[i], "t", "integrator.t")
+            event_formulas[i] = replace_reactionid_with_math(event_formulas[i], model_SBML)
+
             # Species typically given in substance units, but formulas in conc. Thus we must account for assignment 
             # formula being in conc., but we are changing something by amount 
             if event_assign_to[i] ∈ keys(model_SBML.species)
@@ -405,11 +421,19 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
             model_dict["derivatives"][product.species] *= "+" * stoichiometry * compartment_scaling * "(" * formula * ")"
         end
     end
-    # For states given in amount but model equations are in conc., multiply with compartment 
+    # For states given in amount but model equations are in conc., multiply with compartment, also handle potential 
+    # reaction identifayers in the derivative
     for (state_id, derivative) in model_dict["derivatives"]
+
+        model_dict["derivatives"][state_id] = replace_reactionid_with_math(model_dict["derivatives"][state_id], model_SBML)
+
         if model_dict["stateGivenInAmounts"][state_id][1] == false
             continue
         end
+        # Here equations should be given in amounts 
+        if model_dict["hasOnlySubstanceUnits"][state_id] == true
+            continue
+        end
         # Algebraic rule (see below)
         if replace(derivative, " " => "")[end] == '~' || replace(derivative, " " => "")[end] == '0'
             continue
@@ -420,7 +444,13 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
 
     # For states given by assignment rules 
     for (state, formula) in model_dict["assignmentRulesStates"]
-        model_dict["derivatives"][state] = state * " ~ " * formula
+        _formula = rewrite_derivatives(formula, model_dict, base_functions, model_SBML; check_scaling=true)
+        # Must track if species is given in amounts or conc.
+        if state ∈ keys(model_SBML.species) && model_SBML.species[state].only_substance_units == false
+            cmult = model_dict["stateGivenInAmounts"][state][1] == true ? " * " * model_SBML.species[state].compartment : ""
+            _formula = "(" * _formula * ")" * cmult
+        end
+        model_dict["derivatives"][state] = state * " ~ " * _formula
         if state ∈ non_constant_parameter_names
             delete!(model_dict["states"], state)
             delete!(model_dict["parameters"], state)
@@ -509,7 +539,17 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
         if divide_with_compartment == false
             continue
         end
-        model_dict["derivatives"][specie] = specie * " ~ (" * model_dict["states"][specie] * ") / " * c
+        if model_dict["stateGivenInAmounts"][specie][1] == true
+            model_dict["derivatives"][specie] = specie * " ~ (" * model_dict["states"][specie] * ") / " * c
+        else
+            # Must account for the fact that the compartment can change in size, and thus need to 
+            # account for its initial value 
+            if c ∈ keys(model_dict["parameters"])
+                model_dict["derivatives"][specie] = specie * " ~ (" * model_dict["states"][specie] * ")"
+            else
+                model_dict["derivatives"][specie] = specie * " ~ (" * model_dict["states"][specie] * ")" * " * " * string(model_dict["states"][c]) * " / " * c 
+            end
+        end
     end
 
     # In case the model has a conversion factor 
@@ -541,6 +581,35 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
         non_constant_parameter_names = filter(x -> x != id, non_constant_parameter_names)
     end
 
+    #=
+        Sometimes the volume might change over time but the amount should stay constant, as we have a boundary condition
+        In this case it follows that amount n (amount), V (compartment) and conc. are related 
+        via the chain rule by - I need to change my ODE:s and add state
+        dn/dt = d(n/V)/dt*V + n*dV/dt 
+    =#
+    for specie in keys(model_SBML.species)
+
+        if !(model_dict["stateGivenInAmounts"][specie][1] == true && 
+             model_SBML.species[specie].compartment ∈  rate_rules_names && 
+             model_dict["isBoundaryCondition"][specie] == true &&
+             specie ∈ rate_rules_names && 
+             model_dict["hasOnlySubstanceUnits"][specie] == false)
+            continue
+        end
+
+        # Derivative and inital values for concentratin species
+        compartment = model_SBML.species[specie].compartment
+        specie_conc = "__" * specie * "__conc__"
+        model_dict["states"][specie_conc] = model_dict["states"][specie] * " / " * compartment
+        i_start = findfirst(x -> x == '~', model_dict["derivatives"][specie]) + 1
+        i_end = findlast(x -> x == ')', model_dict["derivatives"][specie])
+        model_dict["derivatives"][specie_conc] = "D(" * specie_conc  * ") ~ " * model_dict["derivatives"][specie][i_start:i_end]
+
+        # Rebuild derivative for amount specie
+        itmp1, itmp2 = findfirst(x -> x == '~', model_dict["derivatives"][specie_conc]) + 1, findfirst(x -> x == '~', model_dict["derivatives"][compartment]) + 1
+        model_dict["derivatives"][specie] = "D(" * specie * ") ~ " * model_dict["derivatives"][specie_conc][itmp1:end] * "*" *  compartment * " + " * specie * "*" * model_dict["derivatives"][compartment][itmp2:end] * " / " * compartment
+    end
+
     model_dict["numOfParameters"] = string(length(keys(model_dict["parameters"])))
     model_dict["numOfSpecies"] = string(length(keys(model_dict["states"])))
     model_dict["non_constant_parameter_names"] = non_constant_parameter_names
@@ -676,11 +745,12 @@ function create_ode_model(model_dict, path_jl_file, model_name, write_to_file::B
             s_index += 1
         end
     end
-    for key in keys(model_dict["assignmentRulesStates"])
+    for (key, formula) in model_dict["assignmentRulesStates"]
+        what_write = key ∈ keys(model_dict["derivatives"]) ? model_dict["derivatives"][key] : key * " ~ " * model_dict["assignmentRulesStates"][key]
         if s_index != 1
-            dict_model_str["derivatives"] *= ",\n    " * key * " ~ " * model_dict["assignmentRulesStates"][key]
+            dict_model_str["derivatives"] *= ",\n    " * what_write
         else
-            dict_model_str["derivatives"] *= "    " * key * " ~ " * model_dict["assignmentRulesStates"][key]
+            dict_model_str["derivatives"] *= "    " * what_write
             s_index += 1
         end
     end
@@ -923,3 +993,12 @@ function _SBML_math_to_str(math::SBML.MathConst)
         return math.id, false
     end
 end
+
+
+function replace_reactionid_with_math(formula::T, model_SBML)::T where T<:AbstractString
+    for (reaction_id, reaction) in model_SBML.reactions
+        reaction_math = SBML_math_to_str(reaction.kinetic_math)
+        formula = replace_whole_word(formula, reaction_id, reaction_math)
+    end
+    return formula
+end