# Analytical solution for diffusion of CO2 into the intercellular air space
mytime = seq(0, 10, 0.125) # this gives our time series with start, end and interval
CO2i_initial = 300 # set the initial intracellular CO2 to 300 (the units are parts per million, ppm)
CO2outside = 360 # set the constant atmoshperic CO2 to 360 ppm
conductance = 1 # set the conductance, determined by size, number and opening/closing of stomata
# integrated rate equation to evaluate intracellular CO2
CO2inside = CO2outside+(CO2i_initial-CO2outside)*exp(-1.56*conductance*mytime)
# Plot CO2i vs time as a line plot, and give the plot a title
plot(mytime, CO2inside, type="l", main="Diffusion into cell")
---------------------------------------------
# Euler method for numerical integration using an example of a simple exponential equation
# We need to do the following things:
# 1) Create the function that calculates the value of our dependent variable, Y, at each time point,
# based on the previous value and a formula for how the value is changing (slope or first
# differential).
# 2) Assign values for independent variables and initial values
# 3) Call our function (with parameters), and assign the result to a variable
# 4) Plot the result, and compare the numerical result with the analytical solution
# Y is set as a function that requires two parameters:
# "yinit" is the initial value of Y at time 0. Since our formula will only calculate a new value
# of Y based on the previous value, we need to provide an initial value
# "x" is our independent variable, time, which is a vector of values containing each of the time
# points at which we will calculate our dependent variable
# The function "y" does four things:
# 1) assigns the initial value to (the result) "y"
# 2) calculates how many times it needs to calculate a value for "y", and how big the time step is
# 3) loops through the calculation of "y" at each time point by first calculating the change in y,
# as "dy", and then adding the change to the previous value of y, "y[i-1]".
# 4) returns the resulting vector of y values
y = function(yinit,x){
y = c(yinit) # Initial value passed to function
# Divide x interval into n steps of length dx
# Doing this here, inside the function, will allow us to evaluate y easily over different
# intervals by passing different vectors of time points to the function
n = length(x) # calculate how many times "y" will need to be evaluated
dx=(x[n]-x[1])/(n-1) # calculate the size of the time step, this assumes that all time-steps
# are the same size
# Set-up a loop to calculate "y" at each time point. We don't need to calculate it at the first
# time point, t=0, because we provided an initial value, so we start at the second time point
for (i in 2:n) {
dy = f(x[i-1])*dx # the function "f" is defined below for the differential rate equation
y[i] = y[i-1]+dy # calculate "y" as HAVE = HAD + CHANGE
}
return(y)
}
f = function(x) {x^2} # This is the differntial rate equation (the slope, or first derivative,
# at time xi), that describes the rate of change in y given a value of x
# Create a vector of time points
xinit = 0 # start time
xmax = 2 # end time
dx = .01 # time-step
x = seq(xinit,xmax,dx) # vector of time points
yinit = 0 # initial value of the dependent variable
x1 = x # I don't know why this assignment was necessary
y1 = y(yinit,x) # call our function "y" with it's needed parameters, and assign the result to "y1"
plot(x1,y1, type="l") # Plot the result as a line plot, (again, I don't know why "x1" is necessary)
lines(spline(x, x^3/3, n=201), col="red") # Add a plot of (x cubed over 3) to the plot, which is the
# integrated version of this differential. Since our
# differential equation (x squared) didn't include any
# variable quantities, we can integrate it analytically
# and compare it to our numerically integrated result.
# Now try changing the time-step, "dx", and see how well
# the numerical and analytical solutions compare.
---------------------------------------------
# Euler method version of photosynthesis model - uses finite difference equations and quantity
# values at one time-point to estimate the change in each quantity that will have occured by the
# next time-point
# Constants
CO2out = 360; g = 1; kdiff = 1.56*g; PAR = 300
dt = .005; tmax = 10 # time-step interval, and ending time
t = seq(0,tmax,dt); n = tmax/dt+1 # create a vector of time points, and calculate the number of
# time points at which our quantities will be estimated
# Create a matrix which will contain all of our dependant variables at each time point
# row 1 = x[1,] =CO2i
# row 2 = x[2,] = PGA
# row 3 = x[3,] = G3P
# row 4 = x[4,] = sucrose
x = matrix(rep(0,4*n),nrow=4,ncol=n)
x[1,1] = 300; x[2,1] = 0; x[3,1] = 0 ;x[4,1] = 0 # Initial values
for(i in 2:n){ # don't need to calculate our intital values we just assigned so we loop up from 2
# evaluate each of the finite difference equations representing the change in our four
# dependant variables
dx1 = 1.56*g*(CO2out-x[1,i-1])*dt - ((100*(x[1,i-1]-40.4))/(x[1,i-1]+300*(1.7)))*dt
dx2 = ((100*(x[1,i-1]-40.4))/(x[1,i-1]+300*(1.7)))*dt - ((0.064*PAR*(x[1,i-1]-40.4))/(x[1,i-1]+80.8))*dt
dx3 = ((0.064*PAR*(x[1,i-1]-40.4))/(x[1,i-1]+80.8))*dt - min(50, x[3,i-1])*dt
dx4 = min(50, x[3,i-1])*dt # take the smaller value of 50 or the amount of G3P on the last step
# now determine the new value of each quantity as HAVE = HAD + CHANGE
x[1,i] = x[1,i-1]+dx1 # CO2in
x[2,i] = x[2,i-1]+dx2 # PGA
x[3,i] = x[3,i-1]+dx3 # G3P
x[4,i] = x[4,i-1]+dx4 # sucrose
}
plot(t, x[1,], type="l", lty=1, main="Photosynthesis Time Course", xlab="t/s", ylab="Conc/M",
ylim=c(0,375), bty="l")
lines(t, x[2,], lty=2)
lines(t, x[3,], lty=3)
lines(t, x[4,], lty=4)
text(8, 325, "CO2in")
text(8, 185, "PGA")
text(8, 25, "G3P")
text(8, 115, "Sucrose")
# flux equations from the Stella model
# diffusion - 1.56*g*(Ca-Ci)
# fixation - (100*(Ci-40.4))/(Ci+300*(1.7))
# reduction - (0.064*PAR*(Ci-40.4))/(Ci+80.8)
# sucrose formation - MIN(50,G3P)
---------------------------------------------
# Improved Euler Method, or Runge-Kutta 2 -
# Use the rate of change and quantity value at prior time-point to estimate the change in
# quantity, as in the Euler method above
# then calculate the change based on the estimated quantity at the end of the time-step
# use the average of the two estimated changes, at the start and end of the interval, as the
# change over the interval
# Constants
CO2out = 360; g = 1; kdiff = 1.56*g; PAR = 300
dt = .05; tmax = 10 # time-step interval, and ending time (larger steps)
t = seq(0,tmax,dt); n = tmax/dt+1 # create a vector of time points, and calculate the number of
# time points at which our quantities will be estimated
# Create a matrix which will contain all of our dependant variables at each time point
x = matrix(rep(0,4*n),nrow=4,ncol=n)
x[1,1] = 300; x[2,1] = 0; x[3,1] = 0 ;x[4,1] = 0 # Initial values
for(i in 2:n){
# Calculate the estimated change in each quantity using the finite difference equations with
# the values of the quantities at the prior time-point
dx1a = 1.56*g*(CO2out-x[1,i-1])*dt - ((100*(x[1,i-1]-40.4))/(x[1,i-1]+300*(1.7)))*dt
dx2a = ((100*(x[1,i-1]-40.4))/(x[1,i-1]+300*(1.7)))*dt - ((0.064*PAR*(x[1,i-1]-40.4))/(x[1,i-1]+80.8))*dt
dx3a = ((0.064*PAR*(x[1,i-1]-40.4))/(x[1,i-1]+80.8))*dt - min(50, x[3,i-1])*dt
dx4a = min(50, x[3,i-1])*dt # take the smaller value of 50 or the amount of G3P on the last step
# Calculate the Euler estimated quantities at the current time-point
x1b = x[1,i-1]+dx1a # CO2in
x2b = x[2,i-1]+dx2a # PGA
x3b = x[3,i-1]+dx3a # G3P
x4b = x[4,i-1]+dx4a # sucrose
# Calculate the change in quantities based on the Euler estimated quantities
dx1b = 1.56*g*(CO2out-x1b)*dt - ((100*(x1b-40.4))/(x1b+300*(1.7)))*dt
dx2b = ((100*(x1b-40.4))/(x1b+300*(1.7)))*dt - ((0.064*PAR*(x1b-40.4))/(x1b+80.8))*dt
dx3b = ((0.064*PAR*(x1b-40.4))/(x1b+80.8))*dt - min(50, x3b)*dt
dx4b = min(50, x3b)*dt # take the smaller value of 50 or the amount of G3P on the last step
# average the two rates of change to determine the improved Euler (RK2) estimate of the change
dx1 = (dx1a+dx1b)/2
dx2 = (dx2a+dx2b)/2
dx3 = (dx3a+dx3b)/2
dx4 = (dx4a+dx4b)/2
# calculate the improved Euler (RK2) estimate of each quantity at the new time-point
x[1,i] = x[1,i-1]+dx1 # CO2in
x[2,i] = x[2,i-1]+dx2 # PGA
x[3,i] = x[3,i-1]+dx3 # G3P
x[4,i] = x[4,i-1]+dx4 # sucrose
}
plot(t, x[1,], type="l", lty=1, main="Photosynthesis Time Course", xlab="t/s", ylab="Conc/M",
ylim=c(0,375), bty="l")
lines(t, x[2,], lty=2)
lines(t, x[3,], lty=3)
lines(t, x[4,], lty=4)
text(8, 325, "CO2in")
text(8, 185, "PGA")
text(8, 25, "G3P")
text(8, 115, "Sucrose")
---------------------------------------------
# Using the "rk4" Runge-Kutta 4 solver function in the "odesolve" package
psmodel <- function(t, x, parms) {
ci <- x[1] # CO2 internal
pga <- x[2] #
g3p <- x[3] # glucose 3 phoshpate
suc <- x[4] # sucrose
with(as.list(parms),{
import <- approx(signal$times, signal$import, t)$y
dci <- 1.56*g*(CO2out-x1b)*dt - ((100*(x1b-40.4))/(x1b+300*(1.7)))*dt
dpga <- c*s*p - d*k*p
dg3p <- e*p*k - f*k
dsuc <- e*p*k - f*k
res<-c(ds, dp, dk)
list(res)
})
}
parms = c(CO2out = 360; g = 1; kdiff = 1.56*g; PAR = 300)