This R Markdown is highlighting data used in my Biomathmatics class modeling wall calcium fluctuations in plant pollen tubes.
rm(list=ls(all=TRUE))
library(deSolve)
state <- c(WCA=1)
Unlike the majority of the work I have done in R, this model requires no input of data, and is completely contained within R. To achieve this, three main variables must be imputed into the script for the model to function: t, state, and parameters (these variables must be imported in a certain way for R to understand). The first variable imputed here is state, which lists an initial condition for R to begin the model at (this initial condition may be changed later).
For reference, the variables used below are listed here: WCA = free calcium that enters the cell [and is not taken up by the cell wall] FCA = free calcium found in extra cellular matrix RPD = rate of demethylized pectin delivery [to the cell wall] Alpha = the binding constant of calcium to pectin D = diffusion constant of calcium from the extra cellular matrix into the cell wall
CalciumFlux <- function(t, state, parameters){
with (as.list(c(t, state, parameters)),{
DWCA <- D*(FCA-WCA)-RPD*(FCA*Alpha/(1+FCA*Alpha))
list(c(DWCA))
})
}
The next series of commands imputed create the function. The function used above is based off of a predator-prey model that has been modified to simulate calcium flux. make sure to note the difference between brackets and parentheses during these steps!
times= seq(0, 100, by= .1)
Now we wish to give R the time parameters we wish to use for this function, which are imputed into R using the previous command.
for(P.d in (1) ){
RPD = 2
Alpha = .3
D = 1
parameters = c(Alpha=Alpha, RPD=RPD, D=D, FCA=100)
out=ode(y=state, times=times, func=CalciumFlux, parms=parameters)
}
This last step sets the other variables the function uses, as well as completing the ode command that R requires to interpret the function. The last thing to do is simply plot the function:
plot(out, xlab="time",ylab="wall calcium", col="green")
As you can see, this model shows that wall calcium plateaus quickly, just as a prey system model with a carrying capacity. To better understand how this system would oscillate however, we can rework the commands again:
This time, we change rate of pectin export (P.d) from one condition into ten:
for (P.d in (1:10) ){
alpha =.1
D=.05
parameters= c(alpha = alpha, P.d=P.d, D = D, Ca.ex = 100)
out = ode(y=state, times=times, func=Calcium, parms=parameters)
if (P.d ==1) plot(out[,1],out[,2],type="l",col=P.d,xlab="Time",
main=paste("D= ",D, " alpha = ",alpha," with different P.d"),ylab = "Ca.w")
if (P.d > 1) lines(out[,1],out[,2], col=P.d)
text(90,out[1000,2],P.d)
}
After plotting the function with different rates of pectin export, we see a different “carrying capacity” of wall calcium. This is very important, as is gives us an insight into how this model could be perfected to show an oscillation. Because this model only has one function, which is based off of the prey aspect of a Lotka-Volterra model, the incorporation of a second function would make this system oscillate. For example, if we added