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

**URL:** <https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000>\
**Category:** Modelling & Simulations\
**Created:** [January 15, 2023, 8:31pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000 "2023-01-15T20:31:25Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [January 15, 2023, 8:31pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/1 "2023-01-15T20:31:25Z")

</div>

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 ` **in** `your 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:

```julia

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

cd(raw"/pathway")

SBMLToolkit.checksupport_file("file_sbml.xml");
mdl = readSBML("file_sbml.xml", doc -> begin
    set_level_and_version(3, 2)(doc)
    convert_simplify_math(doc)
        end);

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. 
lb=zeros(n,);
ub=100.0*ones(n,);

push!(ub,1.0); 
push!(lb,1.0); 

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)
     [sol[177,2]] 
 end

n_samples=100;

#Sobol
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);

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 16, 2023, 3:47am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/2 "2023-01-16T03:47:45Z")

</div>

> [@Nick1](#):
>
> 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.

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?

> [@Nick1](#):
>
> - 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?

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

> [@Nick1](#):
>
> 1. 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 ` **in** `your model specification or the` true` solution is unstable`

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.

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [January 16, 2023, 4:26am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/3 "2023-01-16T04:26:39Z")

</div>

> [@Nick1](#):
>
> 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?

You can modify your `f1` function such that the `remake` takes the indexes in`p` that correspond to the initial condition. Assuming all your `u0`s 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])`.

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [January 18, 2023, 5:50pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/4 "2023-01-18T17:50:12Z")

</div>

> [@Vaibhavdixit02](#):
>
> You can modify your `f1` function such that the `remake` takes the indexes in`p` that correspond to the initial condition. Assuming all your `u0`s 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.

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [January 18, 2023, 5:52pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/5 "2023-01-18T17:52:43Z")

</div>

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

> [@Nick1](#):
>
> 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.

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

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [January 19, 2023, 1:17am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/6 "2023-01-19T01:17:35Z")

</div>

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

```julia
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
     end
     prob1 = remake(prob_ode; p = pmap, u0 = u0map)
     sol = solve(prob1,alg_hints=[:auto],save_everystep=false)
     [sol[177,2]] 
 end

```

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 19, 2023, 2:38am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/7 "2023-01-19T02:38:59Z")

</div>

> [@Nick1](#):
>
> Since I am new to Julia, would you suggest using ODAEProblem or ODEProblem? Which are the cases to use each one appropriately?

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

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [January 20, 2023, 7:21pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/8 "2023-01-20T19:21:02Z")

</div>

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

> ODAEProblem

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,

> ODEProblem

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 20, 2023, 7:33pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/9 "2023-01-20T19:33:55Z")

</div>

> [@Nick1](#):
>
> 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.

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.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [January 20, 2023, 7:46pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/10 "2023-01-20T19:46:25Z")

</div>

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

```julia
using ModelingToolkit, Symbolics

```

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [January 20, 2023, 7:51pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/11 "2023-01-20T19:51:19Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [January 21, 2023, 12:11am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/12 "2023-01-21T00:11:57Z")

</div>

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

Another option that might work is if you define `getval` like

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

```

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [January 21, 2023, 1:43am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/13 "2023-01-21T01:43:59Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [January 21, 2023, 1:51am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/14 "2023-01-21T01:51:55Z")

</div>

@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…

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [January 21, 2023, 2:07am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/15 "2023-01-21T02:07:46Z")

</div>

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`):

```julia
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
end
oprob2 = remake(oprob; p=pmap, u0=u0map)
@assert oprob2.u0[1] == 10.0

```

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [January 21, 2023, 2:48am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/16 "2023-01-21T02:48:24Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [January 21, 2023, 2:59am UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/17 "2023-01-21T02:59:12Z")

</div>

I opened an issue in the SBMLToolkit repo, [Symbolic formulas for initial conditions being dropped · Issue #106 · SciML/SBMLToolkit.jl · GitHub](https://github.com/SciML/SBMLToolkit.jl/issues/106), let’s see what @anandj or @paulflang says.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [January 22, 2023, 4:33pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/18 "2023-01-22T16:33:38Z")

</div>

On Slack, @paulflang said that in

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

```

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

```julia
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)
                     sbml_promote_expand(doc)
                 end)

```

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.

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [January 22, 2023, 8:23pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/19 "2023-01-22T20:23:26Z")

</div>

Thanks for your help @isaacsas!

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [January 22, 2023, 8:35pm UTC](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000/20 "2023-01-22T20:35:21Z")

</div>

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.

[Next page](https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000.md?page=2)