I am having a hard time figuring out how to declare dependent_lags that depend
on one of the state variables in my DDE model (using DelayDiffEq.jl). It seems
like it should be easy since the state variable is tracking the length of the
delay as it changes with time. So I feel like it should be
d_lags = ((u,p,t) -> u,) prob = DDEProblem(DynDDE,u0,h,tspan; constant_lags = tauM, dependent_lags = d_lags)
However, when I do this the simulation becomes very slow (usually asks for larger maxiters), and I thought declaring lags was supposed to improve performance.
If I don’t declare the dependent lags, the result is reasonable, fast, and seems to be essentially
identical to when they are declared.
I’m fairly new to Julia so hopefully I’ve just made a rookie mistake that someone
can spot. It’s also possible that I’ve misunderstood exactly what is happening when
lags are declared.
I’ve provided below the simplest working version of the problem as I could without
having to spend a lot of time modifying the original. Sorry it is so long. Briefly,
it models a host-parasite system in which the length of the developmental delay period for one
of the larval stages depends on environmental temperatures, which in turn vary
continuously with time (e.g.,shorter development time in the summer, longer in
the winter). In the original, many of the other parameters are functions of time
as well, and the issue above becomes even more pronounced.
I’d greatly appreciate any help!
using QuadGK using Roots using DelayDiffEq using Plots; #Parameters const Temp0K = 15+273.15; const k = 8.62 * 10^(-5); const H = 0.075;const S = 4000000.0 const b_H = 1/5292.5; const b_S = 0.0037; const lambda = 488.0;const mu_W = 0.0186 const Phi_a = 0.024; const Phi_E = 0.46; const mu_P = 1.5e-7; const mu_M = 0.0056 const tauM = 166.82; const mu_U = 0.0001; const mu_I = mu_U; const tauU_a = 25.1; const tauU_E = 0.87; const tauU_EL = 5.0*tauU_E; const tauU_TL = 282.0218; const Alpha = 1.4752e-6;const phi_H = 2.7728e-7;const c_K = 263.0; const d_K = 22.6; const t0 = 105.6; #A couple functions that describe how parameters are sensitive to temprature #Temperature as a function of time Temp(t) = c_K + d_K * sin((t-t0)*2*pi/365); TempKstop = 273.15-0.0; #Development rate, as a function of Temperature function g_U(T) if T > TempKstop::Float64 g_U = 1/(tauU_a * exp((tauU_E/k)*((1/T)-(1/Temp0K))) * (1 + exp((tauU_EL/k)*((1/T)-(1/tauU_TL))))); #Development rate else g_U = 1/(tauU_a * exp((tauU_E/k)*((1/TempKstop::Float64)-(1/Temp0K))) * (1 + exp((tauU_EL/k)*((1/TempKstop::Float64)-(1/tauU_TL)))));#Development rate if it's below freezing end end gt(t) = g_U(Temp(t)) #Infection rate, as a function of temperature function Phi(T) if T > TempKstop::Float64 Phi = Phi_a * exp((-Phi_E/k)*((1/T)-(1/Temp0K))); #Infection rate else Phi = Phi_a * exp((-Phi_E/k)*((1/TempKstop::Float64)-(1/Temp0K))); #Infection rate if it's too cold end return(Phi) end function DynDDE(du,u,h,p,t) #State variables #W,U,I,M,and P track different stages of a parasite life-cycle #DU1 gives the proportion of U that survive the variabl developmental delay period #DM1 give proportion of M that survive constant delay period #tauU tracks the length of the developmental delay, which changes continuously with time W = u; U = u; I = u; M = u; P = u; DU1 = u; DM1 = u; tauU = u; #History calculations W_past = h(p, t-tauU); #W at time t-tauU U_past = h(p, t-tauU); #U at time t-tauU I_past = h(p, t-tauM); #I at time t-tauM P_past = h(p, t-tauM); #P at time t-tauM #"Constant" Survival functions DM2 = exp(-mu_M*tauM)*exp(-b_H*tauM); DU2 = exp(-b_S*tauU)*exp(-phi_H*H*tauU); #Model du = lambda*P - mu_W*W - Phi(Temp(t))*S*W; #W) du = Phi(Temp(t))*S*W - U*(b_S + phi_H*H) - mu_U*U - Phi(Temp(t-tauU))*S*W_past*DU1*DU2*(g_U(Temp(t))/g_U(Temp(t-tauU))); #U) du = Phi(Temp(t-tauU))*S*W_past*DU1*DU2*(g_U(Temp(t))/g_U(Temp(t-tauU))) - I*(mu_I+b_S) - phi_H*H*I; #I) du = phi_H*H*I - M*(mu_M+b_H) - Alpha*P*M/H - phi_H*H*I_past*DM1*DM2; #M) du = phi_H*H*I_past*DM1*DM2 - P*(mu_P+b_H) - Alpha*P; #P) du = DU1*(mu_U*U_past/S*(g_U(Temp(t))/g_U(Temp(t-tauU))) - mu_U*U/S); #DU1) du = DM1*(Alpha*P_past/H - Alpha*P/H); #DM1) du = 1-(g_U(Temp(t))/g_U(Temp(t-tauU))); #tauU) Change in delay duration end #Set up simulation function DySimulate(tstart,tfin) birthday = tstart; gint(u) = quadgk(gt, birthday-u, birthday) #integrate over gt from days[jj]-u to days[jj] initTauU = [fzero(u -> gint(u) - p, [0.001, 5000]) for p in [1.0]]; #Initial lag duration #Initial survival probabilities mUint = exp(-mu_U*initTauU); bSint = exp(-b_S*initTauU); pHint = exp(-phi_H*H*initTauU); DUinit = mUint*bSint*pHint; DUhist = DUinit; DM2 = exp(-mu_M*tauM)*exp(-b_H*tauM); #Initial values u0 = [0.0,0.0,0.0,0.0,1.0,DUinit,DM2,initTauU]; h(p,t) = [0.0,0.0,0.0,0.0,0.0,DUhist,DM2,initTauU]; #Empty system when t<0 tspan = (birthday,tfin); #Declare dependent lags. What should they be? It works well when d_lags is not declared. d_lags = ((u,p,t) -> u,) prob = DDEProblem(DynDDE,u0,h,tspan; constant_lags = tauM, dependent_lags = d_lags) solDynDDE = solve(prob,alg; isoutofdomain=(u,p,t) -> any(x -> x < 0.0, u),abstol = 1e-6,reltol = 1e-4) end alg = MethodOfSteps(Rosenbrock23()); out = DySimulate(1.0,365.0*10) plot(out, vars = (0,[1,2,3,4,5])) plot(out, vars = (0,)) #Visualize change in delay period over time (dtauU/dt)