Created
September 27, 2020 07:22
-
-
Save ChrisRackauckas/26b97f963c5f8ca46da19959a9bbbca4 to your computer and use it in GitHub Desktop.
State-dependent delay diffeq model from "Modelling the Effects of Temperature on Temperature Mosquito Seasonal Abundance"
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # -------------------------------------------------------------------------------- | |
| # Ewing model | |
| # translation by: [email protected] (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") |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # -------------------------------------------------------------------------------- | |
| # Ewing model | |
| # translation by: [email protected] (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") |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # -------------------------------------------------------------------------------- | |
| # Ewing model | |
| # translation by: [email protected] (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