A general method for parameter estimation in light-response models.

Selecting appropriate initial values is critical for parameter estimation in nonlinear photosynthetic light response models. Failed convergence often occurs due to wrongly selected initial values when using currently available methods, especially the kind of local optimization. There are no reliable methods that can resolve the conundrum of selecting appropriate initial values. After comparing the performance of the Levenberg-Marquardt algorithm and other three algorithms for global optimization, we develop a general method for parameter estimation in four photosynthetic light response models, based on the use of Differential Evolution (DE). The new method was shown to successfully provide good fits (R(2) > 0.98) and robust parameter estimates for 42 datasets collected for 21 plant species under the same initial values. It suggests that the DE algorithm can efficiently resolve the issue of hyper initial-value sensitivity when using local optimization methods. Therefore, the DE method can be applied to fit the light-response curves of various species without considering the initial values.

Photosynthesis is one of the most important biological processes involved in plant growth, and the rate of photosynthesis rate can be affected by a list of factors, such as temperature, CO2 concentration and light intensity12. In particular, many nonlinear models have been developed for describing the rate of photosynthesis in response to the change of irradiance345. There are three types of photosynthetic light response models currently available in literature: exponential, rectangular hyperbola and nonrectangular hyperbola.

Exponential Model6:

where is the net photosynthetic rate, I the irradiance, the initial quantum efficiency, Amax the net light saturated photosynthetic rate, and the dark respiration rate.

Rectangular Hyperbolic Model78:

Nonrectangular Hyperbola Model9:

where is a curvature parameter.

Modified Rectangular Hyperbola Model1011:

where β and γ are constants.

Selecting appropriate initial values is crucial for fitting these nonlinear models and estimating their parameters. For biologically meaningful parameters, using the estimates directly from experiments would be a first try for assigning the initial values. However, for model specific parameters that cannot be measured in experiments, choosing appropriate initial values when fitting the model could become problematic. For example, two of the parameters in Eq. (4) are constants without direct biological meanings. Even for those parameters with direct biological meanings, assigning their initial values in the model fitting could still be time consuming. For instance, all the three parameters in the exponential model (Eq. 1) have sound biological meanings, but only Amax can be directly measured; the other two parameters, a and R, require additional procedures to derive.

There are many statistical software packages currently available for fitting nonlinear models; to list a few, SAS, R, GenStat, MatLab, Sigmaplot, and SPSS. The built-in functionality to fit nonlinear models in these software packages normally implement algorithms of local optimization for regression and parameter estimation, such as the Gauss-Newton, Levenberg-Marquardt and Nelder-Mead algorithms. As local optimization methods are highly sensitive to the setting of initial values, failed convergence often occurs due to wrongly assigned initial values. However, there are no reliable methods available for selecting initial values5.

To develop effective methods for resolving the problem of assigning appropriate initial values when using these nonlinear models, we need to go beyond local optimization. Genetic algorithms and simulated annealing algorithms have been widely applied to search global optimal solutions in multi-dimensional space121314151617. In particular, recent comparisons have suggested that the Differential Evolution (DE) algorithm performs better than other global optimization methods such as the genetic algorithm and simulated annealing1819. There are four main steps for implementing the DE algorithms: initialization, mutation, crossover and selection, which is also capable of dealing with non-differentiable problems18. Compared to local optimization methods, global optimization algorithms are less sensitive to initial values2021. However, the application of global optimization algorithms to parameter estimation in photosynthetic light-response models is lacking.

The objectives of this study are to: (1) test the performance of the Levenberg-Marquardt (LM) algorithm, a local optimization method, and three global optimization algorithms: Genetic Algorithm (GA), Generalized Simulated Annealing (GSA) and the Differential Evolution (DE) algorithm; (2) develop a reliable method for fitting photosynthetic light-response models using the global optimization algorithm and simultaneously solve the issue of hyper initial-value sensitivity in local optimization algorithms.

Results and Discussion

Model comparison

According to calculated values of AIC (Tables 1 and 2), the modified rectangular and the nonrectangular hyperbolic models performed better than the exponential and rectangular hyperbolic models. The nonrectangular model is more flexible for its additional curvature parameter (; when θ = 0, it becomes the rectangular model522. For example, the estimates of for the photosynthetic datasets of Brassica rapa var. chinensis (L.) Kitam., Camellia sinensis (L.) Kuntze and Impatiens balsamina L. were all close to zero, the other three parameters of the nonrectangular model were also equal to those estimated from the rectangular model (Supplementary Tables S1 and S2).

Table 1

Akaike’s information criterion (AIC) and R values of each model using the Differential Evolution algorithm.

Species	Exponential		Rectangular		Nonrectangular		Modified Rectangular
	R square	AIC	R square	AIC	R square	AIC	R square	AIC
1	0.999	−21.42	0.998	−13.58	0.999	−23.14	0.999	−15.01
	0.998	−25.65	0.995	−14.58	0.998	−21.34	0.998	−22.42
2	0.987	−17.50	0.974	−9.90	0.98	−11.25	0.990	−18.15
	0.961	−14.18	0.953	−12.15	0.96	−11.12	0.970	−14.98
3	0.979	−10.42	0.995	−25.89	0.995	−23.89	0.994	−22.09
	0.991	−19.81	0.998	−34.71	0.998	−32.71	0.997	−29.49
4	0.965	0.28	0.985	−9.69	0.985	−7.69	0.985	−7.55
	0.988	−6.87	0.994	−14.62	0.994	−12.62	0.994	−12.10
5	0.978	−18.15	0.960	−11.71	0.99	−21.75	0.975	−14.85
	0.994	−18.38	0.999	−37.09	1.00	−35.06	0.999	−37.72
6	0.990	−36.34	0.992	−37.97	1.00	−41.01	0.989	−33.37
	0.997	−34.44	0.992	−25.48	1.00	−29.27	0.998	−35.79
7	0.999	−25.26	0.998	−16.01	1.00	−32.13	0.999	−22.95
	0.997	−11.62	0.998	−13.86	1.00	−12.09	0.998	−13.15
8	0.986	−20.87	0.996	−38.73	1.00	−36.62	0.996	−35.76
	0.973	−7.89	0.984	−14.25	0.98	−12.19	0.982	−10.89
9	0.996	−26.03	0.990	−18.43	0.99	−20.07	0.999	−37.43
	0.998	−28.81	0.998	−25.38	1.00	−34.55	1.000	−38.14
10	0.999	−23.71	0.998	−17.03	0.999	−22.01	0.999	−21.12
	0.995	−5.62	0.994	−2.39	0.996	−5.40	0.995	−3.93

Table 2

Akaike’s information criterion (AIC) and R values of each model using the Differential Evolution algorithm.

Species	Exponential		Rectangular		Nonrectangular		Modified Rectangular
	R square	AIC	R square	AIC	R square	AIC	R square	AIC
11	0.991	−19.48	0.975	−8.64	0.988	−14.65	0.993	−21.20
	0.958	−10.46	0.919	−3.32	0.965	−10.52	0.984	−19.03
12	0.992	−8.97	0.996	−16.51	0.996	−14.83	0.995	−13.51
	0.998	−17.39	0.999	−25.17	0.999	−25.08	0.998	−20.70
13	0.999	−10.09	1.000	−16.75	1.00	−15.45	1.000	−14.75
	0.999	−19.33	0.999	−14.16	1.00	−14.69	0.998	−7.84
14	0.997	−33.31	0.999	−41.85	1.00	−44.65	0.992	−23.31
	0.992	−22.77	0.994	−24.92	0.99	−23.36	0.996	−26.88
15	0.998	−4.91	0.996	1.26	1.00	−5.05	0.998	−1.05
	0.998	1.95	0.993	9.79	1.00	−1.34	0.996	6.48
16	0.995	−3.05	0.984	10.86	0.99	0.32	0.997	−7.49
	0.995	−11.99	0.988	−0.58	0.99	−5.62	0.994	−7.58
17	0.995	−28.86	0.994	−26.81	0.996	−30.64	0.995	−27.36
	0.997	−30.19	0.995	−22.64	0.999	−37.67	0.997	−27.04
18	0.994	−30.15	0.994	−31.26	0.995	−31.42	0.994	−28.23
	0.982	−28.23	0.984	−29.42	0.984	−27.41	0.982	−26.21
19	0.996	−18.87	0.998	−27.78	1.00	−36.33	0.999	−28.39
	1.000	−32.57	0.996	−6.78	1.00	−23.79	1.000	−28.34
20	0.987	−10.90	0.996	−24.77	0.996	−22.77	0.993	−17.47
	0.998	−31.17	0.997	−27.28	0.999	−45.01	0.998	−28.24
21	0.991	−14.11	0.968	−3.82	0.99	−11.83	0.985	−7.63
	0.997	−10.88	0.989	−0.85	1.00	−21.03	0.995	−5.20

Comparison of local and global optimization algorithms

Not surprisingly, failed convergence often occurred from wrongly assigned initial values when only using the LM algorithm, due to its hyper sensitivity to initial values (Table 3). In particular, the nonrectangular hyperbolic model failed to converge for 1260 times due to inappropriate initial values in the simulation (Table 3, Supplementary Table S3). There were also many compromised fits with R < 0.90 even when the LM algorithm converged successfully (Fig. 1). By contrast, the DE algorithm was less sensitive to initial values, which provided good fits with R approximately 0.98 for all the 42 datasets of 21 species (Tables 1 and 2) when using the same initial values for each photosynthetic model, as listed in Table 4. The DE also performed the best among the three global optimization methods, followed by the GSA (Fig. 2). In comparison, parameter estimates from the GA were not reliable with even negative R(Fig. 3). The speed of convergence for the GA algorithm was also slower than that of the DE and GSA algorithms. Although the best optimum solutions of the 100 simulation tests from the LM were similar to those from the DE (Fig. 4), however, it was not recommended to only use the LM as the large amount of failures occurred due to wrongly sampled initial values (Table 3). After setting the initial values to the parameter estimates from the DE algorithm, the value of R2 estimated from the LM algorithm only slightly increased (<0.05) for the rectangular and modified rectangular models; however, the values of R almost remained the same for the exponential and nonrectangular models (Fig. 5). Although local optimization methods are sensitive to the input of initial values, local optimization algorithms are generally faster and more accurate than global optimization method in a local space202324. In addition, there was no failures of convergence when the parameters of rectangular and modified rectangular models obtained from the DE algorithm were set as the initial values for running the LM algorithm. Therefore, we recommend to combine the DE and LM algorithms for estimating parameters of the rectangular and modified rectangular models if higher accuracy is needed, instead of increasing the number of iterations when running the DE algorithm.

Table 3

Performance of the Levenberg-Marquardt in the simulation test.

Models	Times of failed convergence
Exponential Model	358
Rectangular Hyperbolic Model	251
Nonrectangular Hyperbola Model	1260
Modified Rectangular Hyperbola Model	820

Figure 1

Distribution of R when using the Levenberg–Marquardt algorithm with successful convergence of the modified rectangular (A), nonrectangular (B), rectangular (C) and exponential (D) models.

Table 4

Lower and upper bounds of parameters for each model optimized by the Genetic Algorithm, Generalized Simulated Annealing and the Differential Evolution algorithms.

Models	Parameters	Lower bound	Upper bound
Exponential Model	a	0	1
	Amax	0	100
	Rd	0	100
Rectangular Hyperbolic Model	a	0	1
	Amax	0	100
	Rd	0	100
Nonrectangular Hyperbola Model	θ	0	1
	a	0	1
	Amax	0	100
	Rd	0	100
Modified Rectangular Hyperbola Model	a	0	1
	β	0	1
	γ	0	1
	Rd	0	100

The maximum number of iteration is 500.

Figure 2

Difference of R when using the Differential Evolution and Generalized Simulated Annealing algorithm for the modified rectangular (A), nonrectangular (B), rectangular (C) and exponential (D) models.

Figure 3

Values of R when using the Genetic Algorithm for the modified rectangular (A), nonrectangular (B), rectangular (C) and exponential (D) models.

Figure 4

Difference of R when using the Differential Evolution and the best optimum from the Levenberg–Marquardt algorithm for the modified rectangular (A), nonrectangular (B), rectangular (C) and exponential (D) models.

Figure 5

Difference of R when using the Differential Evolution, combined with the Levenberg–Marquardt algorithm for the modified rectangular (A), nonrectangular (B), rectangular (C) and exponential (D) models.

In conclusion, the application of the DE algorithm can not only effectively resolve the issue of hyper sensitivity of initial values in local optimization but also provide good fits for the photosynthetic models. The lower and upper bounds of initial values (Table 4) for each photosynthetic light-response model can be regarded as a standard setting when using the DE algorithm. In addition, the estimated parameters (Supplementary Tables S1 and S2) of the mainstream photosynthetic models for both herbaceous and woody plant species can be used as a reference for selecting initial values when using the traditional local optimization method.

Methods

A total of 42 datasets from 21 species (Supplementary Table S4), including both herbaceous and woody plants, were used for evaluating the performance of three global optimization algorithms: Genetic Algorithm (GA), Generalized Simulated Annealing (GSA) and the Differential Evolution (DE) algorithms, when fitting the photosynthetic light response models. For comparison, the Levenberg-Marquardt (LM) algorithm, a derivative based local optimization method, was also applied to estimate parameters of the photosynthetic models. The raw data of species 2, 5–9, 14–17, 20 and 22 (see Supplementary Table S4) were extracted from figures in publications252627282930313233343536 using the GetData Graph Digitizer 2.26 (http://getdata-graph-digitizer.com). After 500 iterations, the estimates from the three optimization methods were used to evaluate their performance. To avoid biasing the results, we used the same lower and upper bounds for parameter estimation when using the three global optimization (Table 4). For each dataset, 100 random values were drew from a uniform distribution with the same lower and upper bounds for the global optimization algorithms as initial values for testing the LM algorithm. A total of 4200 simulations were repeated for each photosynthetic model. Parameters of all the photosynthetic light-response models were estimated by the LM, GA, GSA and DE algorithms, implemented in the minpack.lm, GA, GenSA, and DEoptim packages of R37. The coefficient of determination and Akaike’s information criterion (AIC)38 were used as the goodness-of-fit index and for model comparison, respectively.

Additional Information

How to cite this article: Chen, L. et al. A general method for parameter estimation in light-response models. Sci. Rep. 6, 27905; doi: 10.1038/srep27905 (2016).