Chris, here are some more details.

```
function myCallback(u, t, integrator)
C = @view(u[1:1:2]) # Calcium
Ce = @view(u[3:1:4]) # Calcium in ER
I = @view(u[5:1:6]) # IP3
h = @view(u[7:1:8])
J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, p)
end
```

I am returning several variables. There are 8 equations, two equations for C, two for Ce, etc.

I setup the callbacks:

```
saved_values = SavedValues(Float64, Vector::Any)
cb = SavingCallback(callback, saved_values)
prob = ODEProblem(Cresswell!, u0, tspan, pars) # d is a parameter dictionary
sol = solve(prob, Tsit5(), callback=cb)
```

I thought that Vector::Any was completely general in the sense that Vector elements could be any type. But that is not the case according to the first line of the error message. Here is a zoomed view:

Each of the return elements is an array of two elements (could be more in general). I want the capability of adding more variables without having to change the return type. So I am concerned with usability, not efficiency, since this is a debugging tool for me at this stage.

I did check out https://docs.juliadiffeq.org/release-5.0/features/callback_library.html

(I cannot copy/paste this error. The part in read is not copiable. Why would that be?)

I would like to also know how bring in my parameter array “p” into the myCallback() function. I could access it as a global variable, but that is poor programming.

Thanks!

Here is the full code.

```
using DifferentialEquations, Plots, BenchmarkTools
function globalParametersEvan()
# Try to duplicate Evan's model as much as possible
d = Dict()
d[:Rgas] = 8.31 # joule / kelvin / mole
d[:D_I] = 5.3e-10 # Hz * meter^2
d[:D_C] = 5.3e-10 # Hz * meter^2
d[:D_CE] = 5.3e-10 # Hz * meter^2
d[:Cm] = 1.e-6 # Farad * meter^-2
d[:e] = 1.602176634e-19 # Amp * second
d[:N_A] = 6.02214076e23 # mole^-1
d[:k_B] = 1.38064852e-23 # meter^2 * Kg * Hz^2 / kelvin
d[:F] = 96485.3329 # Amp * second / mole = Cb / mole
d[:V_T] = 0.026 # voltage_cubic
d[:P_Ca] = 4.46e-15 # meter * Hz
d[:Crest] = 0.073e-6 # Molar
d[:Vrest] = -80. # Volt
d[:Ce] = 1.800 # mole/meter**3]
d[:Ce] = 1.e-4 # mole/meter**3] # WHICH ONE TO CHOOSE?
d[:O_delta] = 0.6e-3 # Hz * mole/meter^3
d[:k_delta] = 1.5e-3 # mole/meter^3
d[:K_delta] = 0.1e-3 # mole/meter^3
d[:K_D] = 0.7e-3 # mole/meter^3
d[:K_3] = 1.0e-3 # mole/meter^3
d[:O_3K] = 4.5e-3 # Hz*mole/meter^3
d[:o_5P] = 0.05 # Hz
d[:K_5P] = 10.0e-3 # mole/meter^3
d[:o_delta] = 0.6e-3 # mole/meter^3
d[:k_delta] = 1.5e-3 # Hz * mole/meter^3
d[:K_delta] = 0.1e-3 # mole/meter^3
d[:Lambda] = 2100e-18 # meter^3
d[:rho_A] = 0.18 # ER-to-cytoplasm volume ratio ?
d[:rho] = d[:rho_A] / (1. + d[:rho_A]) # # ER-to-cytoplasm volume ratio ?
# Multiply value of amplitudes by cytosolic volume (changed rho_A to rho)
# It does not look as if I did the multiplication
#self.o_delta = 0.6 # * umolar * Lambda * (1-rho) / second
d[:o_delta] = 0.6e-3 # * mole/meter**3 * Lambda * (1-rho) / second
d[:o_3K] = 4.5e-3 # * mole/meter**3 * Lambda * (1-rho) / second
d[:Omega_L] = 0.1 #/second # Maximal rate of Ca^2+ leak from the ER
d[:Omega_2] = 0.2 #/second # IP_3R binding rate for Ca^2+ inhibition
d[:Omega_u] = 0.2-3 #/second # uptake/dissassociation constant
# --- IP_3R kinectics
d[:d_1] = 0.13e-3 #*mole/meter**3 # IP_3 binding affinity
d[:d_2] = 1.05e-3 #*mole/meter**3 # Ca^2+ inactivation dissociation constant
d[:O_2] = 0.2e3 #/(mole/meter**3)/second # IP_3R binding rate for Ca^2+ inhibition
d[:d_3] = 0.9434e-3 #*mole/meter**3 # IP_3 dissociation constant
d[:d_5] = 0.08 #*mole/meter**3 # Ca^2+ activation dissociation constant
d[:p_open] = 1
d[:P_r] = 1e-3 #*mole/meter**3/second
d[:P_CE] = 1e-3 #*mole/meter**3/second
d[:eta_p] = 1
d[:d_ER] = 1e-6 #* meter
d[:dv_ER] = 10e-3
d[:K_P] = 0.05e-3
d[:L] = 8e-6 # umeter
d[:R] = 3e-6 # meter
d[:r] = 2e-6 # meter
d[:O_delta] = 0.6e-21 # mole/meter**3 * Lambda * (1-rho) / second
d[:O_5P] = 0.05e-21 # mole/meter**3 * Lambda * (1-rho) / second
d[:F] = 0.1 # ALREADY DEFINED?
d[:D] = 0.05
# Evan parameters
d[:vIPR] = 6 # sec^-1
d[:vSERCA] = 4.4e-3 # mole m^-3 s^-1 (1 Molar = 1 Mole/liter = 10^3 mole/meter^3)
d[:vLeak] = 0.11 # sec^-1
d[:d1] = 0.13e-3 # mole m^-3
d[:d5] = 0.08234e-3 # mole m^-3
# open channels
d[:d2] = 1.049e-3 # mole m^-3
d[:d1] = 0.13e-3 # mole m^-3
d[:d3] = 0.9434e-3 # mole m^-3
d[:a2] = 0.20e-3 # mole m^-3 s^-1
# IP3 (J_δ)
d[:o_δ] = .025e-3 # Mole m^-3 s^-1
d[:K_δ] = 0.5e-3 # Mole m^-3, Ca2+ affinity of PLCδ
d[:k_δ] = 1.e-3 # Mole m^-3, IP3 inhibiting affinity // Oschmann2009: 1.5 | COMPGLIO pg130: 1.
# End Evan Parameters
return d
end
lJp = []
lJr = []
lJ1 = []
lJleak = []
lJSERCA = []
lJIPR = []
lJβ = []
lCount = []
lTime = []
function getCurrents(C, Ce, I, h, t, p)
minf = @. Hill(I,p[:d1],1) * Hill(C,p[:d5],1)
Jleak = @. p[:vLeak] * (Ce-C)
JIPR = @. p[:vIPR] * (minf*h)^3 * (C-Ce)
JSERCA = @. p[:vSERCA] * Hill(C,p[:d5],2)
J_β = sin(2*π*t)^10 #*mole/meter**3/second : mole/meter**3/second
J_δ = @. p[:o_δ] * Hill(p[:k_δ], I, 1) * Hill(C, p[:K_δ], 2)
J5P = @. p[:o_5P] * Hill(I, p[:K_5P], 1) #: mole/second # o_5P, K_5P
J3K = @. p[:o_3K] * Hill(C, p[:K_D], 4) * Hill(I, p[:K_3], 1) #: mole/second # o_3K, K_D, K_3
Q2 = @. p[:d2] * (I + p[:d1]) / (I + p[:d3])
OmegaH = @. p[:a2] * (Q2 + C)
hinf = @. p[:d_2] * (I + p[:d_1]) / (p[:d_2]*(I + p[:d_1]) + (I + p[:d_3])*C) #: 1 # d_2, d_1, d_3
return J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2
end
function Cresswell!(du, u, p, t)
C = @view(u[1:1:2]) # Calcium
Ce = @view(u[3:1:4]) # Calcium in ER
I = @view(u[5:1:6]) # IP3
h = @view(u[7:1:8]) # fraction of open channels
push!(lCount, 1)
# rhs_C, etc. Are arrays, created via implicit vectorization involving numpy arrays
dR2 = p[:R]^2 - p[:r]^2 #: meter**2
s = p[:R] + p[:r]
V = @. p[:Vrest] + (p[:F]*s/p[:Cm]) * (C - p[:Crest]) #: volt
VVT = -2. * V/p[:V_T] #: 1 # value of 0, which leads to a zero denominator
J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, p)
p[:J_δ] = J_δ
push!(lJleak, Jleak)
push!(lJIPR, JIPR)
push!(lJSERCA, JSERCA)
push!(lJβ, J_β)
push!(lTime, t)
dh = @. OmegaH * (hinf - h) #: 1
dI = @. (J_β + J_δ - J3K - J5P)
# GE: divde Jp by Lambda (2020-03-09)
dC = @. (JIPR - JSERCA + Jleak) #: mole / meter**3
dCE = @. -(1. / p[:rho_A]) * (JIPR - JSERCA + Jleak) #: mole / meter**3
dI = @. (1. *J_β + 1. * J_δ - 1. *J3K - 1. *J5P) # / (p[:Lambda]*(1-p[:rho])) #: mole/meter**3
dh = @. OmegaH * (hinf - h) #: 1
du[1:2] = dC
du[3:4] = dCE #[0,0]
du[5:6] = dI # [0,0] # blows up!
du[7:8] = dh
end
nb_neurons = 2
eqs_per_neuron = 4
# Initial Conditions
function initialConditions()
c0 = rand(nb_neurons)
ce0 = rand(nb_neurons)
h0 = rand(nb_neurons)
I0 = rand(nb_neurons)
# Initial condition vector
u0 = vcat(c0,ce0,h0,I0)
end
function callback(u, t, integrator)
C = @view(u[1:1:2]) # Calcium
Ce = @view(u[3:1:4]) # Calcium in ER
I = @view(u[5:1:6]) # IP3
h = @view(u[7:1:8])
# How to I get my parameter array as a local variable?
J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, pars)
end
getParams = globalParameters_mod.globalParametersEvan
saved_values = SavedValues(Float64, Vector::Any)
cb = SavingCallback(callback, saved_values)
u0 = initialConditions()
pars = getParams()
tspan = (0., 300.) # make these real
#u0 = [.1, .2, .3, .4]``
# We have 8 equations. How to collect them?
prob = ODEProblem(Cresswell!, u0, tspan, pars) # d is a parameter dictionary
sol = solve(prob, Tsit5(), callback=cb)
```

Here is a modified version of the callback unction:

```
function callback(u, t, integrator)
C = @view(u[1:1:2]) # Calcium
Ce = @view(u[3:1:4]) # Calcium in ER
I = @view(u[5:1:6]) # IP3
h = @view(u[7:1:8])
# How to I get my parameter array as a local variable?
J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, pars)
cur = Dict()
cur[:J5P] = J5P
cur[:J3K] = J3K
cur[:J_β] = J_β
cur[:J_δ] = J_δ
cur[:Jleak] = Jleak
cur[:JIPR] = JIPR
cur[:JSERCA] = JSERCA
cur[:hinf] = hinf
cur[:OmegaH] = OmegaH
cur[:minf] = minf
cur[:Q2] = Q2
return cur
end
saved_values = SavedValues(Float64, Dict{Any,Any})
cb = SavingCallback(callback, saved_values)
```

and this “seems” to work. I do not yet know how to use this callback once the solver has completed.

But my immediate comment is that my callback is now substantially longer than it was and more unwieldy. I could have populated a dictionary in the getCurrents() method, but then I would have to access these dictionary values in the ODE solver. I assume that accessing dictionary values is very cheap?

I have another question. If I had an array:

```
vars = [:J5P, :J3K]
```

for example, would it be possible to create a macro to execute:

cur[:J5P] = J5P

cur[:J3K] = J3K

?

Thanks,