I am trying to implement a discrete stochastic simulation using tools from DifferentialEquations.jl. It works great for earlier time-points but does not do any jumps later on, inexplicably.
using DifferentialEquations
using Plots
#initial conditions
far1_ic = 0
cln_ic = 100
ste12_ic = 1000
far1cln_ic = 0
u0 = [far1_ic, cln_ic, ste12_ic, far1cln_ic]
#paramaters
#Far1
ks_F = 10 #1
kd_F = 0.01 #2
kd_F_C = 0.01 #3
kA_C = 200 #4
ks_F_S = 200 #5
#Cln
ks_C = 10 #6
kd_C = 0.1 #7
kA_F = 500 #8
ks_C_C = 100 #9
#Ste12
ks_S_alpha = 80 #10
ks_S = 1 #11
ks_S_S = 100 #12
kd_S = 0.01 #13
kA_S = 500 #14
kd_S_F = 0.005 #15
#Far1Cln
ks_FC = 100 #16
kd_FC = 0.05 #17
krev_FC = 600 #18
#misc
alpha = 17 #19
power_alpha = 1 #20
alpha_halfmax = 50 #21
p = [ks_F, kd_F, kd_F_C, kA_C, ks_F_S, ks_C, kd_C, kA_F, ks_C_C,
ks_S_alpha, ks_S, ks_S_S, kd_S, kA_S, kd_S_F, ks_FC, kd_FC,
krev_FC, alpha, power_alpha, alpha_halfmax];
s_far1(u,p,t) = (p[5] * ((u[3]^2)/(u[3]^2+p[14]^2)) + p[1]) * (p[19]^p[20]/(p[19]^p[20]+p[21]^p[20]))
d_far1(u,p,t) = p[2]*u[1] + p[3]*u[1]*u[2]
s_cln(u,p,t) = p[6] + p[9]*((u[2]^2)/(u[2]^2+p[4]^2))
d_cln(u,p,t) = p[7]*u[2]
s_ste12(u,p,t) = p[11]+(p[10]+p[12]*(u[3]^2/(u[3]^2+p[14]^2)))*(p[19]^p[20]/(p[19]^p[20]+p[21]^p[20]))
d_ste12(u,p,t) = p[13]*u[3] + p[15]*u[1]*u[3]
s_far1cln(u,p,t) = p[16]*u[1]*u[2]
d_far1cln(u,p,t) = p[17]*u[4]
r_far1cln(u,p,t) = p[18]*u[4]
function syn_f(integrator)
integrator.u[1] += 1
end
function deg_f(integrator)
integrator.u[1] -= 1
end
function syn_c(integrator)
integrator.u[2] += 1
end
function deg_c(integrator)
integrator.u[2] -= 1
end
function syn_s(integrator)
integrator.u[3] += 1
end
function deg_s(integrator)
integrator.u[3] -= 1
end
function syn_fc(integrator)
integrator.u[4] += 1
integrator.u[1] -= 1
integrator.u[2] -= 1
end
function deg_fc(integrator)
integrator.u[4] -= 1
end
function rev_fc(integrator)
integrator.u[4] -= 1
integrator.u[1] += 1
integrator.u[2] += 1
end
jump1 = ConstantRateJump(s_far1, syn_f)
jump2 = ConstantRateJump(d_far1, deg_f)
jump3 = ConstantRateJump(s_cln, syn_c)
jump4 = ConstantRateJump(d_cln, deg_c)
jump5 = ConstantRateJump(s_ste12, syn_s)
jump6 = ConstantRateJump(d_ste12, deg_s)
jump7 = ConstantRateJump(s_far1cln, syn_fc)
jump8 = ConstantRateJump(d_far1cln, deg_fc)
jump9 = ConstantRateJump(r_far1cln, rev_fc)
p[19] = 27 # reset parameter alpha
tspan = (0.0,600.0) #time span
prob = DiscreteProblem(u0,tspan,p) # variables must be discrete. Provide IC, span, paramaters.
jump_prob = JumpProblem(prob,DirectFW(),jump1,jump2,jump3,jump4,jump5,jump6,jump7,jump8,jump9,save_positions=(false,false))
sol = solve(jump_prob,SSAStepper(),saveat=1);
I have tried saving all positions instead of only at every integer. I have tried using Direct instead of DirectFW. I have tried using default instead of SSAStepper. I checked the ODESolution object and it really does just repeat the same values at each timestep later on; it's not just the plot.
Really appreciate any advice!
This was a bug that has been fixed. Run ]up and your code will work.
