Ah sorry, p
is just a vector of parameters. 'p[6] = 3, 'p[7] = 6
, and p[8] = 20
.
To answer the question why I’m collecting all my initial condition:
The documentation for how to use the ensemble problem is here. In particular one needs to define an EnsembleProblem
which is created using a constructor. prob_func
is part of this constructor – it takes the originally defined ODE and reformulates it. I was going to define all of my intial conditions beforehand, then index into the array of intial conditions using i
in prob_func
. Full code below:
using DifferentialEquations
#################################################
##### ODE params
# initialize parameter vector
p = fill(NaN,8);
# Time Constants for Homeostasis -- (g_Ca, g_K, nHalf)
p[1] = 500; p[2] = 500; p[3] = 500;
# C_T & Delta
p[4] = 20; p[5] = 5;
# Sigmoid Ranges -- (g_Ca, g_K, nHalf)
p[6] = 3; p[7] = 6; p[8] = 20;
# Calcium Handling
A = -1; k = -1/.1;
I_app = 0;
C = 1; g_L = .5;
E_L = -50; E_Ca = 100; E_K = -70;
#################################################
## Initial Condition Testing -- Tests for Multistability
IC_V = collect(-75.0:15.0:75.0);
IC_n = collect(0:.2:1);
IC_gCa = collect(0:p[6]/5:p[6]);
IC_gK = collect(0:p[7]/5:p[7]);
IC_Ca = collect(0.0:50.0:300.0);
IC_nHalf = collect(0:p[8]/5:p[8]);
AllICs = collect(Iterators.product(IC_V,IC_n,IC_gCa,IC_gK,IC_Ca,IC_nHalf));
#################################################
# Define functions
sigmoid(x) = 1/(1+exp(-x));
tau(V, halfN) = 3/cosh(((V-halfN)-10)/29);
Ninf(V, halfN) = sigmoid(((V-halfN)-10)/7.25);
# Define Currents
Ica(g_Ca,V) = g_Ca*(sigmoid((V+1)/7.5)+.1)*(V-E_Ca);
Ik(g_K,n,V) = g_K*n*(V-E_K);
Il(V) = g_L*(V-E_L);
# u[1] - V cell 1, u[2] - slow var cell, u[3] - g_Ca, u[4] - g_K, u[5] - [Ca], u[6] - halfN
# p[1] - tau_L.gCa , p[2] - tau_L.gK, p[3] - tau_L.nHalf, p[4] - C_T, p[5] - Delta, p[6] - g_Ca, p[7] - g_K, p[8] - nHalf
function f(du,u,p,t)
du[1] = (-1*Il(u[1]) + -1*Ica(u[3],u[1]) + -1*Ik(u[4],u[2],u[1]) + I_app)/C
du[2] = (Ninf(u[1],u[6])-u[2])/tau(u[1],u[6])
du[3] = (p[6]*sigmoid((p[4]-u[5])/p[5]) - u[3])/p[1]
du[4] = (p[7]*sigmoid((u[5]-p[4])/p[5]) - u[4])/p[2]
du[5] = -1*k*(A*Ica(u[3],u[1]) - u[5])
du[6] = (p[8]*(sigmoid((p[4]-u[5])/p[5])) - u[6])/p[3]
end
tspan = (0.0f0,8000.0f0);
prob = ODEProblem(f,ICs,tspan,p);
function prob_func(prob,i,repeat)
remake(prob,u0=AllICs[i])
end
ensemble_prob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)
@time sim = solve(ensemble_prob,Tsit5(),EnsembleThreads(),trajectories=1000)
Maybe there’s a better way to do this? I don’t know.