material from katrine

lecturers/katrine/POPweb.R

+
+
+# Visual description of parameters relevant for distribution
+# persistent organic pollutants (POPs) in an aquatic foodweb 
+#
+# 	Foodweb linkages	(T)
+#	Lipid content	(L)
+#	Sediment contact	(S)
+#
+# A circular diagram of foodweb linkages with arrows weighted in relation to diet importance (T)
+# Individual foodweb compartments with border thickness proportional to lipid content (L)
+# Each compartment symbol divided vertically according to fraction sediment contact (S)
+# Foodweb compartment names are are taken from the column names of T
+#
+# Tom Andersen, 2000; Converted to R 2014
+
+library(shape)
+
+draw.flows <- function(T) {
+
+	N <- dim(T)[1]
+
+	# Scale to max. flow
+	T <- T / max(as.numeric(T))
+
+	# Place compartment nodes on a N-sided regular polygon
+	t <- (1:N) * (2*pi/N)
+	x <- cos(t)
+	y <- sin(t)
+
+	R <- sin(pi / N) / 2
+
+	# Find non-zero internal flows, calculate connection lines
+	k <- which(T > 0, arr.ind=TRUE)	
+	w <- T[k] 
+
+	X <- cbind(x[k[, 1]], x[k[, 2]]) 
+	Y <- cbind(y[k[, 1]], y[k[, 2]]) 
+
+	a <- atan2(diff(t(X)), diff(t(Y))) 
+
+	# Draw internal flow diagram lines, scale line width to flow magnitude
+	for (i in (1 : length(a))) {
+   		ax <- cbind(
+			X[i,2] + R * (sin(a[i] - pi)),
+         		X[i,2] + R * (sin(a[i] - pi) +  w[i]     * sin(a[i] - (7/6) * pi)),
+         		X[i,2] + R * (sin(a[i] - pi) + (w[i]/2)  * sin(a[i] - (9/8) * pi)),
+         		X[i,1] + R * (w[i]/2) * sin(a[i] + (1/8) * pi),
+         		X[i,1] + R * (w[i]/2) * sin(a[i] - (1/8) * pi),
+         		X[i,2] + R * (sin(a[i] - pi) + (w[i]/2)  * sin(a[i] - (7/8) * pi)),
+         		X[i,2] + R * (sin(a[i] - pi) +  w[i]     * sin(a[i] - (5/6) * pi)), 
+         		X[i,2] + R * (sin(a[i] - pi)))
+   
+   		ay <- cbind(
+			Y[i,2] + R * (cos(a[i] - pi)),
+         		Y[i,2] + R * (cos(a[i] - pi) +  w[i]     * cos(a[i] - (7/6) * pi)),
+         		Y[i,2] + R * (cos(a[i] - pi) + (w[i]/2)  * cos(a[i] - (9/8) * pi)),
+         		Y[i,1] + R * (w[i]/2) * cos(a[i] + (1/8) * pi),
+         		Y[i,1] + R * (w[i]/2) * cos(a[i] - (1/8) * pi),
+         		Y[i,2] + R * (cos(a[i] - pi) + (w[i]/2)  * cos(a[i] - (7/8) * pi)),
+         		Y[i,2] + R * (cos(a[i] - pi) +  w[i]     * cos(a[i] - (5/6) * pi)),
+         		Y[i,2] + R * (cos(a[i] - pi)))
+   
+   		polygon(ax, ay, col="lightgray", border="darkgray")
+	}
+}
+
+
+draw.compartments <- function(L, S) {
+
+	N <- length(S)
+	R <- sin(pi / N) / 2
+   	r <- (1 - 2 * L) * R	
+	t <- (1:N) * (2*pi/N)
+
+	x <- cos(t)
+	y <- sin(t)
+
+	# Place compartment nodes on a N-sided regular polygon
+	# Make annulus thickness proportional to lipid content
+	# Color node according to sediment contact
+	# Label nodes with compartment names
+
+
+	for (i in (1:N)) {
+
+		plotcircle(mid = c(x[i], y[i]), r = R, lwd = 1, lcol = rgb(0.5,0.5,0.5,0.5),
+			col = rgb(t(col2rgb("slategray3")) / 255, alpha=(1 - S[i])))
+		plotcircle(mid = c(x[i], y[i]), r = R, lwd = 1, lcol = rgb(0.5,0.5,0.5,0.5),
+			col = rgb(t(col2rgb("tan3")) / 255, alpha=S[i]))
+
+		plotcircle(mid = c(x[i], y[i]), r = r[i], lwd = 1, lcol = rgb(0.5,0.5,0.5,0.5),
+			col = rgb(t(col2rgb("slategray2")) / 255, alpha=(1 - S[i])))
+		plotcircle(mid = c(x[i], y[i]), r = r[i], lwd = 1, lcol = rgb(0.5,0.5,0.5,0.5),
+			col = rgb(t(col2rgb("tan2")) / 255, alpha=S[i]))
+   	}
+  
+	text(x, y, rownames(L), cex=0.7)
+}
+
+popweb <- function(WEB) {
+	
+	s <- 1.45 # Scaling factor
+	plot(s * c(-1, 1), s * c(-1, 1), type="n", xlab="", ylab="", axes=FALSE, asp=1)
+
+	draw.flows(t(WEB$Pfw))
+	draw.compartments(WEB$Lip, WEB$fSed)
+}
+
+
+

lecturers/katrine/POS.R

+# Power Of Size (POS) model of hydrophobic contaminats in foodwebs 
+#
+# based on: Hendriks et al (2001)
+#   "The power of size. 1. Rate constants and equilibrium ratios for
+#    accumulation of organic substances related to octanol-water
+#    partition ratio and species weight" Env. Tox. Chem. 20(7): 1399-1420
+#
+# Tom Andersen, NIVA - 2003; Converted to R by TA, 2014
+
+# All functions have the same inputs:
+#	POP	- Hydrophobic organic substance properties
+#	WEB	- Foodweb properties
+#	par	- Model parameters (default stdPOSpar)
+
+# POP: Hydrophobic organic substance property object
+#	KOW - Octanol -water partition coeff. (-)
+#	Cwd - Dissolved conc. in water (ng/L)
+#	Csd - Dissolved conc. in sediment (ng/L)
+
+# WEB: Foodweb property object (n compartments)
+#	Vol	(n x 1) Body volume (L = kg w.w.)
+#	Lip	(n x 1) Lipid volume fraction (-)
+#	Kmet	(n x 1) Metabolic degradation rate constant (1/d)
+#	fsed	(n x 1) Fraction time spent in sediments (-)
+#	Pfw	(n x n) Foodweb fractional diet matrix (-)
+
+# Standard parameters for the Hendriks et al (2001) model:
+stdPOSpar <- c(
+  Ksd = 0.25,	# ( - ) - Size dependency exponent 
+  WE0 = 200,	# (1/d) - Water absorbtion/excretion rate at reference size (1 kg) 
+  FI0 = 0.005,	# (1/d) - Food ingestion rate at reference size (1 kg) 
+  GR0 = 0.0006,	# (1/d) - Growth rate at reference size (1 kg) 
+  EFF = 0.75,	# ( - ) - Assimilation efficiency 
+  WRS = 2.8E-3,	# ( d ) - Water layer diffusion resistance, surface (1.4 - 4.1)E-3
+  WRG = 1.1E-5,	# ( d ) - Water layer diffusion resistance, gut (0.0 - 3.9)E-5 
+  LPR = 68)		# ( d ) - Lipid layer permeation resistance (30 - 110) 
+
+POSrates <- function(POP, WEB, par=stdPOSpar) {
+# Rate constant vectors and matrices for given POP and WEB 
+#
+#	Kin   - Water uptake rate constant vector (L/kg/d = 1/d)
+#	Kout  - Water release rate constant vectors (1/d)
+#	Kfeed - Diet uptake rate constant matrix (1/d)
+
+	with(as.list(c(POP, WEB, par)), {
+
+		Psize <- Vol^(-Ksd)	# Size dependency of rate processes (eq. 1)
+
+		fLip  <- Pfw %*% Lip	# Average lipid content of food
+		feed  <- which(fLip > 0)# Feeding compartment (fLip == 0 -> non-feeding)
+
+		# Rate coefficients for exchange with environment (water/sediment)
+		KWU <- Psize / (WRS + (LPR / KOW) + (1 / WE0))	# Uptake from water (eq. 3 + 5)	
+		KWL <- KWU / (((KOW - 1) * Lip) + 1)		# Loss to water (eq. 4)
+
+		# Compute last term of eq. 6 & 9 (avoid division by zero for non-feeding compartments)		
+		FF <- 0.0 * Lip   # Allocate vector of zeros, only feeding compartments will be changed
+		FF[feed] <- Psize[feed] / (WRG + (LPR / KOW) + (1 / (fLip[feed] * KOW * (1 - EFF) * FI0)))
+
+		# Rate coefficients for uptake from food
+		KFL <- (1 / (Lip * (KOW - 1) + 1)) * FF         # Loss by egestion (eq. 9)
+		KFU <- (EFF / (1 - EFF)) * KFL      		# Uptake from food (eq. 6)
+		KGD <- GR0 * Psize					# Growth dilution	(eq. 10)
+
+		# Combined uptake, loss, and feeding rates
+		return(list(
+			Kin   = KWU,					# Uptake from environment 
+			Kout  = cbind(KWL, KGD, KFL),			# Losses to environment (excluding metabolic loss)
+			Kfeed = sweep(Pfw, 1, KFU, FUN="*")))	# Uptake from foodweb allocated by fractional diet
+	})
+}
+
+POSmatrices <- function(POP, WEB, par=stdPOSpar) {
+# Linear DE model matrices for given POP and WEB
+#
+#	dx / dt = A x + B z
+#
+# where x is foodweb POP concentrations 
+#   and z is environmental POP concentrations
+
+	m <- POSrates(POP, WEB, par)
+	Ksum <- rowSums(m$Kout)
+
+	return(list(
+		A=m$Kfeed - diag(as.numeric(Ksum + WEB$Kmet)), 
+		B=diag(as.numeric(m$Kin))))	
+}
+
+POSequ <- function(POP, WEB, par=stdPOSpar) {
+# Equilibrium foodweb concentrations for given POP and WEB
+
+	# Dissolved environmental concentration for a given compartment
+	# is proportional to fraction of time spent in sediments
+	Cenv <- (1 - WEB$fSed) * POP["Cwd"] + WEB$fSed * POP["Csd"]
+
+	m <- POSmatrices(POP, WEB, par)
+
+	# Equilibrium concentration (ng / kg)
+	return(with(m, solve(-A, B %*% Cenv)))
+}
+
+POSflow <- function(POP, WEB, par=stdPOSpar) {
+# Equilibrium flows for given POP and WEB
+
+	m <- POSrates(POP, WEB, par)
+
+	# Equilibrium concentrations
+	Cequ <- POSequ(POP, WEB, par)
+
+	# Equilibrium flow components (ng/kg/d)
+	FWU <- m$Kin * (1 - WEB$fSed) * POP["Cwd"]	# Uptake from open water
+	FSU <- m$Kin * WEB$fSed * POP["Csd"]		# Uptake from pore water
+	FFU <- m$Kfeed %*% Cequ					# Uptake from food
+	FWL <- m$Kout[, 1] * Cequ				# Loss to water
+	FGL <- m$Kout[, 2] * Cequ				# Loss by growth dilution
+	FEL <- m$Kout[, 3] * Cequ				# Loss by egestion
+	FML <- WEB$Kmet * Cequ					# Loss by metabolic degradation
+
+	F <- cbind(FWU, FSU, FFU, FWL, FGL, FEL, FML)
+	colnames(F) <- c("FWU", "FSU", "FFU", "FWL", "FGL", "FEL", "FML")
+
+	return(F)
+}
+
+library(deSolve)
+
+linsys <- function(Time, X, P) {
+
+	dX.dt <- P$A %*% X + P$b
+
+	return(list(dX.dt))
+}
+
+POSdyn <- function(C0, t, POP, WEB, par=stdPOSpar) {
+# Solve dynamical linear DE model for given POP and WEB 
+# from initial concentrations C0 over time interval t
+
+	# Dissolved environmental concentration for a given compartment
+	# is proportional to fraction of time spent in sediments
+	Cenv <- (1 - WEB$fSed) * POP["Cwd"] + WEB$fSed * POP["Csd"]
+
+	m <- POSmatrices(POP, WEB, par)
+	params <- list(A = m$A, b = m$B %*% Cenv)
+
+	x <- ode(C0, t, linsys, params)
+
+	return(as.data.frame(x))
+}
+

lecturers/katrine/POSGobas.1993.R

+# Example parameter set on PCBs in Lake Ontario, based on Gobas (1993)
+#  "A model for predicting the bioaccumulation of hydrophobic organic 
+#   chemicals in aquatic food-webs: application to Lake Ontario" 
+#   Ecol. Modelling 69: 1-17
+
+KOW <- 10^6.6	# Weighted average KOW for total PCBs
+
+Cwt <- 1.1		# Total PCB conc. in water (ng/l)
+Cst <- 570000	# Total PCB conc. in sediment (ng/g DW)
+
+OMw <- 2.5E-7	# Organic matter content of water    (kg/l <-> 0.25 mg/l)
+OMs <- 2.0E-2	# Organic matter content of sediment (kg/l <-> 2#)
+
+dOM <- 0.9		# Density of organic matter (kg/l)
+
+fOMw <- OMw / dOM	# Volume fraction OM in water
+fOMs <- OMs / dOM	# Volume fraction OM in sediment (??? what about water content ???)
+
+PCB <- c(
+	KOW = KOW,
+	Cwd = (1 / (1 + KOW * fOMw)) * Cwt,	# Dissolved PBC concentration in water
+	Csd = (1 / (1 + KOW * fOMs)) * Cst)	# Dissolved PBC concentration in sediment
+
+
+#                Mysis    Pontop.  Oligoch. Cottus   Alosa    Osmerus  Salmon.
+Vol  <- matrix(c(0.00001, 0.0005,  0.0001,  0.0054,  0.032,   0.016,   2.41), ncol=1) # Body volumes (1 L == 1 kg WW)
+Lip  <- matrix(c(0.05,    0.03,    0.01,    0.08,    0.07,    0.04,    0.16), ncol=1) # Lipid fractions		
+Kmet <- matrix(c(0,       0,       0,       0,       0,       0,       0   ), ncol=1) # Metabolic degradation rates
+fSed <- matrix(c(0,       1,       1,       0,       0,       0,       0   ), ncol=1) # Fraction of time in sediments
+
+# Food preferences, P(i,j) is preference of i for j
+		
+Pfw <- matrix(c(
+   0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,  # Mysis
+   0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,  # Pontoporeia
+   0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,  # Oligochaeta
+   0.18, 0.82, 0.00, 0.00, 0.00, 0.00, 0.00,  # Cottus
+   0.60, 0.40, 0.00, 0.00, 0.00, 0.00, 0.00,  # Alosa
+   0.54, 0.21, 0.00, 0.25, 0.00, 0.00, 0.00,  # Osmerus
+   0.00, 0.00, 0.00, 0.10, 0.50, 0.40, 0.00), # Salmonidae
+	nrow=7, ncol=7, byrow=TRUE)
+
+biota <- c("Mysis", "Pontoporeia", "Oligochaeta", "Cottus", "Alosa", "Osmerus", "Salmonidae")
+
+rownames(Vol) <- rownames(Lip) <- rownames(Kmet) <- rownames(fSed) <- biota
+rownames(Pfw) <- colnames(Pfw) <- biota
+
+lake.Ontario <- list(Vol=Vol, Lip=Lip, Kmet=Kmet, fSed=fSed, Pfw=Pfw)
+
+source("POS.R") 		# Load POS model functions
+source("popweb.R")	# Load foodweb visualisation function
+
+popweb(lake.Ontario)	# Food web structure visualization
+
+(r <- POSrates(PCB, lake.Ontario)) # Rate constants
+
+Kin.w <- r[["Kin"]][,1]
+Kout.w <- r[["Kout"]][,1]
+BCF <- Kin.w / Kout.w	# Bioconcentration factors
+
+par(mfrow=c(1,3))
+dotchart(Kin.w, main="Kin")
+dotchart(Kout.w, main="Kout")
+dotchart(BCF, main="BCF")
+
+(A <- POSmatrices(PCB, lake.Ontario)[["A"]])
+lambda <- eigen(A)$values
+
+par(mfrow=c(2,1))
+plot(lambda, xlab="", ylab="Eigenvalue (1/d)")   #similar to resiliense of the system: is the system.  changing fast or slow. smallest absolute value means ow resilience and long response time.
+
+plot((-1 / lambda) / 365, xlab="", ylab="Response time (year)")
+
+Cinf <- POSequ(PCB, lake.Ontario)[,1] / 1000
+dotchart(Cinf,xlab="Equilibrium concentration (g/L)")
+
+
+C0 <- 0 * Vol # Start all compartments with zero PCB
+x <- POSdyn(C0, t=0:2000, PCB, lake.Ontario)
+
+t <- x[, 1] / 365 	# Days -> Years
+C <- x[, -1] / 1000	# ng/kg -> g/kg
+C1 <- C[length(t), ]	# Final concentration
+
+matplot(t, C, type="l", lty=1, lwd=3, col=rainbow(7),
+	xlab="Time (years)", ylab="PCB (g/kg)", xlim=c(0,6))
+text(5.5, C1, biota, pos=4, cex=0.7, col=rainbow(7)) 
+
+
+# Make a new PCB-free environment...
+no.PCB <- PCB
+no.PCB["Cwd"] <- 0.0
+no.PCB["Csd"] <- 0.0
+
+# Start all compartments at equilibrium in a zero PCB environment
+x <- POSdyn(Cinf, t=0:2000, no.PCB, lake.Ontario)
+
+t <- x[, 1] / 365 	# Days -> Years
+C <- x[, -1] / 1000	# ng/kg -> g/kg
+C1 <- C[length(t), ]	# Final concentration
+
+matplot(t, C, type="l", lty=1, lwd=3, col=rainbow(7),
+	xlab="Time (years)", ylab="PCB (g/kg)", xlim=c(0,6))
+text(5.5, C1, biota, pos=4, cex=0.7, col=rainbow(7))
+