|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
25
|
3 suppressPackageStartupMessages(library("optparse"))
|
|
|
4
|
|
26
|
5 options_list <- list(
|
|
29
|
6 make_option(c("-i", "--input"), action="store", help="Input dataset"),
|
|
22
|
7 make_option(c("-o", "--output"), action="store", help="Output dataset")
|
|
0
|
8 )
|
|
|
9
|
|
13
|
10 parser <- OptionParser(usage="%prog [options] file", options_list)
|
|
|
11 args <- parse_args(parser, positional_arguments=TRUE)
|
|
|
12 opt <- args$options
|
|
|
13
|
|
25
|
14 #########################################
|
|
|
15 daylength=function(L){
|
|
|
16 # from Forsythe 1995
|
|
|
17 p=0.8333
|
|
|
18 dl<-NULL
|
|
|
19 for (i in 1:365) {
|
|
|
20 theta<-0.2163108+2*atan(0.9671396*tan(0.00860*(i-186)))
|
|
|
21 phi<-asin(0.39795*cos(theta))
|
|
|
22 dl[i]<-24-24/pi*acos((sin(p*pi/180)+sin(L*pi/180)*sin(phi))/(cos(L*pi/180)*cos(phi)))
|
|
|
23 }
|
|
|
24 dl # return a vector of daylength in 365 days
|
|
|
25 }
|
|
|
26 #########################################
|
|
13
|
27
|
|
25
|
28 #########################################
|
|
|
29 # source("daylength.R")
|
|
|
30 hourtemp=function(L,date){
|
|
|
31 # L=37.5 specify this in main program
|
|
|
32 threshold<-12.7 # base development threshold for BMSB
|
|
|
33 # threshold2<-threshold/24 degree hour accumulation
|
|
|
34 #expdata<-tempdata[1:365,11:13] # Use daily max, min, mean
|
|
|
35 dnp<-expdata[date,2] # daily minimum
|
|
|
36 dxp<-expdata[date,3] # daily maximum
|
|
|
37 dmean<-0.5*(dnp+dxp)
|
|
|
38 #if (dmean>0) {
|
|
|
39 #dnp<-dnp-k1*dmean
|
|
|
40 #dxp<-dxp+k2*dmean
|
|
|
41 #} else {
|
|
|
42 #dnp<-dnp+k1*dmean
|
|
|
43 #dxp<-dxp-k2*dmean
|
|
|
44 #}
|
|
|
45 dd<-0 # initialize degree day accumulation
|
|
|
46
|
|
|
47 if (dxp<threshold) {dd<-0} else
|
|
|
48 {
|
|
|
49 dlprofile<-daylength(L) # extract daylength data for entire year
|
|
|
50 T<-NULL # initialize hourly temperature
|
|
|
51 dh<-NULL #initialize degree hour vector
|
|
|
52 # date<-200
|
|
|
53 y<-dlprofile[date] # calculate daylength in given date
|
|
|
54 z<-24-y # night length
|
|
|
55 a<-1.86 # lag coefficient
|
|
|
56 b<-2.20 # night coefficient
|
|
|
57 #tempdata<-read.csv("tempdata.csv") #import raw data set
|
|
|
58 # Should be outside function otherwise its redundant
|
|
|
59 risetime<-12-y/2 # sunrise time
|
|
|
60 settime<-12+y/2 # sunset time
|
|
|
61 ts<-(dxp-dnp)*sin(pi*(settime-5)/(y+2*a))+dnp
|
|
|
62 for (i in 1:24){
|
|
|
63 if (i>risetime && i<settime) {
|
|
|
64 m<-i-5 # number of hours after Tmin until sunset
|
|
|
65 T[i]=(dxp-dnp)*sin(pi*m/(y+2*a))+dnp
|
|
|
66 if (T[i]<8.4) {dh[i]<-0} else
|
|
|
67 {dh[i]<-T[i]-8.4}
|
|
|
68 } else
|
|
|
69 if (i>settime){
|
|
|
70 n<-i-settime
|
|
|
71 T[i]=dnp+(ts-dnp)*exp(-b*n/z)
|
|
|
72 if (T[i]<8.4) {dh[i]<-0} else
|
|
|
73 {dh[i]<-T[i]-8.4}
|
|
|
74 } else
|
|
|
75 {
|
|
|
76 n<-i+24-settime
|
|
|
77 T[i]=dnp+(ts-dnp)*exp(-b*n/z)
|
|
|
78 if (T[i]<8.4) {dh[i]<-0} else
|
|
|
79 {dh[i]<-T[i]-8.4}
|
|
|
80 }
|
|
|
81 }
|
|
|
82 dd<-sum(dh)/24
|
|
|
83 }
|
|
|
84 return=c(dmean,dd)
|
|
|
85 return
|
|
|
86 }
|
|
|
87 #########################################
|
|
|
88
|
|
|
89
|
|
|
90 #########################################
|
|
|
91 mortality.egg=function(temperature){
|
|
|
92 if (temperature<12.7) {
|
|
|
93 mort.prob=0.8} else
|
|
|
94 {mort.prob=0.8-temperature/40
|
|
|
95 if (mort.prob<0) {mort.prob=0.01}
|
|
|
96 }
|
|
|
97 return=mort.prob
|
|
|
98 return
|
|
|
99 }
|
|
|
100 #########################################
|
|
|
101
|
|
|
102
|
|
|
103 #########################################
|
|
|
104 mortality.nymph=function(temperature){
|
|
|
105 if (temperature<12.7) {
|
|
|
106 mort.prob=0.03} else
|
|
|
107 {mort.prob=temperature*0.0008+0.03}
|
|
|
108 return=mort.prob
|
|
|
109 return
|
|
|
110 }
|
|
|
111 #########################################
|
|
|
112
|
|
|
113 #########################################
|
|
|
114 mortality.adult=function(temperature){
|
|
|
115 if (temperature<12.7) {
|
|
|
116 mort.prob=0.002} else
|
|
|
117 {mort.prob=temperature*0.0005+0.02}
|
|
|
118 return=mort.prob
|
|
|
119 return
|
|
|
120 }
|
|
|
121 #########################################
|
|
|
122
|
|
|
123 # model initialization
|
|
|
124 # setwd(“/home/lunarmouse/Dropbox/Nelson's project/")
|
|
|
125 # PLEASE CHANGE TO YOUR OWN DIRECTORY!!!
|
|
|
126 # PLEASE LOAD BSMB FUNCTIONS FIRST!!!
|
|
|
127
|
|
|
128 n<-1000 # start with 1000 individuals
|
|
|
129 # Generation, Stage, DD, T, Diapause
|
|
|
130 vec.ini<-c(0,3,0,0,0)
|
|
|
131 # overwintering, previttelogenic,DD=0, T=0, no-diapause
|
|
|
132 vec.mat<-rep(vec.ini,n)
|
|
|
133 vec.mat<-t(matrix(vec.mat,nrow=5)) # complete matrix for the population
|
|
|
134 L<-35.58 # latitude for Asheville NC
|
|
|
135 ph.p<-daylength(L) # complete photoperiod profile in a year, requires daylength function
|
|
|
136
|
|
|
137 #load("asheville2014.Rdat") # load temperature data@location/year
|
|
|
138 load(opt$input) # load temperature data@location/year
|
|
|
139 tot.pop<-NULL # time series of population size
|
|
|
140 gen0.pop<-rep(0,365) # gen.0 pop size
|
|
|
141 gen1.pop<-rep(0,365)
|
|
|
142 gen2.pop<-rep(0,365)
|
|
|
143 S0<-S1<-S2<-S3<-S4<-S5<-rep(0,365)
|
|
|
144 g0.adult<-g1.adult<-g2.adult<-rep(0,365)
|
|
|
145 N.newborn<-N.death<-N.adult<-rep(0,365)
|
|
|
146 dd.day<-rep(0,365)
|
|
|
147
|
|
|
148 ptm <- proc.time() # start tick
|
|
|
149
|
|
|
150 for (day in 1:365) { # all the day
|
|
|
151 photoperiod<-ph.p[day] # photoperiod in the day
|
|
|
152 temp.profile<-hourtemp(L,day)
|
|
|
153 mean.temp<-temp.profile[1]
|
|
|
154 dd.temp<-temp.profile[2]
|
|
|
155 dd.day[day]<-dd.temp
|
|
|
156 death.vec<-NULL # trash bin for death
|
|
|
157 birth.vec<-NULL # new born
|
|
|
158 #n<-length(vec.mat[,1]) # population size at previous day
|
|
|
159
|
|
|
160 for (i in 1:n) { # all individual
|
|
|
161 vec.ind<-vec.mat[i,] # find individual record
|
|
|
162
|
|
|
163 # first of all, still alive?
|
|
|
164 if(vec.ind[2]==0){ # egg
|
|
|
165 death.prob=mortality.egg(mean.temp)
|
|
|
166 } else if (vec.ind[2]==1 | vec.ind[2]==2) {
|
|
|
167 death.prob=mortality.nymph(mean.temp)
|
|
|
168 } else if (vec.ind[2]==3 | vec.ind[2]==4 | vec.ind[2]==5) { # for adult
|
|
|
169 if (day<120 && day>270) {death.prob=0.33*mortality.adult(mean.temp)
|
|
|
170 } else {
|
|
|
171 death.prob=mortality.adult(mean.temp)} # reduce adult mortality after fall equinox
|
|
19
|
172 }
|
|
25
|
173
|
|
|
174
|
|
|
175 #(or dependent on temperature and life stage?)
|
|
|
176 u.d<-runif(1)
|
|
|
177 if (u.d<death.prob) {
|
|
|
178 death.vec<-c(death.vec,i)} else # aggregrate index of dead bug
|
|
|
179 {
|
|
|
180 # event 1 end of diapause
|
|
|
181 if (vec.ind[1]==0 && vec.ind[2]==3) { # overwintering adult (previttelogenic)
|
|
|
182 if (photoperiod>13.5 && vec.ind[3]>77 && day<180) { # add 77C to become fully reproductively matured
|
|
|
183 vec.ind<-c(0,4,0,0,0) # transfer to vittelogenic
|
|
29
|
184 vec.mat[i,]<-vec.ind
|
|
|
185
|
|
25
|
186 } else {
|
|
|
187 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
|
|
|
188 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
|
|
29
|
189 vec.mat[i,]<-vec.ind
|
|
25
|
190 }
|
|
|
191 }
|
|
|
192
|
|
|
193 if (vec.ind[1]!=0 && vec.ind[2]==3) { # NOT overwintering adult (previttelogenic)
|
|
|
194 current.gen<-vec.ind[1]
|
|
|
195 if (vec.ind[3]>77) { # add 77C to become fully reproductively matured
|
|
|
196 vec.ind<-c(current.gen,4,0,0,0) # transfer to vittelogenic
|
|
29
|
197 vec.mat[i,]<-vec.ind
|
|
25
|
198 } else {
|
|
|
199 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
|
|
|
200 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
|
|
29
|
201 vec.mat[i,]<-vec.ind
|
|
25
|
202 }
|
|
|
203 }
|
|
|
204
|
|
|
205
|
|
|
206
|
|
|
207 # event 2 oviposition -- where population dynamics comes from
|
|
|
208 if (vec.ind[2]==4 && vec.ind[1]==0 && mean.temp>10) { # vittelogenic stage, overwintering generation
|
|
|
209 if (vec.ind[4]==0) { # just turned in vittelogenic stage
|
|
|
210 n.birth=round(runif(1,2,8))} else{
|
|
|
211 p.birth=0.01 # daily probability of birth
|
|
|
212 u1<-runif(1)
|
|
|
213 if (u1<p.birth) {n.birth=round(runif(1,2,8))}
|
|
|
214 }
|
|
|
215 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
|
|
|
216 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
|
|
29
|
217 vec.mat[i,]<-vec.ind
|
|
25
|
218 if (n.birth>0) { # add new birth -- might be in different generations
|
|
29
|
219 new.gen<-vec.ind[1]+1 # generation +1
|
|
|
220 new.ind<-c(new.gen,0,0,0,0) # egg profile
|
|
|
221 new.vec<-rep(new.ind,n.birth)
|
|
|
222 new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile
|
|
|
223 birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day
|
|
|
224 }
|
|
25
|
225 }
|
|
|
226
|
|
|
227 # event 2 oviposition -- for gen 1.
|
|
|
228 if (vec.ind[2]==4 && vec.ind[1]==1 && mean.temp>12.5 && day<222) { # vittelogenic stage, 1st generation
|
|
|
229 if (vec.ind[4]==0) { # just turned in vittelogenic stage
|
|
|
230 n.birth=round(runif(1,2,8))} else{
|
|
|
231 p.birth=0.01 # daily probability of birth
|
|
|
232 u1<-runif(1)
|
|
|
233 if (u1<p.birth) {n.birth=round(runif(1,2,8))}
|
|
|
234 }
|
|
|
235 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
|
|
|
236 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
|
|
29
|
237 vec.mat[i,]<-vec.ind
|
|
25
|
238 if (n.birth>0) { # add new birth -- might be in different generations
|
|
29
|
239 new.gen<-vec.ind[1]+1 # generation +1
|
|
|
240 new.ind<-c(new.gen,0,0,0,0) # egg profile
|
|
|
241 new.vec<-rep(new.ind,n.birth)
|
|
|
242 new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile
|
|
|
243 birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day
|
|
|
244 }
|
|
25
|
245 }
|
|
|
246
|
|
|
247
|
|
19
|
248
|
|
25
|
249 # event 3 development (with diapause determination)
|
|
|
250 # event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg)
|
|
|
251 if (vec.ind[2]==0) { # egg stage
|
|
|
252 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
|
|
|
253 if (vec.ind[3]>=68) { # from egg to young nymph, DD requirement met
|
|
29
|
254 current.gen<-vec.ind[1]
|
|
|
255 vec.ind<-c(current.gen,1,0,0,0) # transfer to young nym stage
|
|
|
256 } else {
|
|
|
257 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
|
|
|
258 }
|
|
|
259 vec.mat[i,]<-vec.ind
|
|
25
|
260 }
|
|
|
261
|
|
|
262 # event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause)
|
|
|
263 if (vec.ind[2]==1) { # young nymph stage
|
|
|
264 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
|
|
|
265 if (vec.ind[3]>=250) { # from young to old nymph, DD requirement met
|
|
29
|
266 current.gen<-vec.ind[1]
|
|
|
267 vec.ind<-c(current.gen,2,0,0,0) # transfer to old nym stage
|
|
|
268 if (photoperiod<13.5 && day > 180) {vec.ind[5]<-1} # prepare for diapausing
|
|
|
269 } else {
|
|
|
270 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
|
|
|
271 }
|
|
25
|
272 vec.mat[i,]<-vec.ind
|
|
|
273 }
|
|
|
274
|
|
|
275
|
|
|
276 # event 3.3 old nymph to adult: previttelogenic or diapausing?
|
|
|
277 if (vec.ind[2]==2) { # old nymph stage
|
|
|
278 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
|
|
|
279 if (vec.ind[3]>=200) { # from old to adult, DD requirement met
|
|
|
280 current.gen<-vec.ind[1]
|
|
|
281 if (vec.ind[5]==0) { # non-diapausing adult -- previttelogenic
|
|
|
282 vec.ind<-c(current.gen,3,0,0,0)
|
|
|
283 } else { # diapausing
|
|
|
284 vec.ind<-c(current.gen,5,0,0,1)
|
|
29
|
285 }
|
|
25
|
286 } else {
|
|
|
287 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
|
|
|
288 }
|
|
|
289 vec.mat[i,]<-vec.ind
|
|
|
290 }
|
|
|
291
|
|
|
292 # event 4 growing of diapausing adult (unimportant, but still necessary)##
|
|
|
293 if (vec.ind[2]==5) {
|
|
|
294 vec.ind[3]<-vec.ind[3]+dd.temp
|
|
|
295 vec.ind[4]<-vec.ind[4]+1
|
|
|
296 vec.mat[i,]<-vec.ind
|
|
|
297 }
|
|
|
298
|
|
|
299 } # else if it is still alive
|
|
|
300
|
|
|
301 } # end of the individual bug loop
|
|
|
302
|
|
|
303 # find how many died
|
|
|
304 n.death<-length(death.vec)
|
|
|
305 if (n.death>0) {
|
|
|
306 vec.mat<-vec.mat[-death.vec, ]}
|
|
|
307 # remove record of dead
|
|
|
308 # find how many new born
|
|
|
309 n.newborn<-length(birth.vec[,1])
|
|
|
310 vec.mat<-rbind(vec.mat,birth.vec)
|
|
|
311 # update population size for the next day
|
|
|
312 n<-n-n.death+n.newborn
|
|
|
313
|
|
|
314 # aggregate results by day
|
|
|
315 tot.pop<-c(tot.pop,n)
|
|
|
316 s0<-sum(vec.mat[,2]==0) #egg
|
|
|
317 s1<-sum(vec.mat[,2]==1) # young nymph
|
|
|
318 s2<-sum(vec.mat[,2]==2) # old nymph
|
|
|
319 s3<-sum(vec.mat[,2]==3) # previtellogenic
|
|
|
320 s4<-sum(vec.mat[,2]==4) # vitellogenic
|
|
|
321 s5<-sum(vec.mat[,2]==5) # diapausing
|
|
|
322 gen0<-sum(vec.mat[,1]==0) # overwintering adult
|
|
|
323 gen1<-sum(vec.mat[,1]==1) # first generation
|
|
|
324 gen2<-sum(vec.mat[,1]==2) # second generation
|
|
|
325 n.adult<-sum(vec.mat[,2]==3)+sum(vec.mat[,2]==4)+sum(vec.mat[,2]==5) # sum of all adults
|
|
|
326 gen0.pop[day]<-gen0 # gen.0 pop size
|
|
|
327 gen1.pop[day]<-gen1
|
|
|
328 gen2.pop[day]<-gen2
|
|
|
329 S0[day]<-s0
|
|
|
330 S1[day]<-s1
|
|
|
331 S2[day]<-s2
|
|
|
332 S3[day]<-s3
|
|
|
333 S4[day]<-s4
|
|
|
334 S5[day]<-s5
|
|
|
335 g0.adult[day]<-sum(vec.mat[,1]==0)
|
|
|
336 g1.adult[day]<-sum((vec.mat[,1]==1 & vec.mat[,2]==3) | (vec.mat[,1]==1 & vec.mat[,2]==4) | (vec.mat[,1]==1 & vec.mat[,2]==5))
|
|
|
337 g2.adult[day]<-sum((vec.mat[,1]==2 & vec.mat[,2]==3) | (vec.mat[,1]==2 & vec.mat[,2]==4) | (vec.mat[,1]==2 & vec.mat[,2]==5))
|
|
|
338
|
|
|
339
|
|
|
340
|
|
|
341 N.newborn[day]<-n.newborn
|
|
|
342 N.death[day]<-n.death
|
|
|
343 N.adult[day]<-n.adult
|
|
29
|
344 #print(c(day,n,n.adult))
|
|
25
|
345 }
|
|
|
346
|
|
|
347 proc.time() - ptm
|
|
|
348 dd.cum<-cumsum(dd.day)
|
|
29
|
349 save(dd.day,dd.cum,S0,S1,S2,S3,S4,S5,N.newborn,N.death,N.adult,tot.pop,gen0.pop,gen1.pop,gen2.pop,g0.adult,g1.adult,g2.adult,file=opt$output)
|
