The dynamic growth model of lettuce proposed by Van Henten (1994) is modified in this study to numerically investigate the impact of spectral distribution and intensity on lettuce dry mass and yield. Using dry mass as the primary output for the model, this variable is further subdivided into structural dry mass and nonstructural dry mass, which accounts for starch, glucose, and other similar elements. The model assumes that the two categories of dry mass fully define the state of the plant and describes lettuce growth by calculating these sub-variables using the following ordinary differential equations (ODE),

Equations (1) and (2) represent the transient behavior in the structural and non-structural dry mass per unit of area ( g   m − 2 ) in response to photosynthesis. In the above equations, f p h o t = f p h o t (C C O 2 , I , T , X s d m ) is the gross canopy photosynthesis ( g   m − 2   s − 1 ) , f r e s p = f r e s p ( T , X n s d m , X s d m ) is the maintenance respiration ( g   m − 2   s − 1 ) , and r g r = r g r ( T , X n s d m , X s d m ) is the growth rate of structural material ( g   m − 2   s − 1 ) , while c α and c β describe the conversion rate of CO 2   to sugar (CH   2 O) and yield factor which is a measure of non-structural dry mass losses due to respiration and photosynthetic activities, respectively. The value for c α is the molecular weight ratio of CO   2 to CH   2 O and is set to 0.68 . According to Sweeney (1981), c β for lettuce is approximately 0.8 . The growth rate ( r g r ) refers to the rate at which non-structural materials are transformed into structural materials, i.e,

where T is the canopy temperature ( ° C ) , c g r , m a x is the saturation growth rate at 20 ° C , c γ is the growth rate coefficient, and c Q 10 , g r is the measure of growth rate sensitivity to the canopy temperature. Van Holsteijn (1981) approximated the saturation growth rate coefficient to 5 × 10 − 6   s − 1 ; Sweeney (1981) estimated the growth rate coefficient for lettuce to 1.0 . The growth rate sensitivity constant is set to 1.6 , which means that for every 10 ° C increase in the canopy temperature, the growth rate increases by a factor of 1.6 . The maintenance respiration rate is predicted through,

In Equation (4), c r e s p , s h t and c r e s p , r t represent shoot and root maintenance respiration coefficients at 25 ° C and indicate the amount of glucose consumption per structural dry material. Van Henten (1994) estimated shoot and root respiration coefficients as 3.47 × 10 − 7   s − 1 and 1.16 × 10 − 7   s − 1 , respectively. c Q 10 , r e s p is the sensitivity of maintenance respiration to canopy temperature and Van Henten (1994) assigned a value of 2.0 for this coefficient. c τ is the ratio of root dry mass to the overall dry mass of the plant, which can depend on the type of cultivation. Lorenz and Wiebe (1980)reported an average value of 0.15 for lettuce cultivated in soil, while Sakamoto and Suzuki (2015) measured an average value of 0.14 for hydroponic lettuce cultivation. Goudriaan and Monteith (1990) formulated an empirical correlation to estimate gross canopy photosynthesis,

that c K is the extinction coefficient, and for lettuce with planophile characteristics, is set to 0.9 ; c l a r is the structural leaf area ratio and Lorenz and Wiebe (1980) approximated it to 75 × 10 − 3   m 2   g − 1 ; and f p h o t , m a x is the gross CO   2 assimilation rate for a canopy with 1 square meter of effective surface area. Acock et al. (1978) presented an equation to calculate f p h o t , m a x considering the effect of CO   2 concentration and spectral irradiance integrated over wavebands within PAR as well as canopy temperature and photorespiration,

In Equation (6), ε is the light-use efficiency, I is the spectral irradiance integrated over the wavebands within PAR that regulates plant growth, g C O 2 is the conductance of canopy for the diffusion of CO   2 , c ω is the density of CO   2 that has an approximate value of 1.83 × 10 − 3   g   m − 3 (considering greenhouse temperature around 20   ∘ C), C C O 2 is the concentration of CO   2 in the greenhouse, and Γ is the CO   2 compensation point, accounting for the impact of the temperature on photosynthesis rate (Van Henten, 1994). CO   2 compensation is determined based on canopy temperature using the following correlation,

whereas c Γ is the CO   2 compensation point at 20 ° C which is 40   m L   L − 1 , and c Q 10 , Γ is the sensitivity of CO   2 compensation with canopy temperature, which Goudriaan et al. (1985) approximated it as 2.0 . Light-use efficiency is computed considering light level impact on CO   2 compensation and photorespiration (Goudriaan et al., 1985),

In Equation (8), c ε is the quantum use efficiency which is the energy required for a reduction of one mole CO   2 , and Goudriaan et al. (1985) approximated its value to be about 17.0 × 10 − 6   g   J − 1 . In this study, there is an assumption that this parameter is affected by the photon spectral distribution; therefore, its value varies depending on the lighting conditions utilized for lettuce growth. Goudriaan et al. (1985) developed a mathematical correlation for the canopy conductance for CO   2 diffusion, which is derived considering the boundary layer, stomatal, and carboxylation conductance,

where g b n d , g s t m , and g c a r represent the boundary layer, stomatal, and carboxylation conductance, respectively. Stanghellini (1987) estimated the boundary layer conductance to be 0.007   m   s − 1 at a 5°C temperature gradient, 0.1 m   s − 1 wind speed, and leaf with a characteristic length of 0.075   m . For a plant that grows in an environment without stress, Stanghellini (1987) approximated stomatal conductance to be 0.005   m   s − 1 . Carboxylation conductance is a function of canopy temperature and its value (from 5 to 40°C) is determined using the following empirical correlation,

Table 1 provides a summary for the definition of different coefficients and their numerical values within the plant growth model.

A MATLAB code was developed to find a solution for the ODE Equations (1), and (2). The code utilized experimental temperature, spectral irradiance integrated over wavebands from 400 to 750 n m , and CO   2 concentration as inputs to compute the two sub-variable dry masses as outputs. Input and output data were extracted from experiments that investigated the impact of the photon spectrum on production of lettuce ‘Rouxai’ growth by Meng and Runkle (2019); Meng et al. (2019), and Meng et al. (2020). Meng and Runkle (2019) carried out three replications with a PPFD of 100 and 180 μ m o l   m − 2   s − 1 during 0-2 and 2-3 days, respectively. After that, the seedlings were grown under various LED treatments with a 24-hour photoperiod, and the temperature was set to 23 ° C. Fresh and dry mass data were obtained using destructive tests for plants harvested on day 10. Similarly, Meng et al. (2019) performed experiments three times at a total photon flux density of 180 μ m o l   m − 2   s − 1 . The seedlings were transplanted into a hydroponic system with a 20-hour photoperiod, an average air temperature of 21.1   ∘ C, average CO   2 concentration of 402 m L   L − 1 , and relative humidity ranging from 41 % to 70%. Meng et al. (2020) conducted two replications at a temperature of 20 ° C, a total photon flux density of 50 μ m o l   m − 2   s − 1 , and 24-hour photoperiod. The next day, the temperature, the photoperiod, and total photon flux density were set to 22 ° C, 20 hours, and 180 μ m o l   m − 2   s − 1 , respectively. On the fourth day, the seedlings were exposed to nine different light-quality treatments under the same controlled conditions. For the first replication, the average temperature, relative humidity, and CO   2 concentration was 22.4 °C, 410 m L   L − 1 , and 34 % , and in the second replication, these parameters were 22.4 °C, 410 mL L-1 and 35%, respectively. Table 2 represents the experimental data for red-leaf lettuce ‘Rouxai’ under different light treatment conditions, while Figure 1 illustrate the variation of incoming light spectra.

Using integrated spectral irradiance, CO   2 concentration, and temperature, the Van Henten (1994) growth model is used with the parameters in Table 1 to predict dry mass for different light treatment experiments. Figure 2 represents the necessity of considering the impact of spectral distribution on lettuce growth by showing the difference between the experimental data and the plant growth model. For example, the growth model predicted only 6 of the 20 lighting treatments to be within 25 % of the actual dry mass values.

Numerical error is the measure of a difference between the experimental dry mass of lettuce and the model prediction using the suggested value for c ε (Van Henten, 1994), which was 17.0 × 10 − 6   g   J − 1 . With the assumption of unvarying c ε for different experiments, the dynamic growth model does not accurately predict for various light treatment experiments other than typical greenhouse light conditions.

As mentioned earlier, using a constant c ε led to a considerable error in lettuce dry mass prediction; however, as we will show, a varying c ε depending on spectral distribution allows prediction in good accordance with experimental data. Knowing the dry mass of lettuce for different experiments, the solver tries to find a value for c ε that allows a prediction of a state variable in good accordance with experimental data. These values will be used in the next section to develop a model that predicts the impact of the photon spectrum on c ϵ , and eventually on plant growth. Figure 3 compares the lettuce dry mass predicted by the model with the experimental dry mass for the investigated light treatments. It is inferred from Figure 3 that the solver is capable of finding a quantum use efficiency for each experiment that leads to an accurate prediction of the lettuce state variable (dry mass) on the day of harvest.

The aim of this study is to develop a model that predicts the quantum use efficiency ( c ε ) as a function of spectral photon flux density or integrated spectral irradiance within the PAR+FR waveband (i.e., 400-750 nm). A linear regression approach proves to be an invaluable tool for generating a basic model for obtaining weights for different features and establishing a simple mathematical model in the form of c ϵ = ∑ i = 1 4 ( ω i F i ) where ω i and F i correspond to the weight (coefficient) and the value for the i   th input feature, respectively. Since the incoming spectrum is a continuous function ( Figure 1 ), an idea was devised to generate discrete features based on the continuous distribution of the spectrum for the linear regression model, as shown in Figure 4 . This idea involves dividing the incoming spectrum into four segments, which in the context of this study, corresponding to four distinctive wavebands in PAR+FR: blue (400-500 nm), green (500-600 nm), red (600-700 nm), and far-red (700-750 nm). A discrete value is assigned for each segment based on the integral of spectral distribution. Subsection 3.2 provides a detailed description of how these discrete values are obtained for each spectrum.

In addition, it is possible to investigate the interaction between different light wavebands on quantum use efficiency through the regression model in the form of c ¯ ϵ = ∑ i = 1 4 ( ω i F i ) + ∑ j ( ω j   ' F j   ' ) where F j   ' can be defined as the multiplication of two fraction ratios, e.g., the blue and green wavebands. Figure 4 provides an overview of the development of empirical correlation for light-use efficiency based on the continuous spectrum distribution.

As previously stated, the dynamic growth model’s efficiency can be improved by developing a model for c ε that takes into account the impact of spectral distribution. This can be accomplished by creating a function for c ε that is dependent on the photon flux density or integrated spectral irradiance ratio for 100-nm wavebands in PAR and 50-nm FR waveband. Calculation of these ratios based on spectral photon flux distribution is more convenient since Meng et al. (2020) assigned a label based on photon flux density treatments with different wavebands. Considering lighting treatment 18, or “B20G20R120FR20” as an example, the photon flux density for blue and green light is 20   μ m o l   m − 2   s − 1 , red light is 120   μ m o l   m − 2   s − 1 , and far-red light is 20   μ m o l   m − 2   s − 1 . The photon flux density for different wavebands represents the areas under the curve in Figure 5 . Therefore, photon flux density ratios for the blue, green, and far-red wavebands are 20 180 , while for the red waveband is equal to 120 180 . Computation of the ratios for integrated spectral irradiation of the PAR+FR wavebands is different, since spectral irradiance is a measure of the energy carried by a photon. Therefore, integrated spectral irradiance (I) for a specific spectrum of PAR+FR is determined by calculating the energy for each wavelength through multiplication of wavelength energy and its number of photons and integrating those over the specific waveband. For “B20G20R120FR20” as an example, the integrated spectral irradiance for blue, green, red and far-red wavebands are 5.57   W   m − 2 , 4.57   W   m − 2 , 21.74   W   m − 2 , and 3.51   W   m − 2 , respectively. Based on the integrated spectral irradiance values computed for different wavebands, the integrated spectral irradiance intensity ratios for blue, green, red, and far-red are 5.57 35.38 , 4.57 35.38 , 21.74 35.38 , and 3.50 35.38 , or 15.7 : 12.9 : 61.4 : 9.9 respectively.

With the calculation of c ε from the experimental data, the next step is to create regression models to fit polynomial functions over a set of discrete variables and predict the light-use efficiency. As the aim of this study is to establish a framework that can estimate light-use efficiency as a function of incoming spectra, the input features include discrete parameters associated with either PFD or spectral irradiance distribution. The first set of input variables consists of photon flux density ratios for blue, green, red, and far-red wavebands, which are calculated by integrating the photon flux density distribution shown in Figure 1 . On the other hand, the second set of input parameters includes spectral irradiance ratios for the same wavebands, obtained by integrating over spectral irradiance distribution for various lighting distributions. Before investigating various regression models, the properties of the data used for regression are explored. Table 3 represents the mean, standard deviation (std), minimum (min), maximum (max), and percentile values for the investigated features (25% or first quartile, 50% or second quartile or median, 75% or third quartile), the input (F is the fraction ratio which is either based on photon flux density (PFD) or spectral irradiance integrated over wavelengths (I) whereas B, G, R, and FR represent blue, green, red, and far-red wavebands) and the output (c   ε is the light-use efficiency) of the model. Figure 6 visualizes the 3D distribution for photon flux density and spectral irradiance fraction ratios for the studied dataset.

Different polynomial models examined within the aim of this study have a form similar to

Not all of the models have every term presented in Equation (11), e.g., the regression model based on the first-order term is defined as c ¯ ε = a 0 + ∑ i = 1 4 ( a i F i ) . Regression models are defined in a way that includes up to 16 weight coefficients. Within the scope of this study, 22 distinctive terms are investigated, that are provided in Equation (11), and includes ratios ( F i , 4 terms representing each waveband), the square of ratios ( F i 2 , 4 terms), the cubic of ratios ( F i 3 , 4 terms), the quartic of ratios ( F i 4 , 4 terms), and interaction ratios ( F i F j , 6 terms). Therefore, studied regression models within the scope of this study are comprised of linear models with 4 (includes F i terms), 8, 12, 14, and 16 (includes F i , F i 2 , F i 3 , F i 4 terms) weight coefficients. Initially, a regression model with four terms corresponding to the first-order terms was developed and investigated as a baseline model. Higher-order polynomials were then explored to improve the poor performance (R   2 < 0.45%) of this linear model. Polynomials with 8 (combinations of linear, second-order, and interaction terms), 12 (same with third-order terms), 14 (same as 12 with two more terms), and 16 (including fourth-order terms) were considered. It should be noted that all 22 terms (all combinations through fourth-order) were not used because of the small size of the dataset. Results from these model choices will be given in the following section.

Since the dataset is comprised of 20 observations of c ε for different spectral distributions and intensities ( Figure 1 ), increasing the number of coefficients beyond the suggested limit would result in an overfitted model. In other words, this would lead to a model capable of accurate prediction for the studied data; however, evaluating the performance of the model against new data would decrease the fidelity of the model. The term “nonlinear” in this section refers to polynomial models based on Equation (11) in which e i j ≠ 0 for every i and j value. The regression models which estimate the impact of incoming light spectrum on the light-use efficiency are classified into the following categories: 1) models based on PFD ratios that disregard the impact of overall photon flux density; 2) models based onintegrated spectral irradiance ratios that disregard the impact of overall value; 3) models based on PFD ratios, that considered the impact of overall photon flux density; 4) models based on integrated spectral irradiance ratios that considered the impact of overall value. Figures 7A–D , represent the accuracy of categories 1, 2, 3, and 4 using R   2 metric.

To prevent overfitting and ensure the model’s applicability to new data, a validation study is conducted on the dataset. This involves reserving a portion of the data for testing, which is not used during the model training phase. The testing data is used to evaluate the model’s performance on new and unseen data. The goal is to find a model that performs well on both training data and testing data, thereby preventing underfitting or overfitting issues. For instance, for a regression model with 16 weight coefficients, 3 samples are randomly chosen for testing purposes. The remaining 17 samples are used to train the regression model, and then its accuracy is evaluated on the 3 unseen samples. The selected regression model (and the corresponding weights for different terms) is the one that performs well on both training data (17 samples) and testing data (3 samples). This approach ensures that the model is not overfitting and can predict well on new data, making it useful for practical applications. In addition to using unseen data for the validation of the regression model, a regularization penalty (L1 or L2 norm) is introduced into the regression model to decrease the variation caused by the complexity of the model.

In this study, the Scikit-learn built-in function LinearRegression developed by Pedregosa et al. (2011) is used to build a model that minimizes the regularized residual sum of squares (R   2 ) as the criteria for the closest linear functionto the actual data ¹ . The R   2 score is calculated using R   2 = 1 − ∑ i = 1 n ( y i − y ^ i ) 2 ∑ i = 1 n ( y i − y ) 2 , in which y i and y ^ i represented the true and predicted value for the i   t h sample and y is the average of the actual values. In addition to R 2 criteria, the mean absolute percentage error (MAPE) is also calculated for different models and is defined as MAPE = 1 n s a m p l e s ∑ i = 1 n s a m p l e s | y i − y ^ i | m a x   ( ϵ , | y i | ) , where ϵ is an arbitrary non-zero small positive number to ensure that MAPE is defined. Figure 8 compares the actual value, c   ε , s i m u l a t i o n , versus the predicted value, c   ε , p r e d i c t i o n : a straight line would indicate a perfect prediction. In Figure 8A , only terms linear in the features F i are included; in contrast, in Figure 8B , first-order and interaction terms in Equation (11) are used. It is seen that the accuracy increases with more features. Figure 8 provide qualitative insight into the accuracy, whereas R   2 and MAPE provide quantitative metrics, as shown in Figures 7 , 9 . According to Figures 7 , 9 , in spite of simplicity for the first-order and the combination of first and second-order models, these models have poor accuracy in predicting the c ε as a function of spectral distribution with R   2 score less than 0.75 . This is addressed by considering several models created from the general expression in Equation (11): ten different combinations of features are shown, each as a horizontal bar. Not all combinations are shown in the Figures 7 , 9 – 11 ; a representative set for a different number of terms is visualized, but one that includes the best performing model (16 terms).

Normalization of c   ε to accommodate the impact of overall integrated spectral irradiance or photon flux density was examined and found to increase the fidelity of the regression model. Note the improvement for the particular case of a model that includes first-, second-, third-, and fourth-order terms (Fi+Se+Th+Fo bar in Figure 7 ), which we will establish as the best model below. For this case, using normalization resulted in a 16-feature regression model that has an average error of 2%, as compared to the un-normalized error of 15%.

In addition to R   2 and MAPE metrics, the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) are also computed for the studied models, which are the measures of the model’s complexity. The value of AIC and BIC for models with a constant parameter is obtained using the following correlations (Seabold and Perktold, 2010): AIC = − 2 × L L + log   ( N ) × ( k + 1 ) , and BIC = − 2 × L L + 2 × ( k + 1 ) , whereas L L is the log of the likelihood function (how likely it is that the model predicted the actual values), N is the size of training data, and k is the number of features, respectively. It is worth mentioning that a lower value of AIC or BIC indicates that the fitted model is a better fit for the data, as it balances the goodness of fit and the complexity of the model model (Baguley, 2018). Figures 10 , 11 provide visual representations of the impact of the model’s complexity on AIC and BIC value. For the studied dataset, by increasing the accuracy of the model through the addition of new terms, the likelihood function also improves, which ultimately outweighs the negative penalty associated with the higher number of features (k); therefore, for the investigated models, the introduction of new features into the regression model would generally decrease AIC value. Figure 11 demonstrates that the BIC shows a different pattern, where regression models including nonlinear terms and only one of the first-, second-, or third-order terms are better fitted than those consisting ofnonlinear terms and two of the first-, second-, or third-order terms.

It is important to note that while a lower AIC or BIC value suggests a better model, it does not necessarily mean that the model is the best possible fit for the data. In order to determine the best-fitted model, a comprehensive analysis of the performance of the studied models is conducted based on AIC, BIC, MAPE, and R   2 metrics. Based on this comparison, it is evident that the model, which incorporates first-, second-, third-, and fourth-order terms based on the integrated spectral irradiance ratio considering the impact of overall integrated spectral irradiance of light-use efficiency, performs better than the other models with higher accuracy, and relatively lower AIC and BIC values.

Equation (12) demonstrates the high fidelity regression model based on spectral irradiance distribution normalized with integrated spectral irradiance for R180 light treatment based on Table 2 . As shown in Figure 12 , predictions of the suggested model are in good accordance with numerical data (for the best performing model, the testing data include treatment numbers 7, 14, and 19.); however, additional data could prove useful in developing a more comprehensive model.

Equation (12) is constrained by the following condition,

Now that the best-performing model is selected, the impact of adding higher-order terms to the regression model is investigated. The following notation [ a 1   ,   a 2   ,   a 3   ,   a 4 ] is used to evaluate the impact of higher-order term addition on first-order terms, in which a 1 , a 2 , a 3 , and a 4 represent weights for blue, green, red, far-red fraction ratios,respectively. The examined models included a first-order term model with weights of [ − 1.24 × 10 − 5   ,   − 6.96 × 10 − 6   ,   − 3.60 × 10 − 6   ,   − 2.18 × 10 − 7 ] , a combination of first- and second-order terms model with [ 6.95 × 10 − 5   ,   7.34 × 10 − 5   ,   8.02 × 10 − 5   ,   4.29 × 10 − 5 ] , a combination of first- to third-order terms model with [ 6.89 × 10 − 6   ,   − 2.41 × 10 − 5   ,   4.60 × 10 − 5   ,   − 1.60 × 10 − 5 ] , and a combination of first- to fourth-order terms model with [ 1.82 × 10 − 4   ,   2.32 × 10 − 4   ,   − 1.90 × 10 − 4   ,   4.68 × 10 − 4 ] . Comparison between these notations suggests that the addition of higher-order terms to the regression models significantly affects the weights of the first-order terms, highlighting the importance of exploring more complex models with higher-order terms for improved prediction accuracy.

Using Equation (12) and integrated spectral irradiance fraction ratio for the investigated wavebands of the different light treatment cases, the Van Henten (1994) dynamic growth model is utilized to compare the accuracy of the modified growth model that accounts for spectral distribution and intensity. Figure 13 compares the difference between experimental dry mass for lettuce with dry mass prediction of the modified growth model.

Comparing the numerical error of the lettuce growth model using the suggested value for light-use efficiency (labeled constant in Figure 13 ), with that of the lettuce growth model withthe proposed regression model (Equation (11); labeled regression in Figure 13 ), it is inferred that the suggested regression modelof the c ε improves the accuracy of Van Henten (1994) growth model and adequately considers the impact of spectral distribution on plant growth.

In the previous section, the performance of the proposed mathematical model is investigated through integration with the Van Henten (1994) dynamic growth model. To further investigate the accuracy of the light-use efficiency model based on the spectral irradiance ratios of various wavebands, data from Both et al. (1994) on lettuce cultivated in a controlled greenhouse environment is utilized. From October 1992 to March 1993, six controlled light treatments (without the use of supplemental lighting) were conducted for lettuce grown hydroponically. For these treatments, during the first 11 days, the temperature and CO   2 concentration were maintained at 25 ° C and 350 m L L ^-1, respectively. After day 11, the temperature was set to 24 ° C between 7 am and 5 pm and 18.8 ° C for the rest of the day, while CO   2 was enriched to 1000 mL L^-1 . Table 4 represents data used for cross-validation of the suggested regression-based light-use efficiency model. Spectral irradiance intensity and distribution ratios are approximated using reported daily light integral for the different treatments and predicted solar irradiation by Tobiska et al. (2000). Figure 14 displays a comparison of the spectral photon flux density distribution with respect to wavelength for the investigated LED spectrums and natural lighting.

According to Figure 15 , using the proposed regression-based light-use efficiency model, lettuce dry mass predictions were in close agreement (mostly within 2-3 percent error) with the experimental data.