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. 
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