Why the solution of a system of ODEs is different when you run the code a second time?

New to Julia
, ,
1

Hello, I am trying to study the effect of a disturbance on a system of ODEs and when I run the code more than once it gives me a slightly different result compared to when I ran it for the first time. I can’t find the cause of this. When I use callbacks to apply the perturbation instead of solving the systems twice the result is the same as the first run. Can someone please help me figure out what’s going on?

Here is my code, hope it’s not too confusing:

Code 
using DifferentialEquations
using ModelingToolkit
using Plots

# Define ODE systems
module model
    using ModelingToolkit
    @parameters mY n gZ1 theta gZ2 bI g1 bC gC
    @variables t Z1(t) Z2(t) X1(t) XC(t)
    D = Differential(t)

    @named odeFB = ODESystem([D(Z1) ~ mY - (n * Z1 * Z2) - (gZ1 * Z1),
                                D(Z2) ~ (theta * XC) - (n * Z1 * Z2) - (gZ2 * Z2),
                                D(X1) ~    (bI * Z1) - (g1 * X1),
                                D(XC) ~    (bC * X1) - (gC * XC)])

    # External function to obtain xc(t) 
    xc(t) = solFBR(t)[4]
    @register_symbolic xc(t)
        
    @named odeNF = ODESystem([D(Z1) ~ mY - (n * Z1 * Z2) - (gZ1 * Z1),
                                D(Z2) ~ (theta * xc(t)) - (n * Z1 * Z2) - (gZ2 * Z2),
                                D(X1) ~    (bI * Z1) - (g1 * X1),
                                D(XC) ~    (bC * X1) - (gC * XC)])
end

# Set up the parameter values
@parameters mY n gZ1 theta gZ2 bI g1 bC gC
p = Dict([
    bI => 0.06, # min-1
    mY => 1 * 60, # min-1 nM 
    n => 0.01, # nM-1 min-1
    gZ1 => 0, # nM min-1 [0.0301,0.0408]
    gZ2 => 0, # nM min-1 [0.0301,0.0408]
    theta => 0.3, # min-1
    g1 => 0.1, # min-1 
    bC => 1.06*1.05, # min-1 [0.01,0.05565,0.1,0.106]
    gC => 0.1 # nM
]);
pV = collect(pairs(p))

# Set the initial state 
x0 = solve(ODEProblem(model.odeFB,ones(4),1e6,pV); reltol = 1e-12, save_everystep = false).u[2];

# Solve odeFB
solFB1 = solve(ODEProblem(model.odeFB, x0, 100, pV), alg_hints=[:stiff]; reltol = 1e-12);
# Pass the solution to the module that contains the differential equations
@eval model solFBR = Main.solFB1;
# Solve odeNF
solNF1 = solve(ODEProblem(model.odeNF, x0, 100, pV), alg_hints=[:stiff]; reltol = 1e-12);
# Solve odeFB with different initial conditions and tspan 
solFB2 = solve(ODEProblem(model.odeFB, last(solFB1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
@eval model solFBR = Main.solFB2;
# Apply a perturbation to parameter mY
p[mY] *= 1.05
pV = collect(pairs(p))
# Again solve both systems of differential equations
solFB3 = solve(ODEProblem(model.odeFB, last(solFB1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
solNF2 = solve(ODEProblem(model.odeNF, last(solNF1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
p[mY] /= 1.05

# Plot
solFBt = vcat(solFB1.t, solFB3.t .+ 100);
solFBu = vcat(solFB1[4,:], solFB3[4,:]);

solNFt = vcat(solNF1.t, solNF2.t .+ 100);
solNFu = vcat(solNF1[4,:], solNF2[4,:]);

plot(solNFt, solNFu, label = "NF")
plot!(solFBt, solFBu, label = "FB")

Here is the output of the first run:

Screenshot from 2023-06-01 01-29-19
Screenshot from 2023-06-01 01-29-19721×482 20.6 KB

Here is the output of the second run:
Screenshot from 2023-06-01 01-29-37
Screenshot from 2023-06-01 01-29-37716×479 21.1 KB

2

I tried removing the model module, and then I get consistent results. Not sure why it would have that effect, probably something funky going on with replacing the module and using eval that is over my head…

Code 
using DifferentialEquations
using ModelingToolkit
using Plots
using ModelingToolkit

@parameters mY n gZ1 theta gZ2 bI g1 bC gC
@variables t Z1(t) Z2(t) X1(t) XC(t)
D = Differential(t)

@named odeFB = ODESystem([D(Z1) ~ mY - (n * Z1 * Z2) - (gZ1 * Z1),
                            D(Z2) ~ (theta * XC) - (n * Z1 * Z2) - (gZ2 * Z2),
                            D(X1) ~    (bI * Z1) - (g1 * X1),
                            D(XC) ~    (bC * X1) - (gC * XC)])

# External function to obtain xc(t) 
xc(t) = solFBR(t)[4]
@register_symbolic xc(t)
    
@named odeNF = ODESystem([D(Z1) ~ mY - (n * Z1 * Z2) - (gZ1 * Z1),
                            D(Z2) ~ (theta * xc(t)) - (n * Z1 * Z2) - (gZ2 * Z2),
                            D(X1) ~    (bI * Z1) - (g1 * X1),
                            D(XC) ~    (bC * X1) - (gC * XC)])

# Set up the parameter values
@parameters mY n gZ1 theta gZ2 bI g1 bC gC
p = Dict([
    bI => 0.06, # min-1
    mY => 1 * 60, # min-1 nM 
    n => 0.01, # nM-1 min-1
    gZ1 => 0, # nM min-1 [0.0301,0.0408]
    gZ2 => 0, # nM min-1 [0.0301,0.0408]
    theta => 0.3, # min-1
    g1 => 0.1, # min-1 
    bC => 1.06*1.05, # min-1 [0.01,0.05565,0.1,0.106]
    gC => 0.1 # nM
]);
pV = collect(pairs(p))

# Set the initial state 
x0 = solve(ODEProblem(odeFB,ones(4),1e6,pV); reltol = 1e-12, save_everystep = false).u[2];

# Solve odeFB
solFB1 = solve(ODEProblem(odeFB, x0, 100, pV), alg_hints=[:stiff]; reltol = 1e-12);
# Pass the solution to the module that contains the differential equations
solFBR = solFB1;
# Solve odeNF
solNF1 = solve(ODEProblem(odeNF, x0, 100, pV), alg_hints=[:stiff]; reltol = 1e-12);
# Solve odeFB with different initial conditions and tspan 
solFB2 = solve(ODEProblem(odeFB, last(solFB1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
solFBR = solFB2;
# Apply a perturbation to parameter mY
p[mY] *= 1.05
pV = collect(pairs(p))
# Again solve both systems of differential equations
solFB3 = solve(ODEProblem(odeFB, last(solFB1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
solNF2 = solve(ODEProblem(odeNF, last(solNF1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
p[mY] /= 1.05

# Plot
solFBt = vcat(solFB1.t, solFB3.t .+ 100);
solFBu = vcat(solFB1[4,:], solFB3[4,:]);

solNFt = vcat(solNF1.t, solNF2.t .+ 100);
solNFu = vcat(solNF1[4,:], solNF2[4,:]);

plot(solNFt, solNFu, label = "NF")
plot!(solFBt, solFBu, label = "FB")
1 Like
3

Yeah this is really funky. Generally if you’re doing something with @eval you probably want to consider alternative styles to avoid it.

If p is the same vector in both problems, then this will mutate both.

2 Likes
4

I think I’ve figured out what’s happening. The first simulation gets interrupted, and since I use the last values of the solution as the initial state for the second simulation, it gives me a different result. In fact, the code did give me a warning, but I ignored it :stuck_out_tongue: I suppose I can fix the problem by specifying the algorithm.

Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. ODE Solvers · DifferentialEquations.jl).

1 Like