Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
State-dependent delay diffeq model from "Modelling the Effects of Temperature on Temperature Mosquito Seasonal Abundance"
# --------------------------------------------------------------------------------
# Ewing model
# translation by: slwu89@berkeleu.edu (July 2020)
# --------------------------------------------------------------------------------
using DifferentialEquations
using Plots
using LabelledArrays
using StaticArrays
# --------------------------------------------------------------------------------
# constant parameters
# --------------------------------------------------------------------------------
# predation on pupae
pupae_pars = (a=1,h=0.002,r=0.001,V=200)
p0 = pupae_pars.r/pupae_pars.h
p1 = pupae_pars.V/(pupae_pars.a*pupae_pars.h)
parvec = @SLVector (
# temp
:phi, # PHASE
:lambda, # A
:mu, # M
:gamma, # POWER
# photoperiod
:L, # latitude
# oviposition
:max_egg, # max egg raft size, R
# gonotrophic cycle
:q1, # KG
:q2, # QG
:q3, # BG
# egg death
:nu_0E, # U3
:nu_1E, # U4
:nu_2E, # U5
# larvae death
:nu_0L, # U3
:nu_1L, # U4
:nu_2L, # U5
# pupae death
:nu_0P, # U3
:nu_1P, # U4
:nu_2P, # U5
# adult death
:alpha_A, # ALPHA
:beta_A, # BETA
# egg maturation
:alpha_E, # ALPHA
:beta_E, # BETA
# larvae maturation
:alpha_L, # ALPHA
:beta_L, # BETA
# pupae maturation
:alpha_P, # ALPHA
:beta_P, # BETA
# predation on pupae
:p0,
:p1
)
parameters = parvec(
1.4, # phi
6.3, # lambda
10.3, # mu
1.21, # gamma
51, # L
200, # max_egg
# gonotrophic cycle
0.2024, # (q1) KG
74.48, # (q2) QG
0.2456, # (q3) BG
# egg death
0.0157, # (nu_0E) U3
20.5, # (nu_1E) U4
7, # (nu_2E) U5
# larvae death
0.0157, # (nu_0L) U3
20.5, # (nu_1L) U4
7, # (nu_2L) U5
# pupae death
0.0157, # (nu_0P) U3
20.5, # (nu_1P) U4
7, # (nu_2P) U5
# adult death
2.166e-8, # (alpha_A) ALPHA
4.483, # (beta_A) BETA
# egg maturation
0.0022, # (alpha_E) ALPHA
1.77, # (beta_E) BETA
# larvae maturation
0.00315, # (alpha_L) ALPHA
1.12, # (beta_L) BETA
# pupae maturation
0.0007109, # (alpha_P) ALPHA
1.8865648, # (beta_P) BETA
# predation on pupae
p0, # p0
p1 # p1
)
# maximum allowed values for rates (same as FORTRAN)
death_max = 1.0 # 1.0
death_min_a = 0.01 # 0.01
gon_min = 0.0333 # 0.0333
maturation_min = 0.016667 # 0.016667
# --------------------------------------------------------------------------------
# external forcing (temp, photoperiod)
# --------------------------------------------------------------------------------
# temperature as modified cosine function
function temperature(t, pars)
phi = pars.phi # PHASE
lambda = pars.lambda # A
mu = pars.mu # M
gamma = pars.gamma # POWER
temp = 0.0
if t < 0.0
temp = (mu - lambda) + lambda * 2.0 * (0.5 * (1.0 + cos(2.0 * pi * (0.0 - phi) / 365.0)))^gamma
else
temp = (mu - lambda) + lambda * 2.0 * (0.5 * (1.0 + cos(2.0 * pi * (t - phi) / 365.0)))^gamma
end
return temp
end
# photoperiod
function daylight(t, pars)
L = pars.L # latitude (51 in thesis)
# define photoperiod values
EPS = asin(0.39795 * cos(0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (t - 3.5)))))
NUM = sin(0.8333 * pi/ 180.0) + (sin(L * pi / 180.0) * sin(EPS))
DEN = cos(L * pi / 180.0) * cos(EPS)
DAYLIGHT = 24.0 - (24.0 / pi) * acos(NUM / DEN)
return DAYLIGHT
end
# --------------------------------------------------------------------------------
# diapause
# --------------------------------------------------------------------------------
# pp: photoperiod
function diapause_spring(pp)
1.0 / (1.0 + exp(5.0 * (14.0 - pp)))
end
function diapause_autumn(pp)
1.0 / (1.0 + exp(5.0 * (13.0 - pp)))
end
# --------------------------------------------------------------------------------
# oviposition
# --------------------------------------------------------------------------------
# per-capita oviposition rate
# d: diapause
# G: duration of gonotrophic cycle
# pars: parameters
function oviposition(d, G, pars)
max_egg = pars.max_egg
egg_raft = d * max_egg * 0.5
ovi = egg_raft / G
return ovi
end
# --------------------------------------------------------------------------------
# mortality
# --------------------------------------------------------------------------------
# egg mortality
function death_egg_rate(temp, pars)
nu_0E = pars.nu_0E # U3
nu_1E = pars.nu_1E # U4
nu_2E = pars.nu_2E # U5
# calculate egg death rate
egg_d = nu_0E * exp(((temp - nu_1E) / nu_2E)^2)
if egg_d > death_max
egg_d = death_max
end
return egg_d
end
# larvae mortality
function death_larvae_rate(temp, pars)
nu_0L = pars.nu_0L # U3
nu_1L = pars.nu_1L # U4
nu_2L = pars.nu_2L # U5
# calculate egg death rate
larvae_d = nu_0L * exp(((temp - nu_1L) / nu_2L)^2)
if larvae_d > death_max
larvae_d = death_max
end
return larvae_d
end
# pupal mortality
function death_pupae_rate(temp, pars)
nu_0P = pars.nu_0P # U3
nu_1P = pars.nu_1P # U4
nu_2P = pars.nu_2P # U5
# calculate egg death rate
pupal_d = nu_0P * exp(((temp - nu_1P)/nu_2P)^2)
if pupal_d > death_max
pupal_d = death_max
end
return pupal_d
end
# adult mortality
function death_adult_rate(temp, pars)
alpha_A = pars.alpha_A # ALPHA
beta_A = pars.beta_A # BETA
# calculate adult death rate
adult_d = alpha_A * (temp^beta_A)
if adult_d < death_min_a
adult_d = death_min_a
end
return adult_d
end
# --------------------------------------------------------------------------------
# development stage delays
# --------------------------------------------------------------------------------
# G: duration of gonotrophic cycle
function gonotrophic(temp, pars)
q1 = pars.q1 # KG
q2 = pars.q2 # QG
q3 = pars.q3 # BG
# calculate gonotrophic cycle length
if temp < 0.0
grate = 0.0333
else
grate = q1 / (1 + q2*exp(-q3*temp))
end
if grate < gon_min
grate = gon_min
end
return 1.0 / grate
end
# g_E
function egg_maturation_rate(temp, pars)
alpha_E = pars.alpha_E # ALPHA
beta_E = pars.beta_E # BETA
# calculate egg development rate
if temp < 0.0
egg_maturation = 0.016667
else
egg_maturation = alpha_E * (temp^beta_E)
end
if egg_maturation < maturation_min
egg_maturation = maturation_min
end
return egg_maturation
end
# g_L
function larvae_maturation_rate(temp, pars)
alpha_L = pars.alpha_L # ALPHA
beta_L = pars.beta_L # BETA
# calculate larvae development rate
if temp < 0.0
larvae_maturation = 0.016667
else
larvae_maturation = alpha_L * (temp^beta_L)
end
if larvae_maturation < maturation_min
larvae_maturation = maturation_min
end
return larvae_maturation
end
# g_P
function pupae_maturation_rate(temp, pars)
alpha_P = pars.alpha_P # ALPHA
beta_P = pars.beta_P # BETA
# calculate larvae development rate
if temp < 0.0
pupae_maturation = 0.016667
else
pupae_maturation = alpha_P * (temp^beta_P)
end
if pupae_maturation < maturation_min
pupae_maturation = maturation_min
end
return pupae_maturation
end
# --------------------------------------------------------------------------------
# state variable history
# --------------------------------------------------------------------------------
function h(p,t; idxs = nothing)
temp = temperature(t,p)
# history vector
Y = zeros(13)
Y[8] = 1.0 / egg_maturation_rate(temp,p) # tau_E
Y[9] = 1.0 / larvae_maturation_rate(temp,p) # tau_L
Y[10] = 1.0 / pupae_maturation_rate(temp,p) # tau_P
Y[5] = exp(-death_egg_rate(temp,p) * Y[8]) # S_E
Y[6] = exp(-death_larvae_rate(temp,p) * Y[9]) # S_L
Y[7] = exp(-death_pupae_rate(temp,p) * Y[10]) #S_P
temp_L = temperature(t - Y[9],p)
temp_P = temperature(t - Y[10],p)
Y[11] = 1.0 / egg_maturation_rate(temp_L,p) # tau_E(t - tau_L(t))
Y[12] = 1.0 / larvae_maturation_rate(temp_P,p) # tau_L(t - tau_P(t))
temp_LP = temperature(t - Y[10] - Y[12],p)
Y[13] = 1.0 / egg_maturation_rate(temp_LP,p) # tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
idxs === nothing ? Y : Y[idxs]
end
# --------------------------------------------------------------------------------
# system of DDEs
# 13 equations
# 6 delays
# --------------------------------------------------------------------------------
# dxdt
# Y(1) = E
# Y(2) = L
# Y(3) = P
# Y(4) = A
# Y(5) = S_E
# Y(6) = S_L
# Y(7) = S_P
# Y(8) = tau_E
# Y(9) = tau_L
# Y(10) = tau_P
# Y(11) = tau_E(t - tau_L(t))
# Y(12) = tau_L(t - tau_P(t))
# Y(13) = tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# delays
# Z(x,1) = t - tau_E(t)
# Z(x,2) = t - tau_L(t) - tau_E(t - tau_L(t))
# Z(x,3) = t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t))
# Z(x,4) = t - tau_L(t)
# Z(x,5) = t - tau_P(t) - tau_L(t - tau_P(t))
# Z(x,6) = t - tau_P(t)
# when you access values of Z in FORTRAN, corresponds to h in Julia
# Z(4,1) = A(t - tau_E(t))
# Z(4,2) = A(t - tau_L(t) - tau_E(t - tau_L(t)))
# Z(5,4) = S_E(t - tau_L(t))
# Z(4,3) = A(t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# Z(5,5) = S_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# Z(6,6) = S_L(t - tau_P(t))
# Z(2,4) = L(t - tau_L(t))
# DDEs
function ewing_dde(du, u, h, p, t)
# state variables
E = u[1]
LAR = u[2]
PUP = u[3]
ADU = u[4]
SE = u[5]
SL = u[6]
SP = u[7]
DE = u[8] # tau_E(t)
DL = u[9] # tau_L(t)
DP = u[10] # tau_P(t)
DEL = u[11] # tau_E(t - tau_L(t))
DLP = u[12] # tau_L(t - tau_P(t))
DELP = u[13] # tau_E(t - tau_P(t) - tau_L(t - taup_P(t)))
# larval predation parameters
p0 = p.p0
p1 = p.p1
# Z: state variables at each of the 6 lagged times (lags follow same order as Z/BETA in DDE_SOLVER)
# Z
Z1 = h(p, t - DE; idxs = 4) # Z(x,1): t - tau_E(t)
Z2 = h(p, t - DL - DEL; idxs = 4) # Z(x,2): t - tau_L(t) - tau_E(t - tau_L(t))
Z3 = h(p, t - DP - DLP - DELP; idxs = 4) # Z(x,3): t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
Z4 = h(p, t - DL) # Z(x,4): t - tau_L(t)
Z5 = h(p, t - DP - DLP) # Z(x,5): t - tau_P(t) - tau_L(t - tau_P(t))
Z6 = h(p, t - DP) # Z(x,6): t - tau_P(t)
# (lagged) temperature
temp = temperature(t, p)
temp_E = temperature(t - DE, p)
temp_L = temperature(t - DL, p)
temp_P = temperature(t - DP, p)
temp_EL = temperature(t - DL - Z4[8], p)
temp_ELP = temperature(t - DP - Z6[9] - Z5[8], p)
temp_LP = temperature(t - DP - Z6[9], p)
# (lagged) photoperiod
pp = daylight(t, p)
pp_1 = daylight(t - 1, p)
pp_E = daylight(t - DE, p)
pp_EL = daylight(t - DL - Z4[8], p)
pp_ELP = daylight(t - DP - Z6[9] - Z5[8], p)
# (lagged) gonotrophic cycle
gon = gonotrophic(temp, p)
gon_E = gonotrophic(temp_E, p)
gon_EL = gonotrophic(temp_EL, p)
gon_ELP = gonotrophic(temp_ELP, p)
# diapause and birth
if pp > pp_1
dia = diapause_spring(pp)
dia_E = diapause_spring(pp_E)
dia_EL = diapause_spring(pp_EL)
dia_ELP = diapause_spring(pp_ELP)
else
dia = diapause_autumn(pp)
dia_E = diapause_autumn(pp_E)
dia_EL = diapause_autumn(pp_EL)
dia_ELP = diapause_autumn(pp_ELP)
end
birth = oviposition(dia,gon,p)
birth_E = oviposition(dia_E,gon_E,p)
birth_EL = oviposition(dia_EL,gon_EL,p)
birth_ELP = oviposition(dia_ELP,gon_ELP,p)
# (lagged) death
death_egg = death_egg_rate(temp,p)
death_egg_E = death_egg_rate(temp_E,p)
death_larvae = death_larvae_rate(temp,p)
death_larvae_L = death_larvae_rate(temp_L,p)
death_pupae = death_pupae_rate(temp,p)
death_pupae_P = death_pupae_rate(temp_P,p)
death_adult = death_adult_rate(temp,p)
# (lagged) development
larvae_maturation = larvae_maturation_rate(temp,p)
larvae_maturation_L = larvae_maturation_rate(temp_L,p)
larvae_maturation_P = larvae_maturation_rate(temp_P,p)
larvae_maturation_LP = larvae_maturation_rate(temp_LP,p)
egg_maturation = egg_maturation_rate(temp,p)
egg_maturation_E = egg_maturation_rate(temp_E,p)
egg_maturation_L = egg_maturation_rate(temp_L,p)
egg_maturation_EL = egg_maturation_rate(temp_EL,p)
egg_maturation_LP = egg_maturation_rate(temp_LP,p)
egg_maturation_ELP = egg_maturation_rate(temp_ELP,p)
pupae_maturation = pupae_maturation_rate(temp,p)
pupae_maturation_P = pupae_maturation_rate(temp_P,p)
# DDEs describing change in state duration
dDEdt = 1 - egg_maturation/egg_maturation_E
dDLdt = 1 - larvae_maturation/larvae_maturation_L
dDPdt = 1 - pupae_maturation/pupae_maturation_P
dDELdt = (1 - dDLdt) * (1 - egg_maturation_L/egg_maturation_EL)
dDLPdt = (1 - dDPdt) * (1 - larvae_maturation_P/larvae_maturation_LP)
dDELPdt = (1 - dDPdt - dDLPdt) * (1 - egg_maturation_LP/egg_maturation_ELP)
# stage recruitment
R_E = birth * ADU
R_L = birth_E * Z1 * SE * egg_maturation/egg_maturation_E
R_P = birth_EL * Z2 * Z4[5] * SL * larvae_maturation/larvae_maturation_L * (1 - dDELdt)
R_A = birth_ELP * Z3 * Z5[5] * Z6[6] * SP * pupae_maturation/pupae_maturation_P * (1 - dDLPdt) * (1 - dDELPdt)
# maturation rates
M_E = R_L
M_L = R_P
M_P = R_A
# death rates
D_E = death_egg * E
D_L = ((p0*LAR/(p1+LAR)) + death_larvae) * LAR
D_P = death_pupae * PUP
D_A = death_adult * ADU
# DDE system
du_1 = R_E - M_E - D_E # E
du_2 = R_L - M_L - D_L # L
du_3 = R_P - M_P - D_P # P
du_4 = R_A - D_A # A
du_5 = SE * ((egg_maturation * death_egg_E / egg_maturation_E) - death_egg)
du_6 = SL * (((p0*Z4[2] / (p1+Z4[2])) + death_larvae_L) * (1-dDLdt) - (p0*LAR / (p1+LAR)) - death_larvae)
du_7 = SP * ((pupae_maturation * death_pupae_P / pupae_maturation_P) - death_pupae)
du_8 = dDEdt # tau_E(t)
du_9 = dDLdt # tau_L(t)
du_10 = dDPdt # tau_P(t)
du_11 = dDELdt # tau_E(t - tau_L(t))
du_12 = dDLPdt # tau_L(t - tau_P(t))
du_13 = dDELPdt # tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
du[1] = du_1
du[2] = du_2
du[3] = du_3
du[4] = du_4
du[5] = du_5
du[6] = du_6
du[7] = du_7
du[8] = du_8
du[9] = du_9
du[10] = du_10
du[11] = du_11
du[12] = du_12
du[13] = du_13
end
# --------------------------------------------------------------------------------
# initial conditions
# --------------------------------------------------------------------------------
# A0: A(0)
# t0: temp(0); assumed constant for t<0
function calculate_IC(A0, t0, pars)
u0 = zeros(13)
u0[4] = A0
# calculate initial lags first
u0[8] = 1.0 / egg_maturation_rate(t0,pars) # tau_E
u0[9] = 1.0 / larvae_maturation_rate(t0,pars) # tau_L
u0[10] = 1.0 / pupae_maturation_rate(t0,pars) # tau_P
u0[11] = u0[8] # tau_E(t - tau_L(t))
u0[12] = u0[9] # tau_L(t - tau_P(t))
u0[13] = u0[8] # tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# survival probabilities
u0[5] = exp(-u0[8]*death_egg_rate(t0,pars)) # S_E
u0[6] = exp(-u0[9]*death_larvae_rate(t0,pars)) # S_L
u0[7] = exp(-u0[10]*death_pupae_rate(t0,pars)) # S_P
return u0
end
# --------------------------------------------------------------------------------
# lag functions for state dependent delays
# --------------------------------------------------------------------------------
# dependent lags follow same order as DDE_SOLVER "BETA" subroutine
# t - tau_E(t)
function deplag_1(u,p,t)
return u[8]
end
# t - tau_L(t) - tau_E(t - tau_L(t))
function deplag_2(u,p,t)
return u[9] + u[11]
end
# t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t))
function deplag_3(u,p,t)
return u[10] + u[12] + u[13]
end
# t - tau_L(t)
function deplag_4(u,p,t)
return u[9]
end
# t - tau_P(t) - tau_L(t - tau_P(t))
function deplag_5(u,p,t)
return u[10] + u[12]
end
# t - tau_P(t)
function deplag_6(u,p,t)
return u[10]
end
# --------------------------------------------------------------------------------
# simulation with default parameters
# --------------------------------------------------------------------------------
A0 = 12000.0
t0 = 0.0
temp0 = temperature(t0,parameters)
u0 = calculate_IC(A0,temp0,parameters)
t0 = 0.0
times = (t0,t0+365.0*10)
prob = DDEProblem{true}(ewing_dde,u0,h,times,parameters)
solv_alg = MethodOfSteps(Tsit5())
@time soln = solve(prob,solv_alg,abstol=1e-12,reltol=1e-12)
using Profile
@profile soln = solve(prob,solv_alg,abstol=1e-12,reltol=1e-12,dtmax=1.0)
Juno.profiler()
Profile.clear()
# Start
# 62.497892 seconds (679.36 M allocations: 99.183 GiB, 26.59% gc time)
# No tracking
# 10.120807 seconds (130.40 M allocations: 13.106 GiB, 22.32% gc time)
# Tsit5
# 7.837309 seconds (96.16 M allocations: 7.480 GiB, 18.22% gc time)
# No dtmax cap and some scalar interpolation
# 7.266140 seconds (86.18 M allocations: 5.582 GiB, 15.41% gc time)
# in-place
# 5.205283 seconds (68.59 M allocations: 2.103 GiB, 8.82% gc time)
# life stage (1 plot)
plot(
soln,
vars=[1,2,3,4],
title="Life Stages",
legend=:topleft,
labels=["E" "L" "P" "A"]
)
savefig("julia_ewing_stages.pdf")
# survival functions
plot(
soln,
vars=[5,6,7],
title="Survival",
legend=:topleft,
labels=["E" "L" "P"]
)
savefig("julia_ewing_surv.pdf")
# delays
plot(
soln,vars=[8,9,10,11,12,13],
title="Delays",
legend=:bottomleft,labels=["tau_E" "tau_L" "tau_P" "tau_EL" "tau_LP" "tau_ELP"]
)
savefig("julia_ewing_delays.pdf")
# --------------------------------------------------------------------------------
# Ewing model
# translation by: slwu89@berkeleu.edu (July 2020)
# --------------------------------------------------------------------------------
using DifferentialEquations
using Plots
using LabelledArrays
using StaticArrays
# --------------------------------------------------------------------------------
# constant parameters
# --------------------------------------------------------------------------------
# predation on pupae
pupae_pars = (a=1,h=0.002,r=0.001,V=200)
p0 = pupae_pars.r/pupae_pars.h
p1 = pupae_pars.V/(pupae_pars.a*pupae_pars.h)
parvec = @SLVector (
# temp
:phi, # PHASE
:lambda, # A
:mu, # M
:gamma, # POWER
# photoperiod
:L, # latitude
# oviposition
:max_egg, # max egg raft size, R
# gonotrophic cycle
:q1, # KG
:q2, # QG
:q3, # BG
# egg death
:nu_0E, # U3
:nu_1E, # U4
:nu_2E, # U5
# larvae death
:nu_0L, # U3
:nu_1L, # U4
:nu_2L, # U5
# pupae death
:nu_0P, # U3
:nu_1P, # U4
:nu_2P, # U5
# adult death
:alpha_A, # ALPHA
:beta_A, # BETA
# egg maturation
:alpha_E, # ALPHA
:beta_E, # BETA
# larvae maturation
:alpha_L, # ALPHA
:beta_L, # BETA
# pupae maturation
:alpha_P, # ALPHA
:beta_P, # BETA
# predation on pupae
:p0,
:p1
)
parameters = parvec(
1.4, # phi
6.3, # lambda
10.3, # mu
1.21, # gamma
51, # L
200, # max_egg
# gonotrophic cycle
0.2024, # (q1) KG
74.48, # (q2) QG
0.2456, # (q3) BG
# egg death
0.0157, # (nu_0E) U3
20.5, # (nu_1E) U4
7, # (nu_2E) U5
# larvae death
0.0157, # (nu_0L) U3
20.5, # (nu_1L) U4
7, # (nu_2L) U5
# pupae death
0.0157, # (nu_0P) U3
20.5, # (nu_1P) U4
7, # (nu_2P) U5
# adult death
2.166e-8, # (alpha_A) ALPHA
4.483, # (beta_A) BETA
# egg maturation
0.0022, # (alpha_E) ALPHA
1.77, # (beta_E) BETA
# larvae maturation
0.00315, # (alpha_L) ALPHA
1.12, # (beta_L) BETA
# pupae maturation
0.0007109, # (alpha_P) ALPHA
1.8865648, # (beta_P) BETA
# predation on pupae
p0, # p0
p1 # p1
)
# maximum allowed values for rates (same as FORTRAN)
death_max = 1.0 # 1.0
death_min_a = 0.01 # 0.01
gon_min = 0.0333 # 0.0333
maturation_min = 0.016667 # 0.016667
# --------------------------------------------------------------------------------
# external forcing (temp, photoperiod)
# --------------------------------------------------------------------------------
# temperature as modified cosine function
function temperature(t, pars)
phi = pars.phi # PHASE
lambda = pars.lambda # A
mu = pars.mu # M
gamma = pars.gamma # POWER
temp = 0.0
if t < 0.0
temp = (mu - lambda) + lambda * 2.0 * (0.5 * (1.0 + cos(2.0 * pi * (0.0 - phi) / 365.0)))^gamma
else
temp = (mu - lambda) + lambda * 2.0 * (0.5 * (1.0 + cos(2.0 * pi * (t - phi) / 365.0)))^gamma
end
return temp
end
# photoperiod
function daylight(t, pars)
L = pars.L # latitude (51 in thesis)
# define photoperiod values
EPS = asin(0.39795 * cos(0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (t - 3.5)))))
NUM = sin(0.8333 * pi/ 180.0) + (sin(L * pi / 180.0) * sin(EPS))
DEN = cos(L * pi / 180.0) * cos(EPS)
DAYLIGHT = 24.0 - (24.0 / pi) * acos(NUM / DEN)
return DAYLIGHT
end
# --------------------------------------------------------------------------------
# diapause
# --------------------------------------------------------------------------------
# pp: photoperiod
function diapause_spring(pp)
1.0 / (1.0 + exp(5.0 * (14.0 - pp)))
end
function diapause_autumn(pp)
1.0 / (1.0 + exp(5.0 * (13.0 - pp)))
end
# --------------------------------------------------------------------------------
# oviposition
# --------------------------------------------------------------------------------
# per-capita oviposition rate
# d: diapause
# G: duration of gonotrophic cycle
# pars: parameters
function oviposition(d, G, pars)
max_egg = pars.max_egg
egg_raft = d * max_egg * 0.5
ovi = egg_raft / G
return ovi
end
# --------------------------------------------------------------------------------
# mortality
# --------------------------------------------------------------------------------
# egg mortality
function death_egg_rate(temp, pars)
nu_0E = pars.nu_0E # U3
nu_1E = pars.nu_1E # U4
nu_2E = pars.nu_2E # U5
# calculate egg death rate
egg_d = nu_0E * exp(((temp - nu_1E) / nu_2E)^2)
if egg_d > death_max
egg_d = death_max
end
return egg_d
end
# larvae mortality
function death_larvae_rate(temp, pars)
nu_0L = pars.nu_0L # U3
nu_1L = pars.nu_1L # U4
nu_2L = pars.nu_2L # U5
# calculate egg death rate
larvae_d = nu_0L * exp(((temp - nu_1L) / nu_2L)^2)
if larvae_d > death_max
larvae_d = death_max
end
return larvae_d
end
# pupal mortality
function death_pupae_rate(temp, pars)
nu_0P = pars.nu_0P # U3
nu_1P = pars.nu_1P # U4
nu_2P = pars.nu_2P # U5
# calculate egg death rate
pupal_d = nu_0P * exp(((temp - nu_1P)/nu_2P)^2)
if pupal_d > death_max
pupal_d = death_max
end
return pupal_d
end
# adult mortality
function death_adult_rate(temp, pars)
alpha_A = pars.alpha_A # ALPHA
beta_A = pars.beta_A # BETA
# calculate adult death rate
adult_d = alpha_A * (temp^beta_A)
if adult_d < death_min_a
adult_d = death_min_a
end
return adult_d
end
# --------------------------------------------------------------------------------
# development stage delays
# --------------------------------------------------------------------------------
# G: duration of gonotrophic cycle
function gonotrophic(temp, pars)
q1 = pars.q1 # KG
q2 = pars.q2 # QG
q3 = pars.q3 # BG
# calculate gonotrophic cycle length
if temp < 0.0
grate = 0.0333
else
grate = q1 / (1 + q2*exp(-q3*temp))
end
if grate < gon_min
grate = gon_min
end
return 1.0 / grate
end
# g_E
function egg_maturation_rate(temp, pars)
alpha_E = pars.alpha_E # ALPHA
beta_E = pars.beta_E # BETA
# calculate egg development rate
if temp < 0.0
egg_maturation = 0.016667
else
egg_maturation = alpha_E * (temp^beta_E)
end
if egg_maturation < maturation_min
egg_maturation = maturation_min
end
return egg_maturation
end
# g_L
function larvae_maturation_rate(temp, pars)
alpha_L = pars.alpha_L # ALPHA
beta_L = pars.beta_L # BETA
# calculate larvae development rate
if temp < 0.0
larvae_maturation = 0.016667
else
larvae_maturation = alpha_L * (temp^beta_L)
end
if larvae_maturation < maturation_min
larvae_maturation = maturation_min
end
return larvae_maturation
end
# g_P
function pupae_maturation_rate(temp, pars)
alpha_P = pars.alpha_P # ALPHA
beta_P = pars.beta_P # BETA
# calculate larvae development rate
if temp < 0.0
pupae_maturation = 0.016667
else
pupae_maturation = alpha_P * (temp^beta_P)
end
if pupae_maturation < maturation_min
pupae_maturation = maturation_min
end
return pupae_maturation
end
# --------------------------------------------------------------------------------
# state variable history
# --------------------------------------------------------------------------------
function h(p,t)
temp = temperature(t,p)
# history vector
Y = zeros(13)
Y[8] = 1.0 / egg_maturation_rate(temp,p) # tau_E
Y[9] = 1.0 / larvae_maturation_rate(temp,p) # tau_L
Y[10] = 1.0 / pupae_maturation_rate(temp,p) # tau_P
Y[5] = exp(-death_egg_rate(temp,p) * Y[8]) # S_E
Y[6] = exp(-death_larvae_rate(temp,p) * Y[9]) # S_L
Y[7] = exp(-death_pupae_rate(temp,p) * Y[10]) #S_P
temp_L = temperature(t - Y[9],p)
temp_P = temperature(t - Y[10],p)
Y[11] = 1.0 / egg_maturation_rate(temp_L,p) # tau_E(t - tau_L(t))
Y[12] = 1.0 / larvae_maturation_rate(temp_P,p) # tau_L(t - tau_P(t))
temp_LP = temperature(t - Y[10] - Y[12],p)
Y[13] = 1.0 / egg_maturation_rate(temp_LP,p) # tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
return Y
end
# --------------------------------------------------------------------------------
# system of DDEs
# 13 equations
# 6 delays
# --------------------------------------------------------------------------------
# dxdt
# Y(1) = E
# Y(2) = L
# Y(3) = P
# Y(4) = A
# Y(5) = S_E
# Y(6) = S_L
# Y(7) = S_P
# Y(8) = tau_E
# Y(9) = tau_L
# Y(10) = tau_P
# Y(11) = tau_E(t - tau_L(t))
# Y(12) = tau_L(t - tau_P(t))
# Y(13) = tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# delays
# Z(x,1) = t - tau_E(t)
# Z(x,2) = t - tau_L(t) - tau_E(t - tau_L(t))
# Z(x,3) = t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t))
# Z(x,4) = t - tau_L(t)
# Z(x,5) = t - tau_P(t) - tau_L(t - tau_P(t))
# Z(x,6) = t - tau_P(t)
# when you access values of Z in FORTRAN, corresponds to h in Julia
# Z(4,1) = A(t - tau_E(t))
# Z(4,2) = A(t - tau_L(t) - tau_E(t - tau_L(t)))
# Z(5,4) = S_E(t - tau_L(t))
# Z(4,3) = A(t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# Z(5,5) = S_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# Z(6,6) = S_L(t - tau_P(t))
# Z(2,4) = L(t - tau_L(t))
# DDEs
function ewing_dde(u, h, p, t)
# state variables
E = u[1]
LAR = u[2]
PUP = u[3]
ADU = u[4]
SE = u[5]
SL = u[6]
SP = u[7]
DE = u[8] # tau_E(t)
DL = u[9] # tau_L(t)
DP = u[10] # tau_P(t)
DEL = u[11] # tau_E(t - tau_L(t))
DLP = u[12] # tau_L(t - tau_P(t))
DELP = u[13] # tau_E(t - tau_P(t) - tau_L(t - taup_P(t)))
# larval predation parameters
p0 = p.p0
p1 = p.p1
# Z: state variables at each of the 6 lagged times (lags follow same order as Z/BETA in DDE_SOLVER)
# Z
Z1 = h(p, t - DE) # Z(x,1): t - tau_E(t)
Z2 = h(p, t - DL - DEL) # Z(x,2): t - tau_L(t) - tau_E(t - tau_L(t))
Z3 = h(p, t - DP - DLP - DELP) # Z(x,3): t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
Z4 = h(p, t - DL) # Z(x,4): t - tau_L(t)
Z5 = h(p, t - DP - DLP) # Z(x,5): t - tau_P(t) - tau_L(t - tau_P(t))
Z6 = h(p, t - DP) # Z(x,6): t - tau_P(t)
# (lagged) temperature
temp = temperature(t, p)
temp_E = temperature(t - DE, p)
temp_L = temperature(t - DL, p)
temp_P = temperature(t - DP, p)
temp_EL = temperature(t - DL - Z4[8], p)
temp_ELP = temperature(t - DP - Z6[9] - Z5[8], p)
temp_LP = temperature(t - DP - Z6[9], p)
# (lagged) photoperiod
pp = daylight(t, p)
pp_1 = daylight(t - 1, p)
pp_E = daylight(t - DE, p)
pp_EL = daylight(t - DL - Z4[8], p)
pp_ELP = daylight(t - DP - Z6[9] - Z5[8], p)
# (lagged) gonotrophic cycle
gon = gonotrophic(temp, p)
gon_E = gonotrophic(temp_E, p)
gon_EL = gonotrophic(temp_EL, p)
gon_ELP = gonotrophic(temp_ELP, p)
# diapause and birth
if pp > pp_1
dia = diapause_spring(pp)
dia_E = diapause_spring(pp_E)
dia_EL = diapause_spring(pp_EL)
dia_ELP = diapause_spring(pp_ELP)
else
dia = diapause_autumn(pp)
dia_E = diapause_autumn(pp_E)
dia_EL = diapause_autumn(pp_EL)
dia_ELP = diapause_autumn(pp_ELP)
end
birth = oviposition(dia,gon,p)
birth_E = oviposition(dia_E,gon_E,p)
birth_EL = oviposition(dia_EL,gon_EL,p)
birth_ELP = oviposition(dia_ELP,gon_ELP,p)
# (lagged) death
death_egg = death_egg_rate(temp,p)
death_egg_E = death_egg_rate(temp_E,p)
death_larvae = death_larvae_rate(temp,p)
death_larvae_L = death_larvae_rate(temp_L,p)
death_pupae = death_pupae_rate(temp,p)
death_pupae_P = death_pupae_rate(temp_P,p)
death_adult = death_adult_rate(temp,p)
# (lagged) development
larvae_maturation = larvae_maturation_rate(temp,p)
larvae_maturation_L = larvae_maturation_rate(temp_L,p)
larvae_maturation_P = larvae_maturation_rate(temp_P,p)
larvae_maturation_LP = larvae_maturation_rate(temp_LP,p)
egg_maturation = egg_maturation_rate(temp,p)
egg_maturation_E = egg_maturation_rate(temp_E,p)
egg_maturation_L = egg_maturation_rate(temp_L,p)
egg_maturation_EL = egg_maturation_rate(temp_EL,p)
egg_maturation_LP = egg_maturation_rate(temp_LP,p)
egg_maturation_ELP = egg_maturation_rate(temp_ELP,p)
pupae_maturation = pupae_maturation_rate(temp,p)
pupae_maturation_P = pupae_maturation_rate(temp_P,p)
# DDEs describing change in state duration
dDEdt = 1 - egg_maturation/egg_maturation_E
dDLdt = 1 - larvae_maturation/larvae_maturation_L
dDPdt = 1 - pupae_maturation/pupae_maturation_P
dDELdt = (1 - dDLdt) * (1 - egg_maturation_L/egg_maturation_EL)
dDLPdt = (1 - dDPdt) * (1 - larvae_maturation_P/larvae_maturation_LP)
dDELPdt = (1 - dDPdt - dDLPdt) * (1 - egg_maturation_LP/egg_maturation_ELP)
# stage recruitment
R_E = birth * ADU
R_L = birth_E * Z1[4] * SE * egg_maturation/egg_maturation_E
R_P = birth_EL * Z2[4] * Z4[5] * SL * larvae_maturation/larvae_maturation_L * (1 - dDELdt)
R_A = birth_ELP * Z3[4] * Z5[5] * Z6[6] * SP * pupae_maturation/pupae_maturation_P * (1 - dDLPdt) * (1 - dDELPdt)
# maturation rates
M_E = R_L
M_L = R_P
M_P = R_A
# death rates
D_E = death_egg * E
D_L = ((p0*LAR/(p1+LAR)) + death_larvae) * LAR
D_P = death_pupae * PUP
D_A = death_adult * ADU
# DDE system
du_1 = R_E - M_E - D_E # E
du_2 = R_L - M_L - D_L # L
du_3 = R_P - M_P - D_P # P
du_4 = R_A - D_A # A
du_5 = SE * ((egg_maturation * death_egg_E / egg_maturation_E) - death_egg)
du_6 = SL * (((p0*Z4[2] / (p1+Z4[2])) + death_larvae_L) * (1-dDLdt) - (p0*LAR / (p1+LAR)) - death_larvae)
du_7 = SP * ((pupae_maturation * death_pupae_P / pupae_maturation_P) - death_pupae)
du_8 = dDEdt # tau_E(t)
du_9 = dDLdt # tau_L(t)
du_10 = dDPdt # tau_P(t)
du_11 = dDELdt # tau_E(t - tau_L(t))
du_12 = dDLPdt # tau_L(t - tau_P(t))
du_13 = dDELPdt # tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
@SVector [du_1,du_2,du_3,du_4,du_5,du_6,du_7,du_8,du_9,du_10,du_11,du_12,du_13]
end
# --------------------------------------------------------------------------------
# initial conditions
# --------------------------------------------------------------------------------
# A0: A(0)
# t0: temp(0); assumed constant for t<0
function calculate_IC(A0, t0, pars)
u0 = zeros(13)
u0[4] = A0
# calculate initial lags first
u0[8] = 1.0 / egg_maturation_rate(t0,pars) # tau_E
u0[9] = 1.0 / larvae_maturation_rate(t0,pars) # tau_L
u0[10] = 1.0 / pupae_maturation_rate(t0,pars) # tau_P
u0[11] = u0[8] # tau_E(t - tau_L(t))
u0[12] = u0[9] # tau_L(t - tau_P(t))
u0[13] = u0[8] # tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# survival probabilities
u0[5] = exp(-u0[8]*death_egg_rate(t0,pars)) # S_E
u0[6] = exp(-u0[9]*death_larvae_rate(t0,pars)) # S_L
u0[7] = exp(-u0[10]*death_pupae_rate(t0,pars)) # S_P
return u0
end
# --------------------------------------------------------------------------------
# lag functions for state dependent delays
# --------------------------------------------------------------------------------
# dependent lags follow same order as DDE_SOLVER "BETA" subroutine
# t - tau_E(t)
function deplag_1(u,p,t)
return u[8]
end
# t - tau_L(t) - tau_E(t - tau_L(t))
function deplag_2(u,p,t)
return u[9] + u[11]
end
# t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t))
function deplag_3(u,p,t)
return u[10] + u[12] + u[13]
end
# t - tau_L(t)
function deplag_4(u,p,t)
return u[9]
end
# t - tau_P(t) - tau_L(t - tau_P(t))
function deplag_5(u,p,t)
return u[10] + u[12]
end
# t - tau_P(t)
function deplag_6(u,p,t)
return u[10]
end
# --------------------------------------------------------------------------------
# simulation with default parameters
# --------------------------------------------------------------------------------
A0 = 12000.0
t0 = 0.0
temp0 = temperature(t0,parameters)
u0 = calculate_IC(A0,temp0,parameters)
t0 = 0.0
times = (t0,t0+365.0*10)
prob = DDEProblem(ewing_dde,u0,h,times,parameters; dependent_lags = (deplag_1,deplag_2,deplag_3,deplag_4,deplag_5,deplag_6))
solv_alg = MethodOfSteps(Vern6())
soln = solve(prob,solv_alg,abstol=1e-12,reltol=1e-12,dtmax=1.0)
# life stage (1 plot)
plot(
soln,
vars=[1,2,3,4],
title="Life Stages",
legend=:topleft,
labels=["E" "L" "P" "A"]
)
savefig("julia_ewing_stages.pdf")
# survival functions
plot(
soln,
vars=[5,6,7],
title="Survival",
legend=:topleft,
labels=["E" "L" "P"]
)
savefig("julia_ewing_surv.pdf")
# delays
plot(
soln,vars=[8,9,10,11,12,13],
title="Delays",
legend=:bottomleft,labels=["tau_E" "tau_L" "tau_P" "tau_EL" "tau_LP" "tau_ELP"]
)
savefig("julia_ewing_delays.pdf")
# --------------------------------------------------------------------------------
# Ewing model
# translation by: slwu89@berkeleu.edu (July 2020)
# --------------------------------------------------------------------------------
using DifferentialEquations
using Plots
using LabelledArrays
using StaticArrays
# --------------------------------------------------------------------------------
# constant parameters
# --------------------------------------------------------------------------------
# predation on pupae
pupae_pars = (a=1,h=0.002,r=0.001,V=200)
p0 = pupae_pars.r/pupae_pars.h
p1 = pupae_pars.V/(pupae_pars.a*pupae_pars.h)
parvec = @SLVector (
# temp
:phi, # PHASE
:lambda, # A
:mu, # M
:gamma, # POWER
# photoperiod
:L, # latitude
# oviposition
:max_egg, # max egg raft size, R
# gonotrophic cycle
:q1, # KG
:q2, # QG
:q3, # BG
# egg death
:nu_0E, # U3
:nu_1E, # U4
:nu_2E, # U5
# larvae death
:nu_0L, # U3
:nu_1L, # U4
:nu_2L, # U5
# pupae death
:nu_0P, # U3
:nu_1P, # U4
:nu_2P, # U5
# adult death
:alpha_A, # ALPHA
:beta_A, # BETA
# egg maturation
:alpha_E, # ALPHA
:beta_E, # BETA
# larvae maturation
:alpha_L, # ALPHA
:beta_L, # BETA
# pupae maturation
:alpha_P, # ALPHA
:beta_P, # BETA
# predation on pupae
:p0,
:p1
)
parameters = parvec(
1.4, # phi
6.3, # lambda
10.3, # mu
1.21, # gamma
51, # L
200, # max_egg
# gonotrophic cycle
0.2024, # (q1) KG
74.48, # (q2) QG
0.2456, # (q3) BG
# egg death
0.0157, # (nu_0E) U3
20.5, # (nu_1E) U4
7, # (nu_2E) U5
# larvae death
0.0157, # (nu_0L) U3
20.5, # (nu_1L) U4
7, # (nu_2L) U5
# pupae death
0.0157, # (nu_0P) U3
20.5, # (nu_1P) U4
7, # (nu_2P) U5
# adult death
2.166e-8, # (alpha_A) ALPHA
4.483, # (beta_A) BETA
# egg maturation
0.0022, # (alpha_E) ALPHA
1.77, # (beta_E) BETA
# larvae maturation
0.00315, # (alpha_L) ALPHA
1.12, # (beta_L) BETA
# pupae maturation
0.0007109, # (alpha_P) ALPHA
1.8865648, # (beta_P) BETA
# predation on pupae
p0, # p0
p1 # p1
)
# maximum allowed values for rates (same as FORTRAN)
death_max = 1.0 # 1.0
death_min_a = 0.01 # 0.01
gon_min = 0.0333 # 0.0333
maturation_min = 0.016667 # 0.016667
# --------------------------------------------------------------------------------
# external forcing (temp, photoperiod)
# --------------------------------------------------------------------------------
# temperature as modified cosine function
function temperature(t, pars)
phi = pars.phi # PHASE
lambda = pars.lambda # A
mu = pars.mu # M
gamma = pars.gamma # POWER
temp = 0.0
if t < 0.0
temp = (mu - lambda) + lambda * 2.0 * (0.5 * (1.0 + cos(2.0 * pi * (0.0 - phi) / 365.0)))^gamma
else
temp = (mu - lambda) + lambda * 2.0 * (0.5 * (1.0 + cos(2.0 * pi * (t - phi) / 365.0)))^gamma
end
return temp
end
# photoperiod
function daylight(t, pars)
L = pars.L # latitude (51 in thesis)
# define photoperiod values
EPS = asin(0.39795 * cos(0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (t - 3.5)))))
NUM = sin(0.8333 * pi/ 180.0) + (sin(L * pi / 180.0) * sin(EPS))
DEN = cos(L * pi / 180.0) * cos(EPS)
DAYLIGHT = 24.0 - (24.0 / pi) * acos(NUM / DEN)
return DAYLIGHT
end
# --------------------------------------------------------------------------------
# diapause
# --------------------------------------------------------------------------------
# pp: photoperiod
function diapause_spring(pp)
1.0 / (1.0 + exp(5.0 * (14.0 - pp)))
end
function diapause_autumn(pp)
1.0 / (1.0 + exp(5.0 * (13.0 - pp)))
end
# --------------------------------------------------------------------------------
# oviposition
# --------------------------------------------------------------------------------
# per-capita oviposition rate
# d: diapause
# G: duration of gonotrophic cycle
# pars: parameters
function oviposition(d, G, pars)
max_egg = pars.max_egg
egg_raft = d * max_egg * 0.5
ovi = egg_raft / G
return ovi
end
# --------------------------------------------------------------------------------
# mortality
# --------------------------------------------------------------------------------
# egg mortality
function death_egg_rate(temp, pars)
nu_0E = pars.nu_0E # U3
nu_1E = pars.nu_1E # U4
nu_2E = pars.nu_2E # U5
# calculate egg death rate
egg_d = nu_0E * exp(((temp - nu_1E) / nu_2E)^2)
if egg_d > death_max
egg_d = death_max
end
return egg_d
end
# larvae mortality
function death_larvae_rate(temp, pars)
nu_0L = pars.nu_0L # U3
nu_1L = pars.nu_1L # U4
nu_2L = pars.nu_2L # U5
# calculate egg death rate
larvae_d = nu_0L * exp(((temp - nu_1L) / nu_2L)^2)
if larvae_d > death_max
larvae_d = death_max
end
return larvae_d
end
# pupal mortality
function death_pupae_rate(temp, pars)
nu_0P = pars.nu_0P # U3
nu_1P = pars.nu_1P # U4
nu_2P = pars.nu_2P # U5
# calculate egg death rate
pupal_d = nu_0P * exp(((temp - nu_1P)/nu_2P)^2)
if pupal_d > death_max
pupal_d = death_max
end
return pupal_d
end
# adult mortality
function death_adult_rate(temp, pars)
alpha_A = pars.alpha_A # ALPHA
beta_A = pars.beta_A # BETA
# calculate adult death rate
adult_d = alpha_A * (temp^beta_A)
if adult_d < death_min_a
adult_d = death_min_a
end
return adult_d
end
# --------------------------------------------------------------------------------
# development stage delays
# --------------------------------------------------------------------------------
# G: duration of gonotrophic cycle
function gonotrophic(temp, pars)
q1 = pars.q1 # KG
q2 = pars.q2 # QG
q3 = pars.q3 # BG
# calculate gonotrophic cycle length
if temp < 0.0
grate = 0.0333
else
grate = q1 / (1 + q2*exp(-q3*temp))
end
if grate < gon_min
grate = gon_min
end
return 1.0 / grate
end
# g_E
function egg_maturation_rate(temp, pars)
alpha_E = pars.alpha_E # ALPHA
beta_E = pars.beta_E # BETA
# calculate egg development rate
if temp < 0.0
egg_maturation = 0.016667
else
egg_maturation = alpha_E * (temp^beta_E)
end
if egg_maturation < maturation_min
egg_maturation = maturation_min
end
return egg_maturation
end
# g_L
function larvae_maturation_rate(temp, pars)
alpha_L = pars.alpha_L # ALPHA
beta_L = pars.beta_L # BETA
# calculate larvae development rate
if temp < 0.0
larvae_maturation = 0.016667
else
larvae_maturation = alpha_L * (temp^beta_L)
end
if larvae_maturation < maturation_min
larvae_maturation = maturation_min
end
return larvae_maturation
end
# g_P
function pupae_maturation_rate(temp, pars)
alpha_P = pars.alpha_P # ALPHA
beta_P = pars.beta_P # BETA
# calculate larvae development rate
if temp < 0.0
pupae_maturation = 0.016667
else
pupae_maturation = alpha_P * (temp^beta_P)
end
if pupae_maturation < maturation_min
pupae_maturation = maturation_min
end
return pupae_maturation
end
# --------------------------------------------------------------------------------
# state variable history
# --------------------------------------------------------------------------------
function h(p,t)
temp = temperature(t,p)
# history vector
Y = zeros(13)
Y[8] = 1.0 / egg_maturation_rate(temp,p) # tau_E
Y[9] = 1.0 / larvae_maturation_rate(temp,p) # tau_L
Y[10] = 1.0 / pupae_maturation_rate(temp,p) # tau_P
Y[5] = exp(-death_egg_rate(temp,p) * Y[8]) # S_E
Y[6] = exp(-death_larvae_rate(temp,p) * Y[9]) # S_L
Y[7] = exp(-death_pupae_rate(temp,p) * Y[10]) #S_P
temp_L = temperature(t - Y[9],p)
temp_P = temperature(t - Y[10],p)
Y[11] = 1.0 / egg_maturation_rate(temp_L,p) # tau_E(t - tau_L(t))
Y[12] = 1.0 / larvae_maturation_rate(temp_P,p) # tau_L(t - tau_P(t))
temp_LP = temperature(t - Y[10] - Y[12],p)
Y[13] = 1.0 / egg_maturation_rate(temp_LP,p) # tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
return Y
end
# --------------------------------------------------------------------------------
# system of DDEs
# 13 equations
# 6 delays
# --------------------------------------------------------------------------------
# dxdt
# Y(1) = E
# Y(2) = L
# Y(3) = P
# Y(4) = A
# Y(5) = S_E
# Y(6) = S_L
# Y(7) = S_P
# Y(8) = tau_E
# Y(9) = tau_L
# Y(10) = tau_P
# Y(11) = tau_E(t - tau_L(t))
# Y(12) = tau_L(t - tau_P(t))
# Y(13) = tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# delays
# Z(x,1) = t - tau_E(t)
# Z(x,2) = t - tau_L(t) - tau_E(t - tau_L(t))
# Z(x,3) = t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t))
# Z(x,4) = t - tau_L(t)
# Z(x,5) = t - tau_P(t) - tau_L(t - tau_P(t))
# Z(x,6) = t - tau_P(t)
# when you access values of Z in FORTRAN, corresponds to h in Julia
# Z(4,1) = A(t - tau_E(t))
# Z(4,2) = A(t - tau_L(t) - tau_E(t - tau_L(t)))
# Z(5,4) = S_E(t - tau_L(t))
# Z(4,3) = A(t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# Z(5,5) = S_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# Z(6,6) = S_L(t - tau_P(t))
# Z(2,4) = L(t - tau_L(t))
# DDEs
function ewing_dde(u, h, p, t)
# state variables
E = u[1]
LAR = u[2]
PUP = u[3]
ADU = u[4]
SE = u[5]
SL = u[6]
SP = u[7]
DE = u[8] # tau_E(t)
DL = u[9] # tau_L(t)
DP = u[10] # tau_P(t)
DEL = u[11] # tau_E(t - tau_L(t))
DLP = u[12] # tau_L(t - tau_P(t))
DELP = u[13] # tau_E(t - tau_P(t) - tau_L(t - taup_P(t)))
# larval predation parameters
p0 = p.p0
p1 = p.p1
# Z: state variables at each of the 6 lagged times (lags follow same order as Z/BETA in DDE_SOLVER)
# Z
Z1 = h(p, t - DE) # Z(x,1): t - tau_E(t)
Z2 = h(p, t - DL - DEL) # Z(x,2): t - tau_L(t) - tau_E(t - tau_L(t))
Z3 = h(p, t - DP - DLP - DELP) # Z(x,3): t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
Z4 = h(p, t - DL) # Z(x,4): t - tau_L(t)
Z5 = h(p, t - DP - DLP) # Z(x,5): t - tau_P(t) - tau_L(t - tau_P(t))
Z6 = h(p, t - DP) # Z(x,6): t - tau_P(t)
# (lagged) temperature
temp = temperature(t, p)
temp_E = temperature(t - DE, p)
temp_L = temperature(t - DL, p)
temp_P = temperature(t - DP, p)
temp_EL = temperature(t - DL - Z4[8], p)
temp_ELP = temperature(t - DP - Z6[9] - Z5[8], p)
temp_LP = temperature(t - DP - Z6[9], p)
# (lagged) photoperiod
pp = daylight(t, p)
pp_1 = daylight(t - 1, p)
pp_E = daylight(t - DE, p)
pp_EL = daylight(t - DL - Z4[8], p)
pp_ELP = daylight(t - DP - Z6[9] - Z5[8], p)
# (lagged) gonotrophic cycle
gon = gonotrophic(temp, p)
gon_E = gonotrophic(temp_E, p)
gon_EL = gonotrophic(temp_EL, p)
gon_ELP = gonotrophic(temp_ELP, p)
# diapause and birth
if pp > pp_1
dia = diapause_spring(pp)
dia_E = diapause_spring(pp_E)
dia_EL = diapause_spring(pp_EL)
dia_ELP = diapause_spring(pp_ELP)
else
dia = diapause_autumn(pp)
dia_E = diapause_autumn(pp_E)
dia_EL = diapause_autumn(pp_EL)
dia_ELP = diapause_autumn(pp_ELP)
end
birth = oviposition(dia,gon,p)
birth_E = oviposition(dia_E,gon_E,p)
birth_EL = oviposition(dia_EL,gon_EL,p)
birth_ELP = oviposition(dia_ELP,gon_ELP,p)
# (lagged) death
death_egg = death_egg_rate(temp,p)
death_egg_E = death_egg_rate(temp_E,p)
death_larvae = death_larvae_rate(temp,p)
death_larvae_L = death_larvae_rate(temp_L,p)
death_pupae = death_pupae_rate(temp,p)
death_pupae_P = death_pupae_rate(temp_P,p)
death_adult = death_adult_rate(temp,p)
# (lagged) development
larvae_maturation = larvae_maturation_rate(temp,p)
larvae_maturation_L = larvae_maturation_rate(temp_L,p)
larvae_maturation_P = larvae_maturation_rate(temp_P,p)
larvae_maturation_LP = larvae_maturation_rate(temp_LP,p)
egg_maturation = egg_maturation_rate(temp,p)
egg_maturation_E = egg_maturation_rate(temp_E,p)
egg_maturation_L = egg_maturation_rate(temp_L,p)
egg_maturation_EL = egg_maturation_rate(temp_EL,p)
egg_maturation_LP = egg_maturation_rate(temp_LP,p)
egg_maturation_ELP = egg_maturation_rate(temp_ELP,p)
pupae_maturation = pupae_maturation_rate(temp,p)
pupae_maturation_P = pupae_maturation_rate(temp_P,p)
# DDEs describing change in state duration
dDEdt = 1 - egg_maturation/egg_maturation_E
dDLdt = 1 - larvae_maturation/larvae_maturation_L
dDPdt = 1 - pupae_maturation/pupae_maturation_P
dDELdt = (1 - dDLdt) * (1 - egg_maturation_L/egg_maturation_EL)
dDLPdt = (1 - dDPdt) * (1 - larvae_maturation_P/larvae_maturation_LP)
dDELPdt = (1 - dDPdt - dDLPdt) * (1 - egg_maturation_LP/egg_maturation_ELP)
# stage recruitment
R_E = birth * ADU
R_L = birth_E * Z1[4] * SE * egg_maturation/egg_maturation_E
R_P = birth_EL * Z2[4] * Z4[5] * SL * larvae_maturation/larvae_maturation_L * (1 - dDELdt)
R_A = birth_ELP * Z3[4] * Z5[5] * Z6[6] * SP * pupae_maturation/pupae_maturation_P * (1 - dDLPdt) * (1 - dDELPdt)
# maturation rates
M_E = R_L
M_L = R_P
M_P = R_A
# death rates
D_E = death_egg * E
D_L = ((p0*LAR/(p1+LAR)) + death_larvae) * LAR
D_P = death_pupae * PUP
D_A = death_adult * ADU
# DDE system
du_1 = R_E - M_E - D_E # E
du_2 = R_L - M_L - D_L # L
du_3 = R_P - M_P - D_P # P
du_4 = R_A - D_A # A
du_5 = SE * ((egg_maturation * death_egg_E / egg_maturation_E) - death_egg)
du_6 = SL * (((p0*Z4[2] / (p1+Z4[2])) + death_larvae_L) * (1-dDLdt) - (p0*LAR / (p1+LAR)) - death_larvae)
du_7 = SP * ((pupae_maturation * death_pupae_P / pupae_maturation_P) - death_pupae)
du_8 = dDEdt # tau_E(t)
du_9 = dDLdt # tau_L(t)
du_10 = dDPdt # tau_P(t)
du_11 = dDELdt # tau_E(t - tau_L(t))
du_12 = dDLPdt # tau_L(t - tau_P(t))
du_13 = dDELPdt # tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
@SVector [du_1,du_2,du_3,du_4,du_5,du_6,du_7,du_8,du_9,du_10,du_11,du_12,du_13]
end
# --------------------------------------------------------------------------------
# initial conditions
# --------------------------------------------------------------------------------
# A0: A(0)
# t0: temp(0); assumed constant for t<0
function calculate_IC(A0, t0, pars)
u0 = zeros(13)
u0[4] = A0
# calculate initial lags first
u0[8] = 1.0 / egg_maturation_rate(t0,pars) # tau_E
u0[9] = 1.0 / larvae_maturation_rate(t0,pars) # tau_L
u0[10] = 1.0 / pupae_maturation_rate(t0,pars) # tau_P
u0[11] = u0[8] # tau_E(t - tau_L(t))
u0[12] = u0[9] # tau_L(t - tau_P(t))
u0[13] = u0[8] # tau_E(t - tau_P(t) - tau_L(t - tau_P(t)))
# survival probabilities
u0[5] = exp(-u0[8]*death_egg_rate(t0,pars)) # S_E
u0[6] = exp(-u0[9]*death_larvae_rate(t0,pars)) # S_L
u0[7] = exp(-u0[10]*death_pupae_rate(t0,pars)) # S_P
return u0
end
# --------------------------------------------------------------------------------
# lag functions for state dependent delays
# --------------------------------------------------------------------------------
# dependent lags follow same order as DDE_SOLVER "BETA" subroutine
# t - tau_E(t)
function deplag_1(u,p,t)
return u[8]
end
# t - tau_L(t) - tau_E(t - tau_L(t))
function deplag_2(u,p,t)
return u[9] + u[11]
end
# t - tau_P(t) - tau_L(t - tau_P(t)) - tau_E(t - tau_P(t) - tau_L(t - tau_P(t))
function deplag_3(u,p,t)
return u[10] + u[12] + u[13]
end
# t - tau_L(t)
function deplag_4(u,p,t)
return u[9]
end
# t - tau_P(t) - tau_L(t - tau_P(t))
function deplag_5(u,p,t)
return u[10] + u[12]
end
# t - tau_P(t)
function deplag_6(u,p,t)
return u[10]
end
# --------------------------------------------------------------------------------
# simulation with default parameters
# --------------------------------------------------------------------------------
A0 = 12000.0
t0 = 0.0
temp0 = temperature(t0,parameters)
u0 = calculate_IC(A0,temp0,parameters)
t0 = 0.0
times = (t0,t0+365.0*10)
prob = DDEProblem(ewing_dde,u0,h,times,parameters; dependent_lags = (deplag_1,deplag_2,deplag_3,deplag_4,deplag_5,deplag_6))
solv_alg = MethodOfSteps(Vern6())
soln = solve(prob,solv_alg,abstol=1e-12,reltol=1e-12,dtmax=1.0)
# life stage (1 plot)
plot(
soln,
vars=[1,2,3,4],
title="Life Stages",
legend=:topleft,
labels=["E" "L" "P" "A"]
)
savefig("julia_ewing_stages.pdf")
# survival functions
plot(
soln,
vars=[5,6,7],
title="Survival",
legend=:topleft,
labels=["E" "L" "P"]
)
savefig("julia_ewing_surv.pdf")
# delays
plot(
soln,vars=[8,9,10,11,12,13],
title="Delays",
legend=:bottomleft,labels=["tau_E" "tau_L" "tau_P" "tau_EL" "tau_LP" "tau_ELP"]
)
savefig("julia_ewing_delays.pdf")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment