Error building model

[lhdp0110]

Hello, I am using the InformationGeometry.jl package for a project and am encountering an issue at the data model building step. I followed along with the "ODE-based models" advanced example. My project consists of using MLE to model the output from an immune system pathway. The output consists of 5 time points, and the pathway is described by a system of 9 ODEs. I am including my code below. The error given is:

┌ Warning: Interrupted. Larger maxiters is needed. └ @ SciMLBase C:\Users\lou-i.julia\packages\SciMLBase\OHiiA\src\integrator_interface.jl:331 ERROR: "ODE integration failed, maybe try using a lower tolerance value. θ=[0.1, 0.02, 0.02, 0.02, 1.0, 1.0, 0.05, 0.05, 0.05, 0.1]."

I have tried different tol values but keep getting the same error. The AMPs values I've put in are theoretical and should fit the model quite well. Could you please advise? Thank you in advance!

using InformationGeometry, ModelingToolkit, Plots
@parameters t DP cP cI cB cN cNb phiP phiG phiI phiR
@variables P(t) Pb(t) Ib(t) RB(t) RP(t) RN(t) RNb(t) H(t) A(t)
Dt = Differential(t)

V = 1
VN = 0.1
G = 5
m = 2
k = 2
deltaP = 0.01
deltaPb = 0.02
deltaR = 0.01
deltaN = 0.01
cH = 0.05
cA = 0.1
deltaH = 0.05
deltaA = 0.05
rhoIb = 0.01
KR = 2

QP = deltaP
IT = 1
QR = deltaR/V
QH = deltaH/V
QA = deltaA/V

IMDeqs = [ Dt(P) ~ QP - m * cP * exp(-phiP*H) * P^k * G * exp(-phiG*A) - deltaP * P,
        Dt(Pb) ~  cP * exp(-phiP*H) * P^k * G * exp(-phiG*A) - deltaPb + Pb,
        Dt(Ib) ~ cI * exp(-phiI*H) * (IT - Ib) * Pb - rhoIb * Ib,
        Dt(RB) ~  QR/V + (phiR*RNb)/(KR+RNb) - cB * Ib * RB - deltaR*RB,
        Dt(RP) ~ cB * Ib * RB - DP * RP - deltaR * RP,
        Dt(RN) ~ DP * RP / VN - cN * RN + cN * RNb - deltaN * RN,
        Dt(RNb) ~ cN * RN - cNb * RNb,
        Dt(H) ~ QH / V + cH * RNb - deltaH * H,
        Dt(A) ~ QA / V + cA * RNb - deltaA * A]

IMDstates = [P, Pb, Ib, RB, RP, RN, RNb, H, A]
IMDparams = [DP, cP, cI, cB, cN, cNb, phiP, phiG, phiI, phiR]
@named IMDsys = ODESystem(IMDeqs, t, IMDstates, IMDparams)

P0 = 1.0
Pb0 = 0.0
Ib0 = 0.0
RB0 = 1.0
RP0 = 0.0
RN0 = 0.0
RNb0 = 0.0
H0 = 1.0
A0 = 1.0
IMDinitial = [P0, Pb0, Ib0, RB0, RP0, RN0, RNb0, H0, A0]

hours = [0, 120, 240, 480, 720]
AMPs = [0, 2, 7.5, 15, 14]
confInts = [0.5, 0.5, 0.5, 0.5, 0.5]
IMDdataSet = DataSet(hours, AMPs, confInts)

IMDobservables = [7]
IMDdataModel = DataModel(IMDdataSet, IMDsys, IMDinitial, IMDobservables, [0.10, 0.02, 0.02, 0.02, 1, 1, 0.05, 0.05, 0.05, 0.10]; tol=0.1)

[RafaelArutjunjan]

Hello there,

by default, InformationGeometry.jl sets the relative tolerance of ODE models to 1e-7, so tol=0.1 is far too large for the ODE integrator to produce useful solutions.

Also, due to the exponential dependence of your ODE system on the states, your model looks to be relatively "stiff", whereas the default ODE solver choice of InformationGeometry.jl is Tsit5(), which is a "non-stiff" solver (see e.g. DifferentialEquations.jl). In this case, a better choice is given by AutoTsit5(Rosenbrock23()) which automatically detects stiffness and switches between solvers accordingly.

I will try to make this more explicit in the documentation and will probably change the default from Tsit5() to AutoTsit5(Rosenbrock23()) in the future.

However, there is yet another problem: even using only the OrdinaryDiffEq.jl framework to compute numerical solutions to the ODE system you provided without any kind of parameter inference already fails for me.

using OrdinaryDiffEq, Plots
tspan = (0.,710.)
Prob = ODEProblem(ODEFunction(IMDsys), IMDinitial, tspan, [0.10, 0.02, 0.02, 0.02, 1, 1, 0.05, 0.05, 0.05, 0.10])
sol = solve(Prob, AutoTsit5(Rosenbrock23()); reltol=1e-7)
plot(IMDdataSet);       plot!(sol, vars=(0,7))

Using tspan = (0., 710.) produces solutions just fine, but when setting it to tspan = (0., 720.), which is required to include your last time point, the integration appears to fail due to some kind of instability.

Although this too might just be a solver choice issue, I was not really able to get it to work just by playing around with the tolerances and solvers. Are you sure the specification of the ODE system (including its parameters and initial values) you provided is correct? Can you maybe provide a reference of some kind?

[lhdp0110]

Thank you for your prompt and thorough response,

I understand that the tol value in this script is too high, this is one of several that I played around with, ranging from 1e-15 to 0.1 and I encountered the same issue in any case. I am attaching the reference from which I took this model. Refer specifically to page 14 for the ODEs, and page 15 for the parameters (in the table and its caption) and initial values (in the paragraph written paragraph). I estimated the observed values from the graph on page 16; I will use experimental data once I've established that my program works. The paper uses a program written in R, but since I am a complete beginner in computational modeling and I am not familiar with R, I have had a hard time figuring out how it works. However, the ODE system does look the same as in my program. I do have the R file, but am unable to attach it here. I'd be happy to provide it via email or just copy/pasting it here if that would be helpful.

rspb20210786_si_001.PDF

[lhdp0110]

I forgot to mention that the parameters given in the table on page 15 failed to produce the same graph as on page 16 for me. I had to take the values in the R file in order to produce a similar curve, and those values are not found in the paper. Here is a script I wrote to check that I could reproduce the curve using the ODEs and the parameters given:

using DifferentialEquations
using Plots

DP = 0.10
cP = 0.02
cI = 0.02
cB = 0.02
cN = 1
cNb = 1
phiP = 0.05
phiG = 0.05
phiI = 0.05
phiR = 0.10

V = 1
VN = 0.1
G = 5
m = 2
k = 2
deltaP = 0.01
deltaPb = 0.02
deltaR = 0.01
deltaN = 0.01
cH = 0.05
cA = 0.1
deltaH = 0.05
deltaA = 0.05
rhoIb = 0.01
KR = 2

QP = deltaP
IT = 1
QR = deltaR/V
QH = deltaH/V
QA = deltaA/V

function pathwayEq!(du,u,p,t)
    du[1] = QP - m * cP * exp(-phiP * u[8]) * u[1]^k * G * exp(-phiG * u[9]) - deltaP * u[1]
    du[2] = cP * exp(-phiP * u[8]) * u[1]^k * G * exp(-phiG * u[9]) - deltaPb * u[2]
    du[3] = cI * exp(-phiI * u[8]) * (IT - u[3]) * u[2] - rhoIb * u[3]
    du[4] = QR / V + (phiR * u[7])/(KR+u[7]) - cB * u[3] * u[4] - deltaR * u[4]
    du[5] = cB * u[3] * u[4] - DP * u[5] - deltaR * u[5]
    du[6] = DP * u[5] / VN - cN * u[6] + cN * u[7] - deltaN * u[6]
    du[7] = cN * u[6] - cNb * u[7]
    du[8] = QH / V + cH * u[7] - deltaH * u[8]
    du[9] = QA / V + cA * u[7] - deltaA * u[9]
end

P0 = 1.0
Pb0 = 0.0
Ib0 = 0.0
RB0 = 1.0
RP0 = 0.0
RN0 = 0.0
RNb0 = 0.0
H0 = 1.0
A0 = 1.0
u0 = [P0;Pb0;Ib0;RB0;RP0;RN0;RNb0;H0;A0]
tspan = (0,800)

prob = ODEProblem(pathwayEq!,u0,tspan)
sol = solve(prob)
plot(sol, vars=(0,7), xaxis = "Time (hours)", yaxis = "Bound Nuclear Relish", label = false)

[RafaelArutjunjan]

Your last bit of code works for me as well.

In that case there is no problem and it should be relatively easy to make your example work with minimal changes, namely by modifying your definition of pathwayEq! such that it accepts parameters via the p argument instead of using the globally defined variables:

using InformationGeometry, OrdinaryDiffEq, Plots

function pathwayEq!(du,u,p,t)
    DP, cP, cI, cB, cN, cNb, phiP, phiG, phiI, phiR = p

    V = 1
    VN = 0.1
    G = 5
    m = 2
    k = 2
    deltaP = 0.01
    deltaPb = 0.02
    deltaR = 0.01
    deltaN = 0.01
    cH = 0.05
    cA = 0.1
    deltaH = 0.05
    deltaA = 0.05
    rhoIb = 0.01
    KR = 2
    QP = deltaP
    IT = 1
    QR = deltaR/V
    QH = deltaH/V
    QA = deltaA/V

    du[1] = QP - m * cP * exp(-phiP * u[8]) * u[1]^k * G * exp(-phiG * u[9]) - deltaP * u[1]
    du[2] = cP * exp(-phiP * u[8]) * u[1]^k * G * exp(-phiG * u[9]) - deltaPb * u[2]
    du[3] = cI * exp(-phiI * u[8]) * (IT - u[3]) * u[2] - rhoIb * u[3]
    du[4] = QR / V + (phiR * u[7])/(KR+u[7]) - cB * u[3] * u[4] - deltaR * u[4]
    du[5] = cB * u[3] * u[4] - DP * u[5] - deltaR * u[5]
    du[6] = DP * u[5] / VN - cN * u[6] + cN * u[7] - deltaN * u[6]
    du[7] = cN * u[6] - cNb * u[7]
    du[8] = QH / V + cH * u[7] - deltaH * u[8]
    du[9] = QA / V + cA * u[7] - deltaA * u[9]
end

P0 = 1.0
Pb0 = 0.0
Ib0 = 0.0
RB0 = 1.0
RP0 = 0.0
RN0 = 0.0
RNb0 = 0.0
H0 = 1.0
A0 = 1.0
u0 = [P0;Pb0;Ib0;RB0;RP0;RN0;RNb0;H0;A0]

hours = [0, 120, 240, 480, 720]
AMPs = [0, 2, 7.5, 15, 14]
confInts = [0.5, 0.5, 0.5, 0.5, 0.5]
IMDdataSet = DataSet(hours, AMPs, confInts; xnames=["Time (hours)"], ynames=["Bound Nuclear Relish"])

params = [0.10, 0.02, 0.02, 0.02, 1, 1, 0.05, 0.05, 0.05, 0.10]

IMDdataModel = DataModel(IMDdataSet, ODEFunction(pathwayEq!), u0, [7], params)

plot(IMDdataModel)

MLE(IMDdataModel)
# [0.0632, 0.0376, 0.0224, 0.02, 1.0, 0.998, 0.138, 0.098, 0.0177, 0.0604]

It's not even necessary to change the tolerances, so it appears like the whole "stiffness problem" came from some error in the definition, a typo somewhere maybe?

Some remarks:

• It's not necessary to use ModelingToolkit.jl to create ODESystem objects to input into the DataModel constructor, an ODEFunction also works. The only difference is that InformationGeometry.jl cannot infer the symbolic variable names of the concentrations and parameters from the ODEFunction for plotting and such.
• Although the model is calibrated to the given data, it is of course not possible to uniquely fit 10 parameters from 5 data points, so it would arguably be misleading to say that the resulting parameter values constitute a "maximum likelihood estimate".
• It may look weird to have the constants defined inside of pathwayEq! but the compiler is generally pretty good at optimizing superfluous variable definitions inside functions away, so there should be no performance downside. In fact, it should be faster than accessing globally defined variables.

I hope this solves your problem.

[lhdp0110]

Thank you very much for your extremely valuable insights. I was finally able to get my program to work with your help. Much gratitude!

[RafaelArutjunjan]

You're welcome, if you should encounter any more problems, feel free to open another issue.