library (minpack.lm)
library (Metrics)
##########general inputs for temp and vpr#####
T.ppdk<- c (10 ,15 ,20 ,25 ,30 ,35 ,40 )
# From Serrano-Romero and Cousins. 2020 New Phytologist
ppdk<- c (13.2 ,21.6 ,32.4 ,48 ,70.4 ,98.8 ,142 )
# From Serrano-Romero and Cousins. 2020 New Phytologist
ppdk.df<- data.frame (T.ppdk, ppdk); rm (T.ppdk, ppdk)
ppdk.df$ T2<- ppdk.df$ T.ppdk^ 2
ppdk.mod<- lm (log (ppdk) ~ T.ppdk + T2, ppdk.df)
mod.sum<- summary (ppdk.mod)
intercept<- mod.sum$ coefficients[1 ,1 ]
CTleaf.slope<- mod.sum$ coefficients[2 ,1 ]
CTleaf2.slope<- mod.sum$ coefficients[3 ,1 ]
temp.ppdk.df<- data.frame (intercept, CTleaf.slope, CTleaf2.slope)
############################################
# load implementation of C4 ACi curve fitter
# (based on von Caemmerer 2000 and function AciC4 from plantecophys)
############################################
# Km,GammaStar Optionally, provide Michaelis-Menten coefficient for
# Farquhar model, and Gammastar. If not provided, they are calculated
# with a built-in function of leaf temperature
# Rd Day respiration rate (mu mol m-2 s-1), optional (if not provided,
# calculated from Tleaf, Rd0, Q10 and TrefR). Must be a positive value
# (an error occurs when a negative value is supplied
O2= 210 #O2 Mesophyll O2 concentration
FRM= 0.5
#FRM Fraction of day respiration that is mesophyll respiration (Rm)
alpha= 0 #alpha Fraction of PSII activity in the bundle sheath (-)
Q10= 2 #Q10 T-dependence parameter for Michaelis-Menten coefficients
x= 0.4 #x Partitioning factor for electron transport
THETA= 0.7 #THETA Shape parameter of the non-rectangular hyperbola
low_gammastar <- 1.93e-4
# Half the reciprocal for Rubisco specificity (NOT CO2 compensation
# point)
#Vpr=80 #Vpr PEP regeneration (mu mol m-2 s-1)
gbs= 3e-3 #gbs Bundle sheath conductance (mol m-2 s-1)
######enzyme-limited photosynthesis function#####
# enzyme-limited photosynthesis function
A.enzyme.func= a ~ (- (- (((pmin (ci * Vpmax/ (ci + Kp), Vpr)) - Rm +
gbs * ci) + (Vcmax - Rd) + gbs * K + ((alpha/ 0.047 ) *
(low_gammastar * Vcmax + Rd * Kc/ Ko)))) -
sqrt ((- (((pmin (ci * Vpmax/ (ci + Kp), Vpr)) -
Rm + gbs * ci) + (Vcmax - Rd) + gbs * K +
((alpha/ 0.047 ) * (low_gammastar * Vcmax + Rd *
Kc/ Ko))))^ 2 - 4 * a.c * ((Vcmax - Rd) * ((pmin (ci *
Vpmax/ (ci + Kp), Vpr)) - Rm + gbs * ci) -
(Vcmax * gbs * low_gammastar * O2 + Rd *
gbs * K))))/ (2 * a.c)- Rd
#######enzyme-limited predicted photosynthetic function ########
Pre_enz<- function (ci) {
(- (- (((pmin (ci * Vpmax/ (ci + Kp), Vpr)) - Rm +
gbs * ci) + (Vcmax - Rd) + gbs * K + ((alpha/ 0.047 ) *
(low_gammastar * Vcmax + Rd * Kc/ Ko)))) -
sqrt ((- (((pmin (ci * Vpmax/ (ci + Kp), Vpr)) -
Rm + gbs * ci) + (Vcmax - Rd) + gbs * K +
((alpha/ 0.047 ) * (low_gammastar * Vcmax +
Rd * Kc/ Ko))))^ 2 - 4 * a.c * ((Vcmax - Rd) *
((pmin (ci * Vpmax/ (ci + Kp), Vpr)) - Rm +
gbs * ci) - (Vcmax * gbs * low_gammastar *
O2 + Rd * gbs * K))))/ (2 * a.c)- Rd
}
######## light-limited photosynthesis function ##########
A.light.func = a~ (- (- ((x * ((1 / (2 * THETA)) *
(Qp2 + Jmax - sqrt ((Qp2 + Jmax)^ 2 - 4 * THETA *
Qp2 * Jmax)))/ 2 - Rm + gbs * ci) + ((1 - x) *
((1 / (2 * THETA)) * (Qp2 + Jmax - sqrt ((Qp2 +
Jmax)^ 2 - 4 * THETA * Qp2 * Jmax)))/ 3 - Rd) +
gbs * (7 * low_gammastar * O2/ 3 ) + alpha *
low_gammastar/ 0.047 * ((1 - x) * ((1 / (2 *
THETA)) * (Qp2 + Jmax - sqrt ((Qp2 + Jmax)^ 2 -
4 * THETA * Qp2 * Jmax)))/ 3 + 7 * Rd/ 3 ))) -
sqrt ((- ((x * ((1 / (2 * THETA)) * (Qp2 + Jmax -
sqrt ((Qp2 + Jmax)^ 2 - 4 * THETA * Qp2 * Jmax)))/
2 - Rm + gbs * ci) + ((1 - x) * ((1 / (2 * THETA)) *
(Qp2 + Jmax - sqrt ((Qp2 + Jmax)^ 2 - 4 * THETA *
Qp2 * Jmax)))/ 3 - Rd) + gbs * (7 * low_gammastar *
O2/ 3 ) + alpha * low_gammastar/ 0.047 * ((1 - x) *
((1 / (2 * THETA)) * (Qp2 + Jmax - sqrt ((Qp2 +
Jmax)^ 2 - 4 * THETA * Qp2 * Jmax)))/ 3 +
7 * Rd/ 3 )))^ 2 - 4 * a.j * (((x * ((1 / (2 *
THETA)) * (Qp2 + Jmax - sqrt ((Qp2 + Jmax)^ 2 -
4 * THETA * Qp2 * Jmax)))/ 2 - Rm + gbs *
ci) * ((1 - x) * ((1 / (2 * THETA)) *
(Qp2 + Jmax - sqrt ((Qp2 + Jmax)^ 2 -
4 * THETA * Qp2 * Jmax)))/ 3 - Rd)) -
gbs * low_gammastar * O2 * ((1 - x) *
((1 / (2 * THETA)) * (Qp2 + Jmax -
sqrt ((Qp2 + Jmax)^ 2 - 4 * THETA *
Qp2 * Jmax)))/ 3 + 7 * Rd/ 3 ))))/ (2 * a.j)- Rd
#########light-predicted photosynthesis function ######
pre_light = function (ci){
(- (- ((x * ((1 / (2 * THETA)) * (Qp2 + Jmax -
sqrt ((Qp2 + Jmax)^ 2 - 4 * THETA * Qp2 *
Jmax)))/ 2 - Rm + gbs * ci) + ((1 - x) *
((1 / (2 * THETA)) * (Qp2 + Jmax - sqrt ((Qp2 +
Jmax)^ 2 - 4 * THETA * Qp2 * Jmax)))/ 3 -
Rd) + gbs * (7 * low_gammastar * O2/ 3 ) +
alpha * low_gammastar/ 0.047 * ((1 - x) *
((1 / (2 * THETA)) * (Qp2 + Jmax -
sqrt ((Qp2 + Jmax)^ 2 - 4 * THETA * Qp2 *
Jmax)))/ 3 + 7 * Rd/ 3 ))) - sqrt ((- ((x *
((1 / (2 * THETA)) * (Qp2 + Jmax -
sqrt ((Qp2 + Jmax)^ 2 - 4 * THETA *
Qp2 * Jmax)))/ 2 - Rm + gbs * ci) +
((1 - x) * ((1 / (2 * THETA)) * (Qp2 +
Jmax - sqrt ((Qp2 + Jmax)^ 2 - 4 *
THETA * Qp2 * Jmax)))/ 3 - Rd) + gbs *
(7 * low_gammastar * O2/ 3 ) + alpha *
low_gammastar/ 0.047 * ((1 - x) * ((1 / (2 *
THETA)) * (Qp2 + Jmax - sqrt ((Qp2 +
Jmax)^ 2 - 4 * THETA * Qp2 * Jmax)))/ 3 +
7 * Rd/ 3 )))^ 2 - 4 * a.j * (((x * ((1 / (2 *
THETA)) * (Qp2 + Jmax - sqrt ((Qp2 +
Jmax)^ 2 - 4 * THETA * Qp2 * Jmax)))/ 2 -
Rm + gbs * ci) * ((1 - x) * ((1 / (2 *
THETA)) * (Qp2 + Jmax - sqrt ((Qp2 +
Jmax)^ 2 - 4 * THETA * Qp2 * Jmax)))/ 3 -
Rd)) - gbs * low_gammastar * O2 * ((1 - x) *
((1 / (2 * THETA)) * (Qp2 + Jmax - sqrt ((Qp2 +
Jmax)^ 2 - 4 * THETA * Qp2 * Jmax)))/ 3 +
7 * Rd/ 3 ))))/ (2 * a.j)- Rd
}
############################################
# set parameters
############################################
# id_c4 <- 'KBS_Zmay_1' # identifier for C4 curve of interest
b_tresp_R <- 0.1781075
# value for b parameter in temperature response curve for respiration
# (From Smith and Dukes (2017), doi = 10.1111/gcb.13735)
c_tresp_R <- - 0.00179152
# value for c parameter in temperature response curve for respiration
# (From Smith and Dukes (2017), doi = 10.1111/gcb.13735)
############################################
# read in ACi curve
############################################
aci_data<- read.csv ('data/fake-racir.csv' )
aci_data$ a = aci_data$ A
aci_data$ ci = aci_data$ Ci
# set fit parameters (see documentation for function AciC4 from # plantecophys)
############################################
# Tleaf=mean(aci_data$TleafCnd, na.rm =TRUE) # in dataframe, in CHR
# format # mean of measurement_t = 24.7
Tleaf= 25
Rd= 1
# find that article on this. (Wang et al., 2014 ; experimental botany;
# Three distinct biochemical subtypes of C4 photosynthesis)
PPFD= mean (aci_data$ Qin, na.rm = TRUE )
# Michaelis-Menten coefficients for CO2 (Kc, mu mol mol-1) and O
# (Ko, mmol mol-1) and combined (K)
Kc <- 650 * Q10^ ((Tleaf-25 )/ 10 )
Kp <- 80 * Q10^ ((Tleaf-25 )/ 10 )
Ko <- 450 * Q10^ ((Tleaf-25 )/ 10 )
K <- Kc* (1 + O2/ Ko)
Rm <- FRM* Rd # Day leaf respiration, umol m-2 s-1
Qp2 <- PPFD* 0.85 * (1 -0.15 )/ 2
# Non-rectangular hyperbola describing light effect on electron
# transport rate (J)
a.c <- 1 - (alpha* Kc)/ (0.047 * Ko)
a.j <- 1 - 7 * low_gammastar * alpha/ (3 * 0.047 )
Vpr<- exp (temp.ppdk.df$ intercept+
(temp.ppdk.df$ CTleaf.slope* Tleaf)+
(temp.ppdk.df$ CTleaf2.slope* Tleaf** 2 ))
############################################
# plot data to estimate transition point & find irregularities
############################################
ci_trans <- 150 # estimated Ci transition point
# fit the enzyme and light limited portions of the curve separately
############################################
#est.Vcmax<-30
fit_enzyme <- nlsLM (A.enzyme.func, data = aci_data,
start= list (Vcmax= 50 ,Vpmax= 150 ), control= nls.control (maxiter= 500 , minFactor= 1 / 10000 ))
# fit the enzyme limited portion of the curve (Vcmax and Vpmax)
fit_light = nlsLM (A.light.func, data= aci_data,
start = list (Jmax = 180 ), control = nls.control (maxiter = 500 , minFactor = 1 / 10000 ))
# fit the light-limted portion of the curve (Jmax)
# predict(fit_light)
######Summary of enzyme and light####
sum.enzyme<- summary (fit_enzyme)
sum.light<- summary (fit_light)
Vcmax<- sum.enzyme$ coefficients[1 ,1 ]
Vpmax<- sum.enzyme$ coefficients[2 ,1 ]
Jmax<- sum.light$ coefficients[1 ,1 ]
########Enzyme limited dataframe #####
output_enz<- data.frame (A_enz_obs= as.numeric (),
A_enz_pre= as.numeric (), no= as.numeric ())
##########Light limited dataframe #####
output_light<- data.frame (A_light_obs= as.numeric (),
A_light_pre= as.numeric (), no= as.numeric ())
##### Plot graph ######
plot (aci_data$ ci, aci_data$ a, xlab = "Ci (ppm)" , ylab = "A" , pch = 20 , col = 'blue' )
lines (aci_data$ ci, predict (fit_enzyme), col = 'red' , lwd = 2 )
lines (aci_data$ ci, predict (fit_light), col = 'brown' , lwd = 2 )
abline (v = ci_trans)
legend ("bottomright" , c ("Observed" , "enzyme limt" , "ETR limit" ),
pch = c (20 , NA , NA ), col = c ("#3300CC" , "red" , "brown" ), lwd = c (NA , 2 , 2 ))
# make a temp df for storage ####
temp.df<- data.frame (Vcmax, Vpmax,Jmax, ci_trans, Tleaf,Vpr,Kc,Ko,Kp)
temp.df