Issue with Local Sensitivites and dt <= dtmin

Specific Domains Numerics
,
1

Hi,

I am trying to calculate the local sensitivities for a large circulation model (Shi 2006)

I am Importing the model using CellML.jl when I define an ODEproblem and simulate the solution appears fine and also when performing GSA the results are also sound.

Using the code

using ForwardDiff

savetime = LinRange(18, 19, 100)

p = prob.p 
u0 =prob.u0


function circ_local(x)
    newprob = remake(prob; p=x[1:end])
    sol = solve(newprob, Tsit5(); reltol=1e-6, abstol=1e-12,saveat = savetime)
    [mean(sol[@nonamespace sys.LV₊V]), mean(sol[@nonamespace sys.Pat₊Pi])]
end
S = ForwardDiff.jacobian(circ_local,p)

Note: all the solver opts I use are the same for both simulation and performing GSA.

The above returns the following warnings with the NAN’s

Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/ZBye7/src/initdt.jl:121
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/Vg9hW/src/integrator_interface.jl:345
2×154 Matrix{Float64}:
 NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  …  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN
 NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN     NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN

Any help would be very much appreciated :slight_smile:

Cheers,
Harry

2

have you tried using a stiff solver?

3

Indeed, trying a method for stiff equations is the first thing to do in this case. Did you try it?

Follow the guide:

4

Hi @ChrisRackauckas & @Oscar_Smith

Tried a stiff solver and all the other solvers that I could find with increasingly high tolerances too.

image

The Image above is the graph of the solution at the savetime.

I hate to sound exactly like the meme in the guide above just I am cautious to claim model error seen as the simulation and function for GSA works well plus I am loading the model in from CellML.

Also intersting when I try and run the model for Waveform data

function circ_local(params)
    newprob = remake(prob;  p=params[1:end])
    newsol = solve(newprob, KenCarp4(); reltol=1e-12, abstol=1e-12, saveat=savetime)
    [newsol[@nonamespace sys.LV₊V]; newsol[@nonamespace sys.Sas₊Pi] ]
end
S = ForwardDiff.jacobian(circ_local,p)

The outcome is a 0x154 and errors in the same fashion as the above

┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/ZBye7/src/initdt.jl:121
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/Vg9hW/src/integrator_interface.jl:345
0×154 Matrix{Float64}

I also feel I probably should have added this extra point in yesterday however I did not give it much thought when solving the model using

ml = CellModel("/home/harry/Desktop/PhD/Year 1/Sem 2/CellML project/shi_hose_2009/shi_hose_2009-a679cdc2e974/ModelMain.cellml")
prob = ODEProblem(ml, (0.0,20.0))
sol = solve(prob, TRBDF2(), reltol=1e-6, abstol=1e-12)

I recieve the following warning

Warning: Unrecognized keyword arguments found. Future versions will error.
│ The only allowed keyword arguments to `solve` are: 
│ (:dense, :saveat, :save_idxs, :tstops, :d_discontinuities, :save_everystep, :save_on, :save_start, :save_end, :initialize_save, :adaptive, :abstol, :reltol, :dt, :dtmax, :dtmin, :force_dtmin, :internalnorm, :controller, :gamma, :beta1, :beta2, :qmax, :qmin, :qsteady_min, :qsteady_max, :qoldinit, :failfactor, :calck, :alias_u0, :maxiters, :callback, :isoutofdomain, :unstable_check, :verbose, :merge_callbacks, :progress, :progress_steps, :progress_name, :progress_message, :timeseries_errors, :dense_errors, :calculate_errors, :initializealg, :alg, :save_noise, :delta, :seed, :alg_hints, :kwargshandle, :trajectories, :batch_size, :sensealg, :advance_to_tstop, :stop_at_next_tstop, :default_set, :second_time)
│ See https://diffeq.sciml.ai/stable/basics/common_solver_opts/ for more details.
│ 
│ Set kwargshandle=KeywordArgError for an error message and more details.
│ Set kwargshandle=KeywordArgSilent to ignore this message.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/dVZCn/src/solve.jl:346

I checked this morning and this warning is not present when I formulate an ODE problem not using CellML, could this be related?

Cheers

5

That sounds like a problem with the CellML reader then. Open an issue there with the file you’re trying to read. Do you also have the SBML from the SBMLBiomodels repository? That could be a possible way to read it in: the CellML reader has had less work than SBMLToolbox.jl

There was a version of MTK which was passing an argument that it shouldn’t’ve. That got cleaned up and the latest versions should be fine if you upgrade.