LI-6800 光响应曲线与CO2响应曲线的拟合

以 R 语言代码为例

祝介东

北京力高泰科技有限公司(2024.3.18离职)

主要内容

  • 光响应曲线拟合
    • 相关资源
    • 示例代码
  • \(CO_2\) 响应曲线拟合
    • 相关资源
    • C3 植物拟合
    • C4 植物拟合

流程

光响应曲线拟合

光响应曲线拟合所需要的测量参数


光响应曲线拟合仅需使用光合有效辐射和净光合速率两个参数。


数据表格参数名 含义
A 净光合速率
Qin 叶片表面的光合有效辐射

光响应曲线拟合

软件包

  • FitAQ: 基于非直角双曲线模型的软件包,发布于 github
  • photosynthesis:生理生态学模型和分析工具。

不足之处是仅有非直角双曲线模型可以使用。

安装

devtools::install_github("MarkusLoew/FitAQ")
install.packages('photosynthesis')


代码

photosynthesis 包使用示例


# 加载软件包
library(photosynthesis)

# 读取数据
dat = read.csv("data/dat-aq.csv")

# 进行拟合
fit = fit_photosynthesis(
  .data = dat, # 指定拟合数据
  .photo_fun = "aq_response", # 指定光响应曲线拟合
  .vars = list(.A = A, .Q = Qin), # 指定拟合所需数据
)

# 拟合结果的统计汇总
summary(fit)

photosynthesis 包的拟合结果



Formula: .A ~ marshall_biscoe_1980(Q_abs = .data[[".Qabs"]], k_sat, phi_J, 
    theta_J) - Rd

Parameters:
         Estimate Std. Error t value Pr(>|t|)    
k_sat   9.1669265  0.0305156  300.40   <2e-16 ***
phi_J   0.0488003  0.0005419   90.05   <2e-16 ***
theta_J 0.7503171  0.0063474  118.21   <2e-16 ***
Rd      0.4296552  0.0223950   19.18   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.08653 on 600 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 1.49e-08

photosynthesis 计算补偿点


dt_lcp = as.data.frame(t(coef(fit)))
lcp = ((dt_lcp$Rd) * (dt_lcp$Rd * dt_lcp$theta_J - dt_lcp$k_sat) / (dt_lcp$phi_J * (dt_lcp$Rd - dt_lcp$k_sat)))
lcp
[1] 8.912453


MARSHALL and BISCOE (1980):

\[\begin{align} & aP^2 +bP_n+c=0 \\ & a = \theta \\ & b = -(P_{max}+\alpha I -\theta R_d) \\ & c = \alpha I(P_{max} - (1-\theta) R_d) - R_d P_{max} \end{align}\]

photosynthesis 拟合结果作图示例

library(ggplot2) # 加载软件包

# 计算拟合的A,用于作图----------------
b = coef(fit)
pred_Qin = seq(0, 1000, length.out = 100)
pred_A = marshall_biscoe_1980(
    Q_abs = pred_Qin,
    k_sat = b["k_sat"],
    b["phi_J"],
    b["theta_J"]
    ) - b["Rd"]
df_predict = data.frame(A = pred_A, Qin = pred_Qin)
#----------------------------------------------------

ggplot(mapping = aes(x = Qin, y = A)) +
  geom_point(data = dat, color = 'lightgray', size = 2)+
  geom_line(data = df_predict, color = 'blue', linewidth = 2) +
  labs(
    x = expression("Irradiance (" * mu * mol ~ m^{-2} ~ s^{-1} * ")"),
    y = expression(A[net] ~ "(" * mu * mol ~ m^{-2} ~ s^{-1} * ")")
  ) + 
  theme_bw()

photosynthesis 结果作图示例

photosynthesis 光响应曲线拟合结果作图

\(CO_2\) 响应曲线的拟合

\(CO_2\) 响应曲线拟合所需要的参数

二氧化碳响应曲线所需参数较多,以 plantecophys 为例,使用的参数如下表:

数据表格参数名 含义
A 净光合速率
Qin 叶片表面的光合有效辐射
TleafCnd 计算使用的叶片温度
Ci 胞间二氧化碳浓度
Pa 大气压强

\(CO_2\) 响应曲线拟合软件包或代码


C3 二氧化碳响应曲线拟合


安装

install.packages("plantecophys")
install.packages("photosynthesis")

\(CO_2\) 响应曲线拟合软件包或代码


C4 二氧化碳响应曲线拟合


安装

devtools::install_github('eloch216/PhotoGEA')
#install.packages("lattice")

一些参考的代码1

C3 植物的 \(CO_2\) 响应曲线拟合

photosynthesis 软件包的拟合

library(photosynthesis)

# 读入数据
dat = read.csv("data/rice.csv")

# 一定要手动转换为K氏温度
dat$T_leaf <- dat$Tleaf + 273.15

# 拟合
fit = fit_aci_response(dat, 
  varnames = list(
    A_net = "A",
    T_leaf = "T_leaf",
    C_i = "Ci",
    PPFD = "Qin"
  ),
  P = mean(dat$Pa)
  )

# 拟合参数
fit[[1]]

photosynthesis 包的拟合结果



Num V_cmax V_cmax_se J_max J J_se V_TPU V_TPU_se R_d R_d_se cost citransition1 citransition2 V_cmax_pts J_max_pts V_TPU_pts alpha alpha_g gamma_star25 Ea_gamma_star K_M25 Ea_K_M g_mc25 Ea_g_mc Oconc theta_J
78 0 31.7537 0.226247 52.14723 50.53414 0.0235851 1000 NA -1.678633 0.0372401 1.886089 370.7022 729.2324 80 70 0 0.24 0 42.75 37830 718.4 65508.28 0.08701 0 21 0.85

photosynthesis 包的拟合结果作图



fit[[2]]

plantecopys 包的拟合代码

library(plantecophys)

# Read in your data
dat = read.csv("data/rice.csv")

fit = fitaci(
  dat,
  varnames = list(
    ALEAF = "A",
    Tleaf = "TleafCnd",
    Ci = "Ci",
    PPFD = "Qin"
  ),
  Patm = mean(dat$Pa),
  alpha = 0.5,
  theta = 1
)
fit$pars
fit$nlsfit

plantecopys 包的拟合结果

       Estimate Std. Error
Vcmax 27.906434 0.18557176
Jmax  48.109668 0.18259949
Rd     1.715652 0.03378911
Nonlinear regression model
  model: ALEAF ~ acifun_wrap(Ci, PPFD = PPFD, Vcmax = Vcmax, Jmax = Jmax,     Rd = Rd, Tleaf = Tleaf, Patm = Patm, TcorrectVJ = Tcorrect,     alpha = alpha, theta = theta, gmeso = gmeso, EaV = EaV, EdVC = EdVC,     delsC = delsC, EaJ = EaJ, EdVJ = EdVJ, delsJ = delsJ, Km = Km,     GammaStar = GammaStar)
   data: data
 Vcmax   Jmax     Rd 
27.906 48.110  1.716 
 residual sum-of-squares: 3.183

Number of iterations to convergence: 8 
Achieved convergence tolerance: 6.559e-06

plantecophys 包的拟合结果作图

plot(fit)

C4 植物的 \(CO_2\) 响应曲线拟合

PhotoGEA 数据导入失败

PhotoGEA数据导入失败,则需开启 excel 公式的迭代功能:

PhotoGEA 包的拟合代码

library(PhotoGEA)
library(lattice) # 快速作图查看数据的依赖


# 导入数据,时间戳使用'time',非必须
df =  read_gasex_file('data/PhotoGEA.xlsx', 'time')

# 阿伦尼乌斯方程计算 C4 photosynthetic 参数
df =  calculate_arrhenius(df, c4_arrhenius_von_caemmerer)

# 计算总的压强 bar
df <- calculate_total_pressure(df)

# 指定叶肉导度,如未知,则是 inf
df <- set_variable(
  df, 'gmc',
  value = Inf
)


# 计算叶肉 分压 PCm
df <- apply_gm(
  df,
  'C4' # Indicate C4 photosynthesis
)

# 拟合
c4_fit = fit_c4_aci(df, Ca_atmospheric = 400)


# 提取要查看的参数列
columns_for_viewing <-
  c('Rd_at_25', 'Vcmax_at_25', 'Vpmax_at_25')

# 查看参数结果
c4_fit_parameters <-
 c4_fit$parameters[ , columns_for_viewing, TRUE]

print(c4_fit_parameters)

PhotoGEA 包的拟合结果


部分拟合结果

  Rd_at_25 [fit_c4_aci] (micromol m^(-2) s^(-1))
1                                       10.58733
  Vcmax_at_25 [fit_c4_aci] (micromol m^(-2) s^(-1))
1                                          53.58885
  Vpmax_at_25 [fit_c4_aci] (micromol m^(-2) s^(-1))
1                                          186.3173


缩写含义

PhotoGEA 包的拟合结果作图



# 快速作图查看数据
xyplot(
  A + A_fit ~ Ci,
  data = c4_fit$fits$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste('Intercellular CO2 concentration [', c4_fit$fits$units$Ci, ']'),
  ylab = paste('Net CO2 assimilation rate [', c4_fit$fits$units$A, ']')
)

PhotoGEA 包的拟合结果作图

PhotoGEA 包的拟合结果作图



# Plot the C4 A-PCm fits (including limiting rates)
xyplot(
  A + Apc + Ar +A_fit ~ PCm,
  data = c4_fit$fits$main_data,
  type = 'l',
  auto.key = list(space = 'right'),
  grid = TRUE,
  xlab = paste('Mesophyll CO2 pressure [', c4_fit$fits$units$PCm, ']'),
  ylab = paste('Net CO2 assimilation rate [', c4_fit$fits$units$A, ']'),
  par.settings = list(
    superpose.line = list(col = multi_curve_colors()),
    superpose.symbol = list(col = multi_curve_colors())
  )
)

PhotoGEA 包的拟合结果作图

C4 植物手写代码及作图

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

C4 植物拟合代码及作图

     Vcmax    Vpmax     Jmax ci_trans Tleaf    Vpr  Kc  Ko Kp
1 161.9322 117.1201 243.0344      150    25 48.281 650 450 80

参考文献

Duursma, Remko A. 2015. “Plantecophys - An R Package for Analysing and Modelling Leaf Gas Exchange Data.” PLOS ONE 10 (11): e0143346. https://doi.org/10.1371/journal.pone.0143346.
Lochocki, Edward B. 2023. PhotoGEA: Photosynthetic Gas Exchange Analysis.
MARSHALL, B., and P. V. BISCOE. 1980. “A Model for C3 Leaves Describing the Dependence of Net Photosynthesis on Irradiance.” Journal of Experimental Botany 31 (1): 29–39. https://doi.org/10.1093/jxb/31.1.29.
Wickham, Hadley, and Garrett Grolemund. 2023. R for Data Science. " O’Reilly Media, Inc.".
Zhou, Haoran, Erol Akçay, and Brent R. Helliker. 2019. “Estimating C4 Photosynthesis Parameters by Fitting Intensive A/Ci Curves.” Photosynthesis Research 141 (2): 181–94. https://doi.org/10.1007/s11120-019-00619-8.

感谢您的参与

感谢使用 LI-6800

关注公众号后获取更多精彩内容