.Rhistory

"Hosts (susceptible)", #4
"Hosts (infected)", #5
"Total hosts") #6
# select output to plot
out_name <- out_names[5]
# choose your death value you want to plot at the end
death_access <- 0.9
### -------------------------- run code from here to plot
#pdf(paste0(getwd(),"/npsi_model_plot_",ttl,".pdf"),onefile=T,width=10,height=8,paper="a4r")
# then run this part to plot in your live R session
layout(matrix(c(1,2,3), 1, 3, byrow = TRUE)) # set plot window
wastes <- list(out_tibble_wh,out_tibble_wd,out_tibble_ws)
ttl_list <- c("Host waste","Drool waste","Summed waste")
pf <- 1 # plot feature counter
pc <- 1 # plot number counter
colvec<-brewer.pal(length(beta_pars),"RdYlGn") # !!! make sure length of colour palette is > beta_pars
# dp <- 0.2
# bp <- 0.2
for(out_tibble in wastes){
for(beta_access in beta_pars){
outplot <- filter(out_tibble, death == death_access & beta == beta_access)
outplot <- outplot$outs ; outplot <- as.data.frame(outplot) # clean output
outplot$"Total host population" <- outplot[,"S"] + outplot[,"I"] # add sum host population
colnames(outplot) <- out_names
if(out_name==names(outplot[2]) | out_name==names(outplot[3])){ylim=c(0,200)}else{ylim=c(0,50)};ylim # set proper ylim
if(pf==length(beta_pars)+1){pf <- 1} # reset pf to reset colour palette for new plot
plot(outplot[,1],outplot[,out_name],type="l",las=1,bty="n",
xlab="Time (years)",ylab=out_name,
col=colvec[pf],#lty=pf,
ylim=ylim
)
par(new=T)
#dev.off()
pf <- pf + 1
} # end bbb loop
pc <- pc + 1
title(paste0(out_name,"\nbeta=",min(beta_pars)," to ",max(beta_pars),", death = ",death_access,"\n",ttl_list[pc-1]))
text(years-40,ylim[2],paste0("b=",beta_pars[1]),col=colvec[1]) # low beta
text(years-10,ylim[2],paste0("b=",max(beta_pars)),col=colvec[length(colvec)]) # high beta
} # end wastes loop
out_names <- c("Time", #1
"Nutrient biomass", #2
"Plant biomass", #3
"Hosts (susceptible)", #4
"Hosts (infected)", #5
"Total hosts") #6
# select output to plot
out_name <- out_names[2]
# choose your death value you want to plot at the end
death_access <- 0.1
# pdf(paste0(getwd(),"/npsi_model_plot_",ttl,".pdf"),onefile=T,width=10,height=8,paper="a4r")
# then run this part to plot in your live R session
layout(matrix(c(1,2,3), 1, 3, byrow = TRUE)) # set plot window
wastes <- list(out_tibble_wh,out_tibble_wd,out_tibble_ws)
ttl_list <- c("Host waste","Drool waste","Summed waste")
pf <- 1 # plot feature counter
pc <- 1 # plot number counter
colvec<-brewer.pal(length(beta_pars),"RdYlGn") # !!! make sure length of colour palette is > beta_pars
# dp <- 0.2
# bp <- 0.2
for(out_tibble in wastes){
for(beta_access in beta_pars){
outplot <- filter(out_tibble, death == death_access & beta == beta_access)
outplot <- outplot$outs ; outplot <- as.data.frame(outplot) # clean output
outplot$"Total host population" <- outplot[,"S"] + outplot[,"I"] # add sum host population
colnames(outplot) <- out_names
if(out_name==names(outplot[2]) | out_name==names(outplot[3])){ylim=c(0,200)}else{ylim=c(0,100)};ylim # set proper ylim
if(pf==length(beta_pars)+1){pf <- 1} # reset pf to reset colour palette for new plot
plot(outplot[,1],outplot[,out_name],type="l",las=1,bty="n",
xlab="Time (years)",ylab=out_name,
col=colvec[pf],#lty=pf,
ylim=ylim
)
par(new=T)
#dev.off()
pf <- pf + 1
} # end bbb loop
pc <- pc + 1
title(paste0(out_name,"\nbeta=",min(beta_pars)," to ",max(beta_pars),", death = ",death_access,"\n",ttl_list[pc-1]))
text(years-40,ylim[2],paste0("b=",beta_pars[1]),col=colvec[1]) # low beta
text(years-10,ylim[2],paste0("b=",max(beta_pars)),col=colvec[length(colvec)]) # high beta
} # end wastes loop
out_name <- out_names[2]
# choose your death value you want to plot at the end
death_access <- 0.1
### -------------------------- run code from here to plot
#pdf(paste0(getwd(),"/npsi_model_plot_",ttl,".pdf"),onefile=T,width=10,height=8,paper="a4r")
# then run this part to plot in your live R session
layout(matrix(c(1,2,3), 1, 3, byrow = TRUE)) # set plot window
wastes <- list(out_tibble_wh,out_tibble_wd,out_tibble_ws)
ttl_list <- c("Host waste","Drool waste","Summed waste")
pf <- 1 # plot feature counter
pc <- 1 # plot number counter
colvec<-brewer.pal(length(beta_pars),"RdYlGn") # !!! make sure length of colour palette is > beta_pars
# dp <- 0.2
# bp <- 0.2
for(out_tibble in wastes){
for(beta_access in beta_pars){
outplot <- filter(out_tibble, death == death_access & beta == beta_access)
outplot <- outplot$outs ; outplot <- as.data.frame(outplot) # clean output
outplot$"Total host population" <- outplot[,"S"] + outplot[,"I"] # add sum host population
colnames(outplot) <- out_names
if(out_name==names(outplot[2]) | out_name==names(outplot[3])){ylim=c(0,200)}else{ylim=c(0,50)};ylim # set proper ylim
if(pf==length(beta_pars)+1){pf <- 1} # reset pf to reset colour palette for new plot
plot(outplot[,1],outplot[,out_name],type="l",las=1,bty="n",
xlab="Time (years)",ylab=out_name,
col=colvec[pf],#lty=pf,
ylim=ylim
)
par(new=T)
#dev.off()
pf <- pf + 1
} # end bbb loop
pc <- pc + 1
title(paste0(out_name,"\nbeta=",min(beta_pars)," to ",max(beta_pars),", death = ",death_access,"\n",ttl_list[pc-1]))
text(years-40,ylim[2],paste0("b=",beta_pars[1]),col=colvec[1]) # low beta
text(years-10,ylim[2],paste0("b=",max(beta_pars)),col=colvec[length(colvec)]) # high beta
} # end wastes loop
beta_access <- 0.9 # choose your beta value you want to plot at the end
death_access <- 0.1 # choose your death value you want to plot at the end
colvv <- "orange" # choose your plot line colour
# then run this part to plot in your live R session
outplot <- filter(out_tibble, death == death_access & beta == beta_access)
outplot <- outplot$outs ; outplot <- as.data.frame(outplot) # clean output
outplot$"Total host population" <- outplot[,"S"] + outplot[,"I"] # add sum host population
layout(matrix(c(1,2,3,4,5,5), 2, 3, byrow = TRUE)) # set plot window
colnames(outplot) <- c("Time", #1
"Nutrient biomass", #2
"Plant biomass", #3
"Hosts (susceptible)", #4
"Hosts (infected)", #5
"Total hosts") #6
for (name in names(outplot)[c(3:5,2,6)]){ # start plot for product, suscep hosts, infec hosts, nutrients, and total hosts
plot(outplot[,1],outplot[,name],type="l",las=1,bty="n",
xlab="Time (years)",ylab=name,col=colvv,
ylim=c(0,round_any(max(outplot[,name]),10,ceiling))
)
title(paste0(name,"\nbeta = ",beta_access," , death = ",death_access))
} # end plot
out_name <- out_names[5]
# choose your death value you want to plot at the end
death_access <- 0.1
# pdf(paste0(getwd(),"/npsi_model_plot_",ttl,".pdf"),onefile=T,width=10,height=8,paper="a4r")
# then run this part to plot in your live R session
layout(matrix(c(1,2,3), 1, 3, byrow = TRUE)) # set plot window
wastes <- list(out_tibble_wh,out_tibble_wd,out_tibble_ws)
ttl_list <- c("Host waste","Drool waste","Summed waste")
pf <- 1 # plot feature counter
pc <- 1 # plot number counter
colvec<-brewer.pal(length(beta_pars),"RdYlGn") # !!! make sure length of colour palette is > beta_pars
# dp <- 0.2
# bp <- 0.2
for(out_tibble in wastes){
for(beta_access in beta_pars){
outplot <- filter(out_tibble, death == death_access & beta == beta_access)
outplot <- outplot$outs ; outplot <- as.data.frame(outplot) # clean output
outplot$"Total host population" <- outplot[,"S"] + outplot[,"I"] # add sum host population
colnames(outplot) <- out_names
if(out_name==names(outplot[2]) | out_name==names(outplot[3])){ylim=c(0,200)}else{ylim=c(0,100)};ylim # set proper ylim
if(pf==length(beta_pars)+1){pf <- 1} # reset pf to reset colour palette for new plot
plot(outplot[,1],outplot[,out_name],type="l",las=1,bty="n",
xlab="Time (years)",ylab=out_name,
col=colvec[pf],#lty=pf,
ylim=ylim
)
par(new=T)
#dev.off()
pf <- pf + 1
} # end bbb loop
pc <- pc + 1
title(paste0(out_name,"\nbeta=",min(beta_pars)," to ",max(beta_pars),", death = ",death_access,"\n",ttl_list[pc-1]))
text(years-40,ylim[2],paste0("b=",beta_pars[1]),col=colvec[1]) # low beta
text(years-10,ylim[2],paste0("b=",max(beta_pars)),col=colvec[length(colvec)]) # high beta
} # end wastes loop
out_name <- out_names[2]
# choose your death value you want to plot at the end
death_access <- 0.1
# pdf(paste0(getwd(),"/npsi_model_plot_",ttl,".pdf"),onefile=T,width=10,height=8,paper="a4r")
# then run this part to plot in your live R session
layout(matrix(c(1,2,3), 1, 3, byrow = TRUE)) # set plot window
wastes <- list(out_tibble_wh,out_tibble_wd,out_tibble_ws)
ttl_list <- c("Host waste","Drool waste","Summed waste")
pf <- 1 # plot feature counter
pc <- 1 # plot number counter
colvec<-brewer.pal(length(beta_pars),"RdYlGn") # !!! make sure length of colour palette is > beta_pars
# dp <- 0.2
# bp <- 0.2
for(out_tibble in wastes){
for(beta_access in beta_pars){
outplot <- filter(out_tibble, death == death_access & beta == beta_access)
outplot <- outplot$outs ; outplot <- as.data.frame(outplot) # clean output
outplot$"Total host population" <- outplot[,"S"] + outplot[,"I"] # add sum host population
colnames(outplot) <- out_names
if(out_name==names(outplot[2]) | out_name==names(outplot[3])){ylim=c(0,200)}else{ylim=c(0,100)};ylim # set proper ylim
if(pf==length(beta_pars)+1){pf <- 1} # reset pf to reset colour palette for new plot
plot(outplot[,1],outplot[,out_name],type="l",las=1,bty="n",
xlab="Time (years)",ylab=out_name,
col=colvec[pf],#lty=pf,
ylim=ylim
)
par(new=T)
#dev.off()
pf <- pf + 1
} # end bbb loop
pc <- pc + 1
title(paste0(out_name,"\nbeta=",min(beta_pars)," to ",max(beta_pars),", death = ",death_access,"\n",ttl_list[pc-1]))
text(years-40,ylim[2],paste0("b=",beta_pars[1]),col=colvec[1]) # low beta
text(years-10,ylim[2],paste0("b=",max(beta_pars)),col=colvec[length(colvec)]) # high beta
} # end wastes loop
packages <- c("RCurl","RColorBrewer","viridis","deSolve","ggplot2","ggthemes","dplyr","tibble","purrr","reshape2","tidyr","zoo","plyr")
ppp <- lapply(packages,require,character.only=T)
if(any(ppp==F)){cat("\n\n\n ---> Check packages are loaded <--- \n\n\n")}
if(any(ppp==F)){cat("\n\n\n ---> Check packages are loaded <--- \n\n\n")}
# TO DO
# need to balance out nutrient and plant biomass equations (minus terms)
# use parasite-host worm model instead (keep biomass)
# remove waste parts from susc and infected hosts b/c waste doesnt detract from host biomass.
## STATE VARIABLES (units = biomass)
## N = nutrients
## P = plants
## S = susceptible ungulate hosts
## I = infected ungulate hosts
## PARAMETERS
## r = intrinsic growth rate of plants
## K = carrying capacity of plants
## a = rate of nutrient addition
## l = nutrient loss rate
## fp = rate of plant nutrient uptake
## beta = transmission rate
## es = assimilation efficiency of susceptible hosts
## ei = assimilation efficiency of infected hosts
## fs = feeding rate of susceptible hosts
## fi = feeding rate of infected hosts
## d = background host death rate
## v = mortality rate from infection
## ws = rate of waste production from susceptible hosts
## wi = rate of waste production from infected hosts
##########################################################################################
################################# User inputs for model ##################################
##########################################################################################
# set your working directory
# E.g "/Users/malishev/dope_models/my_dope_model/"
wd <- "paste the path to where you saved the model here (with these quotes)"
setwd(wd)
# set parameter ranges (min 0, max 1)
beta_access <- 0.1 # choose your beta value you want to plot at the end
death_access <- 0.9 # choose your death value you want to plot at the end
colvv <- "orange" # choose your plot line colour
## initial conditions
N <- 200 # size of nutrient biomass in env
P <- 200 # initial products in env
S <- 20 # num of susceptible hosts
In <- 20 # num of infected hosts
years <- 100 # number of years to run simulation
time.out <- 0.01 # simulation time step (0.01 = 1 year if years = 100)
##########################################################################################
##################################### Setup simulation model #############################
##########################################################################################
# ---------------------- run the model from here # ----------------------
# load packages
packages <- c("RCurl","RColorBrewer","viridis","deSolve","ggplot2","ggthemes","dplyr","tibble","purrr","reshape2","tidyr","zoo","plyr")
if (require(packages)) {
install.packages(packages,dependencies = T)
require(packages)
}
ppp <- lapply(packages,require,character.only=T)
if(any(ppp==F)){cat("\n\n\n ---> Check packages are loaded <--- \n\n\n")}
# pull plotting function
script <- getURL("https://raw.githubusercontent.com/darwinanddavis/plot_it/master/plot_it.R", ssl.verifypeer = FALSE)
eval(parse(text = script))
display.brewer.all()
# Set global plotting parameters
cat("plot_it( \n0 for presentation, 1 for manuscript, \nset colour for background, \nset colour palette 1. use 'display.brewer.all()', \nset colour palette 2. use 'display.brewer.all()', \nset alpha for colour transperancy, \nset font style \n)")
plot_it(0,"blue","Blues","YlOrRd",1,"mono") # set plot function params
plot_it_gg("white") # same as above for ggplot
# set param space
beta_pars <- seq(0.1,1,0.1) # transmission rate in model
death_pars <- seq(0.1,1,0.1) # death rate in model
# desired outputs
out <- list()
out_master <- list() # NPSI output
out_tibble <- tibble()
outplot <- list()
param_space <- list(beta_pars,death_pars) # summed parameter space
# create empty list
out_master <- rep(
list(structure(list(
pars = numeric(),
outs = list()
),
.Names = c("Parameter", "Output")))
,prod(as.numeric(summary(param_space)[,1]))
)
sc <- 1 # timer in simulation model
##########################################################################################
################################### create simulation model  #################################
# to set pars as individual beta and death values
npsi_func <- function(){ # start npsi_func
# ------- start simulation # -------
for(beta in beta_pars){ # pass through beta values
for(death in death_pars){ # pass through death values
parameters<-c(r=0.2, K=100, a=500, l=5, fp=0.5, beta=beta,
es=0.1, ei=0.05, fs=0.2, fi=0.1,
d=death, v=0.1, ws=0.05, wi=0.09)
state<-c(N=N, P=P, S=S, I=In) # set initial conditions
NPSI<-function(t, state, parameters) {
with(as.list(c(state, parameters)),{
dN.dt <- a - l*N - fp*N*P + (d+ws)*S + (d+v+wi)*I  # nutrients in env
dP.dt <- fp*N*r*P*(1-(P/K)) - P*(fs*S+fi*I) # plants produced
dS.dt <- P*(es*fs*S + ei*fi*I) - beta*S - (d+ws)*S # susceptible host population change
dI.dt <- beta*S - (d+v+wi)*I # infected hosts population change
list(c(dN.dt, dP.dt, dS.dt, dI.dt)) # compile outputs
})
} # end npsi function
# -------  global output # -------
times <- seq(0, years, by=time.out) # set time horizon for simulation (years)
out <- ode(y=state, times=times, func=NPSI, parms=parameters) # run simulation model
out <- data.frame(out)
# save outputs
# out_master[[length(out_master) + 1]] <- out # working with out_master <- list()
out_master[[sc]]$Output <- out # save output for each run
out_master[[sc]]$Parameter[1] <- beta # save beta for each run
out_master[[sc]]$Parameter[2] <- death # save death for each run
sc <- sc + 1
} # end death pars
} # end beta pars
# -------  clean output # -------
# save simulation model to global vector (tibble)
out_tibble <- tibble(
params = map(out_master, "Parameter"),
outs = map(out_master, "Output")
) %>%
mutate(
beta = map(params, 1),
death = map(params, 2)
) %>%
select(beta, death, outs)
# ------- plotting ----------
# start save plot to local dir
pdf(paste0(getwd(),"/npsi_model_plot.pdf"),onefile=T,width=10,height=8,paper="a4r")
outplot <- filter(out_tibble, death == death_access & beta == beta_access)
outplot <- outplot$outs ; outplot <- as.data.frame(outplot) # clean output
outplot$"Total host population" <- outplot[,"S"] + outplot[,"I"] # add sum host population
# plot results
layout(matrix(c(1,2,3,4,5,5), 2, 3, byrow = TRUE)) # set plot window
colnames(outplot) <- c("Time", #1
"Nutrient biomass", #2
"Plant biomass", #3
"Hosts (susceptible)", #4
"Hosts (infected)", #5
"Total hosts") #6
for (name in names(outplot)[c(3:5,2,6)]){ # start plot for product, suscep hosts, infec hosts, nutrients, and total hosts
plot(outplot[,1],outplot[,name],type="l",las=1,bty="n",
xlab="Time (years)",ylab=name,col=colvv,
ylim=c(0,round_any(max(outplot[,name]),10,ceiling))
)
title(paste0(name,"\nbeta = ",beta_access," , death = ",death_access))
} # end plot
# add mean plot
dev.off() # save output to dir
cat(paste0("\n\n\nPlot is saved in \n",getwd(), "\nas npsi_model_plot.pdf\n\n\n"))
return(out_tibble)
} # ------- end npsi_func
### run model function
out_tibble <- npsi_func()
################################### end simulation model  #################################
##########################################################################################
################################### plot results manually  #################################
# set parameter ranges (min 0, max 1)
beta_access <- 0.9 # choose your beta value you want to plot at the end
death_access <- 0.1 # choose your death value you want to plot at the end
colvv <- "orange" # choose your plot line colour
# then run this part to plot in your live R session
outplot <- filter(out_tibble, death == death_access & beta == beta_access)
outplot <- outplot$outs ; outplot <- as.data.frame(outplot) # clean output
outplot$"Total host population" <- outplot[,"S"] + outplot[,"I"] # add sum host population
layout(matrix(c(1,2,3,4,5,5), 2, 3, byrow = TRUE)) # set plot window
colnames(outplot) <- c("Time", #1
"Nutrient biomass", #2
"Plant biomass", #3
"Hosts (susceptible)", #4
"Hosts (infected)", #5
"Total hosts") #6
for (name in names(outplot)[c(3:5,2,6)]){ # start plot for product, suscep hosts, infec hosts, nutrients, and total hosts
plot(outplot[,1],outplot[,name],type="l",las=1,bty="n",
xlab="Time (years)",ylab=name,col=colvv,
ylim=c(0,round_any(max(outplot[,name]),10,ceiling))
)
title(paste0(name,"\nbeta = ",beta_access," , death = ",death_access))
} # end plot
install.packages("metagear")
packages <- c("meme","ggimage","ggplot2","grid")
install.packages(packages, dependencies = T)
install.packages(packages, dependencies = T)
lapply(packages,require,character.only=T)
snail <- "snail_black.png" # img1
snail
dc <- "dc.jpg" # img2
text_upper <- "I don't usually get schisto" # upper text for meme
text_lower <- "but when I do, \n it's hypothesis-driven" # lower text for meme
#################################### create meme
img <- paste0(wd,"/",dc) # set image path
mm <- meme(img,text_upper, text_lower, # create meme
size = 2, color = "white", bgcolor = "purple")
wd <- "/Users/malishev/Documents/Melbourne Uni/Programs/R code/meme" # working dir
snail <- "snail_black.png" # img1
dc <- "dc.jpg" # img2
text_upper <- "I don't usually get schisto" # upper text for meme
text_lower <- "but when I do, \n it's hypothesis-driven" # lower text for meme
#################################### create meme
img <- paste0(wd,"/",dc) # set image path
mm <- meme(img,text_upper, text_lower, # create meme
size = 2, color = "white", bgcolor = "purple")
mm
out_dir <- tempfile("dc_meme",wd,fileext=".png") # save to dir
img <- paste0(wd,"/",snail)
mm <- meme(img)
df <- data.frame(x = sample(20), y = rnorm(20)) # generate some data
height <- 1 # height of img
width <- 1.5 # width of img
# plot
ggplot(df, aes(x, y)) + # plot
geom_line() + # join points with a line
geom_subview(aes(x, y), data=df, subview=mm, width=width, height=height) + # plot img vectors
theme_classic() + # use minimal plot design
ggtitle("Snail (or any image) as data points")
vv <- out_tibble[3][1]$outs[[11]]$N
spectrum(vv,col="white")
periodogram(vv)
spectrum(vv,col="white")
library(deSolve)
Worm_stoich<-function (t, y, parameters){
N = y[1]; A = y[2]; QA = y[3]; H = y[4]; QH = y[5]; P = y[6]
with(
as.list(parameters),
{
ingestion = a_max*(P/H)/(a_h + a_max*P/H)
dNdt = s - l*N - v*N/(hN + N)*(QAx - QA)/(QAx - QAn)*A +
dA*A*QA + dH*H*QH + dP*P*qP + # a*P*QH + (dP+dH+a)*P*qP + a*P^2/H*(k+1)/k*qP +
f*A*H*QA*(1 -ingestion)*(1 - (QHx - QH)/(QHx - QHn)) +
qP*(1 - sigma*f*H/(dE + f*H))*f*A*H*ingestion*min(1,QA/qP)
dAdt = mu_A*(1-A/K)*(1 - QAn/QA)*A - dA*A - f*A*H
dQAdt = v*(N/(hN + N))*(QAx - QA)/(QAx - QAn) - mu_A*(1-A/K)*(1 - QAn/QA)*QA
dHdt = f*A*H*(1 - ingestion)*(1 - QHn/QH) - dH*H
dQHdt = f*A*QA*(1 -ingestion)*(QHx - QH)/(QHx - QHn) - f*A*(1 - ingestion)*(1 - QHn/QH)*QH
dPdt = (sigma*f*H/(dE + f*H))*f*A*H*ingestion*min(1,QA/qP) - dP*P #- a*P^2/H*(k+1)/k
res = c(dNdt,dAdt,dQAdt,dHdt,dQHdt,dPdt)
list(res)
}
)
}
params = c(# nutrient pars
s = 0, l = 0,
# autotroph pars
v = 1, mu_A=1, K=100, hN = 1, QAx = 0.05, QAn = 0.01, dA = 0.01,
# Host pars
f =0.1, QHx = 0.3, QHn = 0.05, dH = 0.05,
# Parasite pars
a_max = 5, a_h = 1, qP = 0.2, dP = 0.1, dE = 0.1, sigma=0.9)
inits = c(N = 10, A = 1, QA = 0.02, H = 1, QH = 0.1, P = 1)
run = lsoda(y = inits, times=0:1000, parms = params, func=Worm_stoich)
plot(run)
Total_N = run[,"N"] + run[,"A"]*run[,"QA"] + run[,"H"]*run[,"QH"] + run[,"P"]*params["qP"]
var(Total_N)
plot(run[,"time"], Total_N, typ="l")
N.out = numeric()
A.out = numeric()
QA.out = numeric()
H.out = numeric()
QH.out = numeric()
P.out = numeric()
for (n in 1:20){
inits["N"] = n
output = lsoda(y = inits, times=0:5000, parms = params, func=Worm_stoich)
N.out[n] = output[5001, 2]
A.out[n] = output[5001, 3]
QA.out[n] = output[5001, 4]
H.out[n] = output[5001, 5]
QH.out[n] = output[5001, 6]
P.out[n] = output[5001, 7]
}
par(mfrow=c(2, 3))
plot(1:20, N.out, typ="l")
plot(1:20, A.out, typ="l")
plot(1:20, QA.out, typ="l")
plot(1:20, H.out, typ="l")
plot(1:20, QH.out, typ="l")
plot(1:20, P.out, typ="l")
plot(run)
plot(run,col='white')
params = c(# nutrient pars
s = 0, l = 0,
# autotroph pars
v = 1, mu_A=1, K=100, hN = 1, QAx = 0.05, QAn = 0.01, dA = 0.01,
# Host pars
f =0.1, QHx = 0.05, QHn = 0.02, dH = 0.05,
# Parasite pars
a_max = 5, a_h = 1, qP = 0.2, dP = 0.1, dE = 0.1, sigma=0.9)
inits = c(N = 10, A = 1, QA = 0.02, H = 1, QH = 0.1, P = 1)
run = lsoda(y = inits, times=0:1000, parms = params, func=Worm_stoich)
plot(run,col='white')
params = c(# nutrient pars
s = 0, l = 0,
# autotroph pars
v = 1, mu_A=1, K=100, hN = 1, QAx = 0.05, QAn = 0.01, dA = 0.01,
# Host pars
f =0.1, QHx = 0.05, QHn = 0.02, dH = 0.05,
# Parasite pars
a_max = 5, a_h = 1, qP = 0.2, dP = 0.1, dE = 0.1, sigma=0.9)
inits = c(N = 10, A = 1, QA = 0.02, H = 1, QH = 0.04, P = 1)
run = lsoda(y = inits, times=0:1000, parms = params, func=Worm_stoich)
plot(run,col='white')