[Solved] Time evaluation issue in a non-autonomous ODE

Hello everyone, I found a weird situation that a time-dependent parameter seemed to be ignored in my calcium dynamics model. After some trimming, I had an MWE here. Any suggestions/corrections are welcome.

Edit: solved by @MatFi I overwrote t in the following code. Sometimes it’s hard to see my own errors. :sweat_smile:

using Parameters, DifferentialEquations

#######################
# Adenylate kinase
# Equilibrium constant of 0.931 from Golding, 1995
#######################

const KEQ_AK = 0.931
const _p = 4KEQ_AK - 1

"Convert energy charge to adenylate fractions"
function adenylates(x)
    # 0 ≤ x ≤ 1 || throw(DomainError(x, "adenylates(x) is only defined for 0 ≤ x ≤ 1."))
    d = (sqrt(1 + 4_p * x * (1 - x))-1)/_p
    t = x - 0.5 * d
    m = 1 - t - d
    return (t = t, d = d, m = m)
end

"Cytosolic calcium oscillator, sine wave"
@with_kw struct CaOsci
    ca_r = 0.09e-3             # Resting cytosolic Ca concentration
    amplitude = 0.25e-3        # Calcium amplitude
    f2   = 2 * inv(60)         # Calcium oscillation frequency
end

"Cytosolic calcium, sine oscillation "
get_ca_c(p::CaOsci, t) = p.ca_r + 0.5 * p.amplitude * (1 + sinpi(p.f2 * t))

# The model
function model(u, p, t)
    ca_m = u

    # ec_c was also a state variable, but I put it here to make a MWE
    # If you remove the following 4 lines , time evaluation works fine
    ec_c = 0.9
    t, d, m = adenylates(ec_c)
    atp_c = t
    adp_c = d

    ca_c = get_ca_c(p, t)

    # Just curious, I added this to see t and ca_c evaluations
    # The result was the same if you comment out the line below
    @show t, ca_c

    # Simplified expressions for calcium dynamics
    jUni = ca_c - ca_m
    jNCLX = ca_m^2

    dca_m = jUni - jNCLX
    return dca_m
end

prob = ODEProblem(model, 0.250e-3, 1000.0, CaOsci())
sol = solve(prob)

The ca_m result was a smooth curve, rather than a sinusoid one like that of ca_c. I added the line @show t, ca_c and saw (t, ca_c) = (0.8252290521896939, 0.00022578878282531064) was kept evaluated but none for t >1.

What was even more interesting was that if you remove the lines indicated below, the evaluations of t and ca_m went back to normal. In this case t was evaluated all the way from 0.0 to 1000.0.

# The ATP part of my model, and trimmed down to the minimum.
ec_c = 0.9
t, d, m = adenylates(ec_c)
atp_c = t
adp_c = d

I’m not sure what caused this difference. Any suggestions and/or corrections are welcome.

System information

OS: Ubuntu 20.04
CPU: Intel(R) Xeon(R) CPU E5-2620 v4 @ 2.10GHz
Julia version: 1.53

Julia packages

[0c46a032] DifferentialEquations v6.15.0
[d96e819e] Parameters v0.12.1

#  Less relevant packages
[91a5bcdd] Plots v1.9.1
[7073ff75] IJulia v1.23.1
[b964fa9f] LaTeXStrings v1.2.0
[2ee39098] LabelledArrays v1.4.0
[d330b81b] PyPlot v2.9.0

from this point on t is constant, so is every parameter depending on t below this line.
Are you sure that you want to overwrite the time variable t in model(u, p, t) with a constant value?

1 Like

Ah, that was an oversight of my own. Thank you @MatFi for taking the time to see my code.

MatFi via JuliaLang <julialang@discoursemail.com> 於 2020年12月8日 週二 22:01 寫道: