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:
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?
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.
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
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);
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.
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
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.
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.
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.
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
end
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?
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.
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.