The following is a reduction of a 40 dimensional system (that’s why some of the syntax is weird). This simulation is crashing and I am trying to understand why. sim.t[end] shows the last time step at 672.6064966897368. However, when I run @show, I see the integrator is taking a couple more steps. I presume these are rejected steps – but I can’t be sure. It looks like the integrator thinks the step size can grow. But when the step size becomes too big I compute the log of a negative number. I’m trying to understand how the error is being computed and why the integrator is increasing its stepsize. I remember seeing the formula for the error floating around in the documentation.
- Does this sound like reasonable strategy to troubleshoot?
- Is there a way to use
@showto report the components of the error formula in the documentation. Maybe this will give me insight on what is failing. - My understanding is that
Tsit5()requires computing up to 5 internal stages, however, when I runsim.k[end], I can see integrator computing 7 stages. Where on github is the code for Tsit5()? I’d like to see what is actually being computed. - I have a hunch that there is something going on with the function
tauNaHbecause when I replace this with a number, the simulations run unimpeded. But I’m having trouble understanding more.
thanks!
######################################################################################
using DifferentialEquations, NaNMath
endd = 20000000.0
#####################################################################
### Parameters
## Reversal Potentials
EK(t) = -80;
EH = -20.0;
ENa = 30.0;
EL = -50.0;
#####################################################################
### Define generic functions
## Generic Activation/Inactication Curves
sigmd(x,translate,center,reactance) = 1/(1+exp(((x-translate) + center)/reactance));
## Define the time constant function for ion channels
tauX(Volt, CT, DT, AT, BT, halfAc) = CT - DT/(1 + exp(((Volt-halfAc) + AT)/BT));
## Some intrinsic currents require a special time constant function
spectau(Volt, CT, DT, AT, BT, AT2, BT2, halfAc) = CT + DT/(exp(((Volt-halfAc) + AT)/BT) + exp(((Volt-halfAc) + AT2)/BT2));
## Define the ionic currents; q is the exponent of the activation variable m
iIonic(g, m, h, q, Volt, Erev) = g*(m^q)*h*(Volt - Erev);
#####################################################################
# Ion Channel Acitvation/Inactivation Curves
NaMinf(V,trns) = sigmd(V,trns, 25.5, -5.29); # m^3
NaHinf(V,trns) = sigmd(V,trns, 48.9, 5.18); # h
CaTMinf(V,trns) = sigmd(V,trns, 27.1, -7.20); # m^3
CaTHinf(V,trns) = sigmd(V,trns, 32.1, 5.50); # h
CaSMinf(V,trns) = sigmd(V,trns, 33.0, -8.1); # m^3
CaSHinf(V,trns) = sigmd(V,trns, 60.0, 6.20); # h
HMinf(V,trns) = sigmd(V,trns, 70.0, 6.0); # m
KdMinf(V,trns) = sigmd(V,trns, 12.3, -11.8); # m^4
KCaMinf(V,trns,IntCa) = (IntCa/(IntCa + 3.0))*sigmd(V,trns, 28.3, -12.6); # m^4
AMinf(V,trns) = sigmd(V,trns, 27.2, -8.70); # m^3
AHinf(V,trns) = sigmd(V,trns, 56.9, 4.90); # h
# Time Constants (ms)
tauNaH(V,trns) = tauX(V, 0.00, -0.67, 62.9, -10.0, trns)*tauX(V, 1.50, -1.00, 34.9, 3.60, trns);
# Calcium Reverse Potential
R = 8.314*1000.0; # Ideal Gas Constant (*10^3 to put into mV)
temp = 10.0; # Temperature in Celcius; Temperature of the Sea lawl
z = 2.0; # Valence of Calcium Ions
Far = 96485.33; # Faraday's Constant
CaOut = 3000.0; # Outer Ca Concentration (uM)
ECa(CaIn) = ((R*(273.15 + temp))/(z*Far))*NaNMath.log(CaOut/CaIn);
#ECa(CaIn) = ((R*(273.15 + temp))/(z*Far))*log(CaOut/CaIn);
# Ionic Currents (mV / ms)
iNa(g,m,h,V) = iIonic(g, m, h, 3.0, V, ENa);
iCaT(g,m,h,V,CaIn) = iIonic(g, m, h, 3.0, V, ECa(CaIn));
iCaS(g,m,h,V,CaIn) = iIonic(g, m, h, 3.0, V, ECa(CaIn));
iH(g,m,V) = iIonic(g, m, 1.0, 1.0, V, EH);
iKd(g,m,V,t) = iIonic(g, m, 1.0, 4.0, V, EK(t));
iKCa(g,m,V,t) = iIonic(g, m, 1.0, 4.0, V, EK(t));
iA(g,m,h,V,t) = iIonic(g, m, h, 3.0, V, EK(t));
iL(V) = iIonic(.01, 1.0, 1.0, 1.0, V, EL);
############################################################################################
function f(du,u,p,t)
du[1] = -iL(u[1])-iNa(u[14],u[2],u[3],u[1])-iCaT(u[15],u[4],u[5],u[1],u[13])-iCaS(u[16],u[6],u[7],u[1],u[13])-iH(u[17],u[8],u[1])-iKd(u[18],u[9],u[1],t)-iKCa(u[19],u[10],u[1],t)-iA(u[20],u[11],u[12],u[1],t)
du[2] = (NaMinf(u[1],u[21]) - u[2])/.1
du[3] = (NaHinf(u[1],u[22]) - u[3])/tauNaH(u[1],u[22])
du[4] = (CaTMinf(u[1],u[23]) - u[4])/6
du[5] = (CaTHinf(u[1],u[24]) - u[5])/40
du[6] = (CaSMinf(u[1],u[25]) - u[6])/18
du[7] = (CaSHinf(u[1],u[26]) - u[7])/100
du[8] = (HMinf(u[1],u[27]) - u[8])/700
du[9] = (KdMinf(u[1],u[28]) - u[9])/5
du[10] = (KCaMinf(u[1],u[29],u[13]) - u[10])/50
du[11] = (AMinf(u[1],u[30]) - u[11])/10
du[12] = (AHinf(u[1],u[31]) - u[12])/30
du[13] = (-.94*(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13]))-u[13]+.05)/20
@show t
end
############################################################################################
#Start
ICs = [-29.337327051014498 0.387305560452983 0.520851480339422 0.15123291866400054 0.7788086343207855 0.18292525224066983 0.40518167654812637 0.40305741894377733 0.09880328448541692 0.1192742729813851 0.17090849597998528 0.4853554113680752 2.455375230646913 79.93706465236964 0.2873998785854699 0.0318308211031974 0.15906997203306664 99.67308490480416 0.8899609102246536 2.6818203693916858 -6.016940524541976 6.0497029078690145 -0.08934230228658809 0.1845828676421055 0.11586968152835098 0.7364356646257856 2.4826107420682773 -6.238217364811495 -1.966199265610444 -2.1624666424135977 2.014861831075435 4.143056495983624e-7 0.999901341448654 0.16812145388254565 0.5889709108552483 0.4212141706080315 0.24284481424824048 -0.03505188961118154 -0.16940390920534412 1.0];
tspan = (0.0,endd);
SaveTimes = (0.0:.1:endd); #in milliseconds
prob = ODEProblem(f,ICs,tspan);
sim = solve(prob,Tsit5(),maxiters=1e10);