Simple Ecological modeling
1: Marine Ecosystem modeling
Ecosystem modeling is the mathematical representation of ecological systems from individuals to the whole biome which is studied to understand the real nature, the models can be analytic or simulations based. It has applications in a wide variety of disciplines. Today, I am going to show one analytical model and one simulation-based model of the ecosystem simply and efficiently.
Consider a conceptual simple marine food web model. The model consists of 6 state variables such as phytoplankton, zooplankton, detritus, ammonium, fish, and bottom detritus. They are connected by 13 flows. There is a forcing agent (sun) and one output variable (chlorophyll). Note: [ I am not considering the spatial scale. I will use the flows that explain the exchange of matter between the state variables].
So, we need an equation that describes conceptual mechanisms. We build our model based on the energy exchange between state variables
For instance, the ammonium ion concentration will increase due to mineralization and excretion of zooplankton. And it also decreases as a result of uptaking by phytoplankton.
Similarly, we can assume that the fish population will be increased by eating more zooplankton and will be reduced due to predation and natural death (it will become detritus later). We are expressing the exchange per day, for example, the ammonium will be in mmol Nm^-3 /day, so all flows are also denoted by the same unit.
So the whole food web can be represented by the following equations. And the law of conservation of mass will be applied to these equations.
(d detritus)/dt=f8+f13+f6+f12-f10-f7
(d phytoplankton)/dt=f1-f2-f8-f9
(d zooplankton)/dt=f2+f3+f8+f12+f13-f7-f10
(d fish)/dt=f5-f6-f12
(d ammonium)/dt=f7+f10+f4-f1
(d bottom_Detritus)/dt=f7+f9-f11 [here= +f denotes increasing rate and -f denotes decreasing rate of mass.]
We are considering that our equations will follow the conservation of mass because total mass should be constant, so
We will use the diagram library for drawing the food web. It is very light and useful for simple graphs based on a transition matrix, utilities to plot flow diagrams.
Code:
library(diagram)
opar <- par()
par(ask=TRUE)
par(mfrow=c(1,1))
figNr <- 1
subtitle <- function()
{
mtext(side=1,outer=TRUE,”Marine foodweb modeling by Hafez Ahmad “,cex=0.7,adj=1,line=-1.5)
mtext(side=1,outer=TRUE,paste(“ Figure:”,figNr,sep=””),cex=0.7,adj=0,line=-1.5)
figNr <<- figNr +1
}par(mfrow=c(1,1))
par(mar=c(0,0,0,0))
names <- c(“Phytoplankton”,expression(NH[4]^”+”),”Zooplankton”,”Detritus”,”Bottom_detritus”,”Fishes”)
M <- matrix(nrow=6,ncol=6,byrow=TRUE,data=c(
# p n z d b f
0,1,0, 0, 0, 0, #p
0,0,4, 10,11,0, #n
2,0,0, 0, 0, 0, #z
8,0,13,0, 0, 12,#d
9,0,0, 7, 0, 0, #b
0,0,5, 0, 0, 0 #f
))
pp<-plotmat(M,pos=c(1,2,1,2),curve=0,name=names,lwd=1,my=0.,cex.txt=0.8,
box.lwd=2,box.size=0.08,box.type=”square”,box.prop=0.5,
arr.type=”triangle”,arr.pos=0.4,shadow.size=0.01,prefix=”f”)# extra arrows: flow 5 to Detritus and flow 2 to detritus
phyto <-pp$comp[1,]
zoo <-pp$comp[3,]
nh3 <-pp$comp[2,]
detritus<-pp$comp[4,]
fish <-pp$comp[6,]# flow6->detritus
m2 <- 0.5*(zoo+fish)
m1 <- detritus
m1[1]<-m1[1]+ pp$radii[4,1]
mid<-straightarrow (from=m2,to=m1,arr.type=”triangle”,arr.pos=0.4,lwd=1)
text(mid[1],mid[2]+0.03,”f6",cex=0.8)# flow3->detritus
m2 <- 0.5*(zoo+phyto)
m1 <- detritus
m1[1] <-m1[1]+ pp$radii[3,1]*0.2
m1[2]<-m1[2] + pp$radii[3,2]
mid<-straightarrow (from=m2,to=m1,arr.type=”triangle”,arr.pos=0.3,lwd=1)
text(mid[1]-0.01,mid[2]+0.03,”f3",cex=0.8)# solar radiation
m1 <- 0.5*(nh3+phyto)
m2 <- c(0.25,0.9)
segments (m1[1],m1[2],m2[1],m2[2],lwd=1,lty=2)
text(m2[1]-0.01,m2[2]+0.03,”solar radiation”,adj=c(0.5,0.5))# chlorophyll
m1 <- phyto
m1[1] <-m1[1]+ pp$radii[1,1]
m2 <- m1
m2[1]<-m2[1]+0.25
segments (m1[1],m1[2],m2[1],m2[2],lwd=1)
textellipse(mid=m2,radx=pp$radii[1,1],rady=pp$radii[1,2],lwd=1,shadow.size=0,lab=”Chlorophyll”)
subtitle()
2: logistic population growth simulation
It is one of the most used models in fisheries science. Because it deals with the population, right? in the logistic model, A population’s per capita growth rate gets smaller and smaller as population size approaches a maximum imposed by limited resources in the environment. Populations of organisms don’t continue to grow exponentially. Instead, the density increases unit reaches a maximum sustainable density called carrying capacity of the environment and it is symbolized as k, was first described in 1845 by the Pierre-Francois and then again in 1920 by Raymond Pearl and Lowell Reed.
There are two versions of the model. Here, I will use the classical model. The logistic growth differential equation. In fact, it has been so influential in population modeling in ecology.
Here (d N)/dt is the rate of change,r=normal net growth rate,k=carrying capacity [the maximum population size that can be supported]
Here
We have set out carrying capacity, k=10, growth rate, r=1,
In the graph, the x-axis is population density with the carrying capacity (k) equal to 10 and the y-axis is the actual rate of population growth per individual. We will simulate 101 times at 0.2 intervals for the 4 initial conditions.
We will utilize desolve library. It has Functions that solve initial value problems of a system of first-order ordinary differential equations (‘ODE’), of partial differential equations (‘PDE’), of differential-algebraic equations (‘DAE’), and of delay differential equations. The package contains routines designed for solving ‘ODEs’ resulting from 1-D, 2-D, and 3-D partial differential equations (‘PDE’) that have been converted to ‘ODEs’ by numerical differencing.
Code:
library(deSolve)
#y<-rN(1-N/k)
r<-1;k<-10;init<-2
derivs<-function(t,N,params)
list(r*N*(1-N/k))times<-seq(0,20,0.2)
times
out<-ode(y=init,times=times,func=derivs,parms=NULL)
plot(out)
matplot(out)init<-12
out1<-ode(y=init,times=times,func=derivs,parms=NULL)init<-6
out3<-ode(y=init,times=times,func=derivs,parms=NULL)
init<-20
out4<-ode(y=init,times=times,func=derivs,parms=NULL)plot(out,out1,out3,out4,main=’logistic population growth simulation’,lwd=2,col=c(“#FF8C00”, “#1E90FF”, “#8B7D6B”, “#00EE00”))
legend(“bottomright”, legend=c(‘inital condition=2’,’int.=12',’ini.=5',’int.=20'),
col=c(“#FF8C00”, “#1E90FF”, “#8B7D6B”, “#00EE00”), bg=”transparent”,lty=1:2, cex=0.8)
I have attached the R [version 4.0.0] code in the GitHub repository. Please check it. Next post will be “how to make an interactive dashboard”
