Performing Global Sensitivity Analysis on a big reaction network imported as a SBML file

Hi everyone,

I have been using Julia language for the last couple of months to simulate a big biochemical reaction network, which is generated by using BioNetGen (around 500 species and 170 parameters but the numbers can go way up) in an SBML format and I use SBMLToolkit to import the model. Most of my equations obey the law of mass action but I have a few custom Hill-type functions for some reactions. My main goal is to perform Global Sensitivity Analysis as a first step to understanding my model but I face a couple of problems and questions:

  1. Some of the free parameters are the initial conditions for the species, but I observed that changing the parameter, does not change the initial condition automatically. That can be fixed by hand when we do not have many initial species but what about when the initial number of species is large?

  2. I observed that ODAEProblem can be faster than ODEProblem but I am not sure if it is the right approach. One issue with that is that using ODAEProblem on models imported from an SBML file does not provide the possibility to use solvers for DAE problems.

  3. Moreover, and most important, it is common to get these 2 errors when I increase the value of the parameters of the model more than an order of magnitude:

a) Warning: dt <= dtmin. Aborting. There is either an error inyour model specification or the true solution is unstable

b) Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

Here is a sample of my code:

using DifferentialEquations
using OrdinaryDiffEq
using Catalyst
using SBMLToolkit
using GlobalSensitivity, QuasiMonteCarlo
using Plots
using Sundials, LSODA


mdl = readSBML("file_sbml.xml", doc -> begin
    set_level_and_version(3, 2)(doc)

rn = ReactionSystem(mdl);
odesys = convert(ODESystem, rn);

T=1.0; #sec
tspan = (0.0, T);
prob_ode = ODEProblem(odesys, Float64[], tspan, Float64[]);

n=length(prob_ode.p)-1; #the last parameter is "cell". equal to 1.0 here. 


bounds=collect(map(collect, zip(lb,ub)))

 f1 = function(p)
     prob1 = remake(prob_ode;p=p)
     sol = solve(prob1,alg_hints=[:auto],save_everystep=false)


sampler = SobolSample(); #choose sampler for Sobol sensitivity analysis
A,B = QuasiMonteCarlo.generate_design_matrices(n_samples,lb,ub,sampler);

sobol_result = gsa(f1,Sobol(),A,B);

But it transforms it into an ODEProblem, which is generally faster to solve. So yes, you cannot use DAE solvers on that anymore, but… why would you?

Can you make an MWE for that? That shouldn’t be the case.

Double check the parameters where that’s happening. Most likely the model isn’t stable at those parameters. You need to define a way that you want to reject said parameters by handling the retcode.

You can modify your f1 function such that the remake takes the indexes inp that correspond to the initial condition. Assuming all your u0s are the first 5 elements in p you would change the remake call to prob1 = remake(prob_ode;u0 = p[1:5], p = p[5:end]).

Thank you for your answer! One problem with the SBMLToolkit is that does not keep the same order for the parameters and the species as the original file (in fact it is random). So it is possible but you have to identify the order and modify the f1 function. It can be feasible in my case but in larger models with many initial conditions it might be difficult and of course, there is always the case to miss something if you do it by hand. I do not know if there is an option in SBMLToolkit to keep the order the same as the original file.

Thanks, Chris for your reply! I was able to identify the parameters that were causing the code to stop.

Since I am new to Julia, would you suggest using ODAEProblem or ODEProblem? Which are the cases to use each one appropriately?

For modifying the parameters, what happens if you use a parameter mapping and repass the initial condition mapping? i.e. something like

getval(v) = getmetadata(v, Symbolics.VariableDefaultValue)
pmap = parameters(odesys) .=> getval.(parameters(odesys))
u0map = states(odesys) .=> getval.(states(odesys))

 f1 = function(p)
     for (idx, param) in enumerate(p)
          pmap[idx] = pmap[idx][1] => param
     prob1 = remake(prob_ode; p = pmap, u0 = u0map)
     sol = solve(prob1,alg_hints=[:auto],save_everystep=false)

One caveat, if ODAEProblem is reducing the number of parameters or states then I’m not sure if what I proposed will still work (since it is using the original system’s full set of states and parameters, and their orderings). You’d have to determine the reduced mappings, in the correct order, that ODAEProblem needs, but I’m not sure how one does that.

It’s very problem dependent and depends on the branching structure of the problem. I’d say, usually just try ODEProblem first.

Thank you for your response Isaac! I tried what you suggested and I get this error:

KeyError: key Symbolics.VariableDefaultValue not found

I am not sure if I need a specific library or if I do something wrong. About the


you are right. I checked it and the order of the parameters and species stays the same but it does not
include the observables and the functions defined in the SBML file as species. On the other hand,


includes them as species and as a result they exist in the initial conditions.

It can eliminate them, though they are still accessible in the solution like sol[extra_variable]), it will just create the extra piece on demand using observed equations.

You need to install and then load the ModelingToolkit and Symbolics packages for the code I provided to work:

using ModelingToolkit, Symbolics

I load them but I still get the same error. I use Julia 1.8.3, v8.41.0 for ModelingToolkit, and v4.14.0 for Symbolics.

Ahh, I guess SBMLToolkit doesn’t use metadata to store the default values.

Another option that might work is if you define getval like

defs = ModelingToolkit.defaults(odesys)
getval(v) = defs[v]

That way does not give any error but still, the initial conditions do not change, they stay equal to 10.0 as I defined them initially. So, the pmap changes but the u0map stay the same as the initial one.

@Nick1 can you check whether in ModelingToolkit.defaults(odesys) the initial conditions that should change are still set to the symbolic parameters, or whether numerical values have already been substituted by SBMLToolkit? If it is substituting in numeric values that would prevent remake from being able to update the initial condition…

The example below works fine for me so I think this may be a problem with SBMLToolkit putting in a numerical value as the initial condition and dropping the symbolic value (which would break remake):

using Catalyst, DifferentialEquations, ModelingToolkit
@parameters β ν S0
@variables t S(t) I(t) R(t)
rx1 = Reaction(β, [S,I], [I], [1,1], [2])
rx2 = Reaction(ν, [I], [R])
defs = Dict([β => .0001, ν => .01, S0 => 999.0, S => S0, I => 1.0, R => 0.0])
@named sir = ReactionSystem([rx1,rx2],t,[S,I,R], [β,ν,S0]; defaults=defs)
odesys = convert(ODESystem, sir)
oprob = ODEProblem(odesys, Float64[], (0.0, 1000.0))
@assert  oprob.u0[1] == 999.0

getval(v) = defs[v]
pmap = parameters(odesys) .=> getval.(parameters(odesys))
u0map = states(odesys) .=> getval.(states(odesys))

# change S0 to 10.0
p = getval.(parameters(odesys))
p[end] = 10.0 

for (idx, param) in enumerate(p)
    pmap[idx] = pmap[idx][1] => param
oprob2 = remake(oprob; p=pmap, u0=u0map)
@assert oprob2.u0[1] == 10.0

They have been substituted in numeric values. I doubled checked, the initial conditions do not change. Is there any way for the SBMLToolkit to keep the symbolic value?

I opened an issue in the SBMLToolkit repo, Symbolic formulas for initial conditions being dropped · Issue #106 · SciML/SBMLToolkit.jl · GitHub, let’s see what @anandj or @paulflang says.

On Slack, @paulflang said that in

mdl = readSBML("my_model.xml", doc -> begin
    set_level_and_version(3, 2)(doc)

the call to convert_simplify_math(doc) is what is substituting out the symbolic expressions, and you should instead do:

const sbml_promote_expand = SBML.libsbml_convert(
    ["promoteLocalParameters", "expandFunctionDefinitions"] .=> Ref(Dict{String,String}()),

mdl = readSBML("my_model.xml", doc -> begin
                     SBML.set_level_and_version(3, 2)(doc)

If that doesn’t work it would probably be best for you to followup on the Julia Slack sciml-sysbio channel, as the SBMLToolkit developers look at messages there and can hopefully help you with further issues.

Thanks for your help @isaacsas!

Let us know if that solves the issue, or else please do follow up on Slack. I think stuff like this is useful feedback for improving SBMLToolkit’s interface.