Modeling leaf area development in soybean (Glycine max L.) based on the branch growth and leaf elongation

Several models have been proposed to simulate the leaf area index (LAI) in soybean (Glycine max L.); however, these models do not directly account for the effect of branch growth. Because the increases in branches and branch node vary with plant density, the evaluation of branch growth is necessary for the application of the LAI model at various plant densities. In this study, we developed an LAI model for soybean, considering the branch growth and leaf elongation at each node. To simplify this model, we estimated the rate of branch and branch node increase based on the rate of main stem node increase. Branch growth was assumed to be restricted when the fraction of canopy radiation interception was increased. Moreover, we calculated the leaf area growth at each node based on leaf elongation at each leaflet. This LAI model was validated using the data of different years and plant densities for model calibration, and it estimated the LAI with a root mean square error of 0.76, which accounted for 92% of the variation in the data. Although further evaluation is needed, the LAI model proposed in this study reveals a high potential for accurate estimation of LAI at various plant densities. ARTICLE HISTORY Received 11 March 2019 Revised 21 August 2019 Accepted 23 November 2019


Introduction
Solar radiation is a major determinant for biomass production. Therefore, for simulating the dynamics of crop growth and yield in the crop models, the accurate estimation of canopy-intercepted radiation is crucial. Canopyintercepted radiation comprises incident solar radiation or photosynthetically active radiation and a fraction of canopy radiation interception (FRI). The FRI is estimated using the leaf area index (LAI) in several crop models (Brisson et al., 2003;Robertson et al., 2002;Sinclair, 1986). At low LAI, FRI changes rapidly with an increase in the LAI values; therefore, the accurate estimation of leaf area development is essential during the early stages of crop growth (Sameshima, 2000;Sinclair, 1984).
Several crop models have been developed for soybean, e.g. SOYMOD (Meyer, Curry, Streeter & Mederski, 1979), Sinclair model (Sinclair, 1986), APSIM (Robertson et al., 2002), CROPGRO (Jones et al., 2003), STICS (Brisson et al., 2003), and SoySim (Setiyono et al., 2010), including ones for the Japanese soybean cultivars (Abe & Momoya, 1987;Sameshima, 2000). The Sinclair, STICS, SoySim, and Sameshima models evaluate the LAI development using the sink-limited approach, which assumes the leaf area development to be limited by temperature, whereas the APSIM and CROPGRO models evaluate the LAI by combining the sink-and source-limited approaches, which assume the leaf area development to be limited by the availability of assimilates. However, all the above models determine LAI based on the increase of the leaf number and leaf area on the main stem only and do not consider the leaf area development in branches directly. One of the reasons is presumably because many models target soybean cultivation in the USA where they cultivate at a high plant density. At high plant density, soybean plants produce considerably fewer branches than at low densities (Tenorio et al., 2017); however, to adapt these crop models to various plant densities, the leaf area development in branches must be considered (Cera et al., 2017). In Japan, soybean is cultivated at half or less plant density than that in the USA (Shiraiwa et al., 2011a). Consequently, the soybean plants produce more branches in the Japanese fields than in the US fields and the leaf area development in branches has a great impact on the LAI in the Japanese soybean fields (Reddy, Timlin & Pachepsky, 1999). Detailed studies on soybean branch growth in Japan indicate that it is a major factor affecting the increase in the soybean crop growth and yield at low plant density (Oizumi, 1962;Torigoe, Shinji & Kurihara, 1981).
In this study, we aimed to develop a leaf area model for soybean to predict the potential performance in terms of leaf area development before seed filling while considering branch growth and to adapt the model to low plant densities. To address these aims, we first estimated the main stem node increase as a function of temperature and the branch and branch node increases based on their coemergence relationship with the main stem node increase (Oizumi, 1962;Torigoe et al., 1981). Second, we estimated the leaf area growth at each node as a function of temperature. Lastly, we developed a leaf area model for the entire plant by combining the estimations of node increase and leaf area growth at each node.

Field experiments and measurements
Field experiments were conducted in the experimental upland field converted from paddy field of the NARO, Tsukubamirai (36°00´N, 140°01´E, Gray Lowland soil), Ibaraki, Japan from 2013 to 2015. In 2013, basal chemical fertilizer (2.5, 10.0, and 10.0 g m −2 of N, P 2 O 5 , and K 2 O, respectively) was applied to the soil before sowing. In 2014 and 2015, basal chemical fertilizer (2.7 and 3.3 g m −2 of N and P 2 O 5 , respectively) was applied to the soil before sowing; K 2 O was not applied based on the soil analysis results. Three determinate Japanese soybean cultivars, Ryuho, Enrei, and Fukuyutaka were planted in June and July each year in 0.7 m rows and 0.15 m interrows (Table 1). To increase plant density, 0.35 m row plots were added at the second sowing date in 2015. The experiments were conducted using the split-split plot design with three replicates where planting dates were the main plots, cultivars were subplots, and plant densities were sub-sub plots. In 2013 and 2014, each plot comprised six, 6 m long rows per 25.2 m 2 area. In 2015, each 0.7 m row plot comprised six, 7 m long rows per 29.4 m 2 area, and each 0.35 m row plot comprised eight, 7 m long rows per 19.6 m 2 area. Intertillage and ridging were conducted only in 0.7 m row plots before flowering, and no irrigation was applied.
A total of 8 plants from the central 4 rows in 0.7 m row plots and 12 plants from the central 6 rows in 0.35 m row plots were destructively sampled during the vegetative, flowering (R2), and initial seed filling (R5) stages. These growth stages were determined according to the method described by Fehr and Caviness (1977). All plants were manually harvested above the cotyledon node. Four representative plants were selected to determine the number of main stem nodes, branches, and branch nodes and to measure the dry weight of each organ (stem, leaf, and petiole) to determine the ratio of the dry weight of each organ. In this study, the unifoliate node was considered as the first node on the main stem, and the node with an unrolled leaf was also counted, as described previously (Tenorio et al., 2017). To determine the branch number, branches with more than two nodes were counted (Sayama et al., 2010). The total dry weight of the whole soybean plant was measured for all samples collected, and the dry weight of each organ was calculated by multiplying the whole plant dry weight with the ratio of the dry weight of each organ. Leaf area was measured from two plants per plot in 2013 and 2014 and one plant per plot in 2015 using the AAC-400 automation area meter (Hayashi Denko, Tokyo, Japan), and leaf area per plot was estimated from the leaf dry weight and specific leaf area. Additionally, images of leaves sampled from 0.7 m row plots in 2015 at the R5 stage were scanned (MFC-J6573CDW, Brother, Nagoya, Japan), and the length and area of each leaflet were measured using image-analysis software (LIA for Win32, https://www.agr.nagoya-u.ac.jp/~shinkan/LIA32). FRI was measured by digital image analysis (Shiraiwa, Kawasaki & Homma, 2011b). Canopy images were periodically recorded using a digital camera (COOLPIX P500, Nikon, Tokyo, Japan), approximately 1.8 m above the ground level. Images of two and three central rows in 0.7 and 0.35 m row plots, respectively, were selected, and FRI was measured from the fraction of the canopy cover using ImageJ software 1.48 h (Abramoff, Magalhaes & Ram, 2004). The upper limit of FRI was defined as 0.9 (Sameshima, 2000). The daily value of FRI after emergence was interpolated using a cubic spline from the measured data.
To determine the rate of individual leaflet elongation, experiments were conducted in the experimental upland field of the NARO, Tsukuba (36°01´N, 140°06´E, low-humic Andosol), Ibaraki, Japan in 2015. Basal chemical fertilizer (2.7 and 3.3 g m −2 of N and P 2 O 5 , respectively) was applied to the soil before sowing; K 2 O was not applied based on the soil analysis results. Two soybean cultivars, Enrei and Fukuyutaka, were planted on four different dates in 2015 (May 14, June 11, July 14, and August 20) in 0.7 m rows and 0.15 m inter-rows. The experiments were conducted using the split plot design with two replicates where planting dates were the main plots and cultivars were subplots. Each plot comprised six, 4 m long rows per 16.8 m 2 area. One representative plant was selected from each plot, and the length of the central leaflet of each trifoliolate leaf on the main stem node was measured periodically before R2. Actively growing leaves were defined as leaves with at least a 10 mm increase in the leaflet length at the successive measurement. Intertillage and ridging were conducted before flowering, and no irrigation was applied.
Daily weather data were obtained from a standard meteorological station located approximately 300 m away from the paddy and upland fields. Statistical analyses were performed using R ver.3.5.1 (R Core Team, 2018).

Model development
Although the main stem node increase is unaffected by plant density, branch and branch node increase is suppressed with increasing plant density (Oizumi, 1962;Torigoe et al., 1981). Therefore, the growth of the main stem node and branch or branch node was separately estimated in this study by dividing the leaf area development into three parts: main stem node increase, branch and branch node increase, and leaf elongation at each node. Moreover, LAI estimates affect the prediction of FRI largely when the LAI is small (Sameshima, 2000;Sinclair, 1984) and the LAI development during seed filling depends on the residue of assimilate available after seed growth (Robertson et al., 2002). Therefore, in this study, our model predicted the leaf area development before seed filling. A schematic illustration of the model concept for describing the node number increase and the replacement of elongating leaflet on each node is shown in Figure 1 and the model equations are presented below.
Description of the node number estimation model Because the increase of the main stem node is highly correlated with air temperature (Tenorio et al., 2017), the rate of main stem node increase was calculated using a three-segment function (Soltani & Sinclair, 2012), as presented below: where, MNN is the main stem node number (plant −1 ); TMN is the trifoliate leaf node number on the main stem (plant −1 ); int () is a function to round down to the nearest integer; RMN is the rate of main stem node increase (node d −1 ); a is a parameter; f () is a temperature function for the rate of main stem node increase; T is the daily mean temperature (ºC); TN b , TN p1 , TN p2 , and TN c are cardinal temperatures of base, lower optimum, upper optimum, and ceiling temperature for the main stem node increase (ºC); and the subscript i is the ith day after emergence. Since the determinate cultivars terminate the stem vegetative growth upon flowering or soon afterward (Bernard, 1972), parameter a was individually determined before and after the R2 stage. The values of TN b , TN p1 , TN p2 , and TN c were set at 10.0°C, 30.0°C, 35.0°C , and 45.0°C, respectively, based on previous reports (Sinclair, 1986;Soltani & Sinclair, 2012).
There is a strong coemergence relationship between the main stem node increase and primary branch growth of branch emergence and branch node increase (Oizumi, 1962;Torigoe et al., 1981). When the fifth trifoliolate leaf is emerging on the main stem, a primary branch emerges from the fourth lower node on the main stem (first trifoliolate leaf node). Therefore, the emergence of primary branches may be expressed as a function of the rate of main stem node increase. Moreover, the number of branches and branch node decreases with the increasing plant density due to differences in the light availability per plant (Agudamu, Yoshihira & Shiraiwa, 2016;Reddy et al., 1999;Torigoe et al., 1981). In this study, the difference in the light environment within the canopy was expressed as the difference in the FRI values, and the effects of FRI started to accumulate upon the formation of branch primordia. Primordia of the next emerging primary branches were assumed to be initiated when the main stem node, with a coemergence relationship with the primary branches, appeared. Therefore, the emergence of the primary branches on the main stem and the increase in the branch node on this branch were calculated as follows: where RBR is the growth rate of the primordium of primary branch (d −1 ); g () is a function of FRI for RBR; k is the radiation extinction coefficient; b is a parameter; IPB is the integrated values of the growth rate of primary branch primordium; EPN is the node number on each primary : Trifoliate leaf node on the main stem with an elongating leaflet : Trifoliate leaf node on the main stem withaout an elongating leaflet : Trifoliate leaf node on primary branches with an elongating leaflet : Trifoliate leaf node on primary branches without an elongating leaflet Figure 1. A schematic illustration of the model concept of the increase in the main stem node, primary branches, and primary branch node with the replacement of elongating leaflet on each node.
The effect of canopy shading on the branch growth has not been considered to illustrate in Figure 1. MNN indicates the main stem node number.
branch (branch −1 ); RBN is the rate of branch node increase (branch −1 d −1 ); and RMNv is the rate of main stem node increase using the parameter a before the R2 stage (node d −1 ). The subscript i represents the ith day after emergence and the subscript n represents the branch order of each primary branch from bottom to top on the main stem. In this study, the function g () was assumed as a linearly decreasing function of FRI, as sufficient data was unavailable to analyze the function. The value of k was set to 0.6, and the upper limit of FRI was set to 0.9 for the Japanese soybean cultivars (Sameshima, 2000;Shiraiwa, Hashikawa, Taka & Sakai, 1994). Moreover, when the value of IPB reaches one, a primary branch emerges on the main stem node. However, if the coemerging relationship fails because the growth rate of the primary branch primordia decreases, the branch emergence is assumed to have failed in the model.
Since only the branches with more than two nodes contributed toward the field measurement of branch number, the number of the primary branches with more than two nodes per branch was estimated from the results of branch node number estimation. Therefore, the numbers of primary branches and branch node for plant were calculated as follows: where CPB is the counter for courting primary branch number; PBR is the primary branch number for plant (plant −1 ); and PBN is the primary branch node number for plant (plant −1 ). Secondary branches also reveal a strong coemergence relationship with the growth of the primary branch node. However, the branch growth rate tends to be delayed, and the coemergence relationship is broken for late-emerging branches (Oizumi, 1962). Therefore, we calculated the secondary branch growth only from the first node on each primary branch to simplify the model. Similar to primary branches, when the sixth trifoliolate leaf was emerging on each primary branch, a secondary branch emerged from the fifth lower node on the branch, and the secondary branch node increased, followed by the coemergence relationship with the primary branch node.
Then, the numbers of branches and branch nodes for plant were calculated by integrating the primary and secondary branches as follows: where WBR is the branch number for plant (plant −1 ); WBN is the branch node number for plant (plant −1 ); SBR is the secondary branch number for plant (plant −1 ); SBN is the branch node number of secondary branches for plant (plant −1 ).

Description of the individual leaflet elongation model
Since the node number estimation is the same as the leaf number estimation in soybean, the leaf area for the whole plant may be calculated by adding estimates of the leaf area at each node. Leaf area of each leaflet may be estimated from the leaflet length because the length and width of a leaflet increase symmetrically (Inaba, 1992), while soybean leaflet shape is genetically controlled (Sawada, 1992;Sayama et al., 2017). Individual leaflet elongation rate is highly correlated with the air temperature (Bunce, 1985). Moreover, fully expanded leaf area increases toward the higher position of the main stem and branches (Cao, Hesketh, Zur & Reid, 1988;Saitoh, Nishimura & Kuroda, 2004;Sameshima, 2000;Sawada, 1992). In addition, the variation of fully expanded leaf area results not from leaf expansion duration but from the leaf expansion rate (Cao et al., 1988). Therefore, the leaflet area of each trifoliolate leaf was calculated including the modification of the maximum leaflet length and leaflet elongation rate with the nodal positions. Moreover, the elongation of each leaflet occurred simultaneously at the upper two nodes, and it stops when the leaf position shifted to the third node on the main stem and branches (Sameshima, 2000). Based on these studies, the equations for calculating the trifoliolate leaf area on the main stem are presented as follows: where FRL1 and FRL2 are reduction factors for maximum leaflet length and leaflet elongation rate for first and second upper leaflet, respectively; RIL1 and RIL2 are the individual leaflet elongation rates for first and second upper leaflet (mm d −1 ), respectively; h () is a temperature function for individual leaflet elongation rate, which is assumed to be the same as f () in Equation  4 according to the response of leaf area growth to temperature (Bunce, 1985;Magalhaes, Peters & Hageman, 1976); MLL is a parameter indicating the maximum individual leaflet length (mm); IIL is a parameter indicating the initial individual leaflet length when the new node emerges (mm); ILL1 and ILL2 are the individual leaflet lengths on first and second upper nodes (mm), respectively; TLA1 and TLA2 are the trifoliolate leaf areas on first and second upper node (mm 2 ), respectively; TLA3 is the total of below third upper node (mm 2 plant −1 ); MLA is the total of trifoliolate leaf area on main stem for plant (mm 2 plant −1 ); c, d, and e are parameters. We conducted the node number measurements based on the timing of the soybean growth stage; therefore, the node number increment as estimated by the model did not match the timing of leaf initiation. Consequently, we calculated the value of IIL based on the relationship between the first and second upper leaflet lengths, with the assumption that the node number estimated by the model increased mid-way through each node number duration. Because each trifoliolate soybean leaf comprised three leaflets, the trifoliolate leaf area at each node was determined as 3-fold greater than the leaflet area. Moreover, based on the measurements, the unifoliate leaves were assumed to appear as two pairs of 70.0 mm long leaflets when the number of the main stem node reached one and the unifoliate leaflet area is calculated in a similar manner as the trifoliolate leaf using the parameter e. The trifoliolate leaf area on each branch of the primary and secondary branches was also calculated in a similar manner as the leaf on the main stem, and the plant leaf area was calculated by integrating the unifoliate leaf area and trifoliolate leaf area on the main stem, primary branches, and secondary branches. Then, LAI was determined by multiplying the whole plant leaf area with plant density.
Although the process of leaf area development was modeled, the estimation of LAI also required determining the decrease in the leaf area due to leaf senescence. Given sufficient soil water and nutrient content, increase in the soybean leaf senescence responsible for the decrease in leaf area primarily occurs due to a reduction in irradiance in the lower canopy layer (Pons & Pearcy, 1994) and nitrogen translocation from leaves to seeds during seed development (Sinclair & de Wit, 1976). The decrease in LAI due to the leaf senescence by canopy shading was estimated using an empirical function (Goudriaan & van Laar, 1994) as presented below because sufficient data were unavailable to prepare a mathematical model of leaf senescence: where, RSL is the senescence rate of LAI (d −1 ), f and g are the parameters. The estimations of the leaf senescence at the main stem and branches were estimated simply by dividing the senescence LAI according to the ratio of LAI at the main stem and branches. Conversely, the decrease in LAI due to nitrogen translocation from the leaves to seeds was not considered in this study because the leaf area development was predicted before seed filling.

Model calibration and validation
Paddy field data from 0.7 m row plots in 2014 and 2015 were used to calibrate the parameters of node estimation and leaf senescence models, and upland field data were used to calibrate the parameters of individual leaf elongation model. To validate these models, results of the model were compared with the paddy field data collected from 0.7 m row plots in 2013 and 0.35 m row plots in 2015. Model performance was evaluated based on root mean square error (RMSE) and the normalized RMSE (NRMSE, %), calculated as the RMSE divided by the mean observed value. The present model includes nine parameters: two for the node number estimation model (a and b), five for the individual leaflet elongation model (c, d, MLL, IIL, and e), and two for the leaf senescence model (f and g). Parameter estimations were separately conducted for each model and parameter optimization was performed using Microsoft Excel Solver (Microsoft Excel 2016, Microsoft, Redmond, WA, US).
First, the parameters for the node estimation model were determined from the paddy field data. To avoid inclusion of the estimation error for phenology and FRI, we used the dates for the R2 stage through direct observation and direct measurements of FRI recorded via digital image analysis. The parameter a was evaluated to minimize the error between the estimated and measured main stem node number. The parameter b for branch and branch node number estimation was evaluated only from the branch node number estimation. This is because the branch node number is equal to the leaf number of the branches, the results of the branch node number estimation directly affected the leaf area development more than those of the branch number estimation.
Second, the parameters for the individual leaflet elongation model were determined from the measurements of leaflet length in the upland field. The values of c and MLL were graphically estimated from the vertical distribution of leaflet length on the main stem ( Figure 2). The value of IIL was evaluated from the relationship between the first and second upper leaflet lengths (Figure 3). For estimating individual leaflet elongation rate, the parameter d was evaluated to minimize the error between the estimated and measured leaflet lengths for an actively growing leaflet on the main stem. In this procedure, the parameters of c and MLL were used from the values determined above, but the values of IIL for each leaflet were determined from the measured values. The parameter e was evaluated from the relationship between the leaflet area and leaflet length at the R5 stages from the paddy field data of 2015 using linear regression.
Finally, the parameters f and g were evaluated to minimize the error between the estimated LAI using the node number estimation and individual leaflet elongation models and LAI measured in the paddy field.

Calibration of the model parameters
The model parameters estimated from the model calibration data are presented in Table 2. The number of main stem node, branches, and branch node, estimated using the measured FRI values, reproduced the calibration data (Table 3). Regression analysis between the estimated and observed values yielded regression lines with slopes not significantly different from unity and with intercepts not significantly different from 0. The model accounted for 81 to 82% of the variation in the calibration data. The NRMSE for branch number was higher than that for branch node number because the parameter for the response of branch and branch node increase to FRI was optimized only for the estimation of branch node number.
The length of the fully expanded leaflet increased with increasing the position from the base to around the 10th trifoliolate leaf node on the main stem and decreased rapidly beyond this position from May to July (Figure 2). Because these measurements were conducted until R2, some leaflets on upper nodes were still elongating and did not reach full expansion. The fully expanded leaflet length in the August sowing was relatively stable across the nodal position, and the dependence pattern of the leaflet length on the nodal position differed from that observed in the other sowing periods. However, sowing in August is uncommon for soybean cultivation in Japan, we excluded the August sowing data when we parameterized the individual leaflet elongation model. The leaflet length of the first trifoliolate leaf was approximately one-third of the maximum leaflet length. To simplify the model, we defined the value of MLL as 160 mm and assumed that the fully expanded leaflet length decreased linearly from the fifth trifoliolate leaf node to one-third of the value at the first trifoliolate leaf node (c = 5 in Equations 16 and 17) excluding the differences among the cultivars and sowing dates. The model responses reproduced the relative change of the fully expanded leaflet length with the nodal position of the main stem. The first and second upper leaflet lengths showed a significant linear relationship (Figure 3). If the increase of a new node was defined to be when the first upper leaflet length reached 20 mm (Sinclair, 1984), the second upper leaflet length was calculated to be approximately 69 mm, indicating that the first upper leaflet length could vary within the range of 20-69 mm. Because our model did not estimate the node number at the time of leaf initiation, we set the initial value of ILL at 45 mm, the median value of the first upper leaflet length.
The parameter c, the potential rate of individual leaflet elongation, was defined as 25.58 mm d −1 . The slope of the regression line between the model-estimated and observed leaflet lengths of the actively growing leaf was close to 1, and the model (Equations 16-21) accounted for 87% of the variation in the calibration data (y = 0.94x + 4.08, R 2 = 0.87, p < 0.001). In addition, the relationship between the square value of the leaflet length and the leaflet area of the trifoliolate leaf showed a significant linear relationship (R 2 = 0.97; p < 0.001), and the parameter e was defined as 0.45.
Eventually, the parameters f and g in the leaf senescence model were defined as 0.044 and 5.14, respectively, by fitting the model-estimated LAI to the observed values of calibration data. The regression equation between model-estimated and observed LAI of the model calibration data was close to the 1:1 line and the model accounted for 89% of the variation in the calibration data (Table 3).

Validation of the leaf area model
We evaluated the accuracy of our leaf area model using the field data collected from 0.7 m row plots in 2013 and 0.35 m row plots in 2015; these data were not used for model construction. The slopes of the regression lines between the model-estimated and observed main stem node number, branch number, branch node number, and LAI were close to 1 ( Table 3). The model accounted  for 92-98% of the variation in the model validation data. The value of NRMSE was the lowest for the main stem node number (5.38%) and highest for the LAI (21.92%). Since the residual errors (estimated − observed) against the observed values did not show consistent trends for all variables, systematic bias was not detected for the validation data (Figure 4).
To evaluate the effect of row width, the model results were shown separately from 0.7 m and 0.35 m row plots. The slopes of the regression lines between modelestimated and observed were close to 1 excluding the branch number and LAI of 0.7 m row plots in 2013 (Table 3). However, the NRMSE values of branch number of 0.7 m row plots were lower than those of 0.35 m row plots. The value of RMSE for LAI estimated from 0.35 m row plots was higher than that from 0.7 m row plots; however, the NRMSE value of 0.35 m row plots was lower than that of 0.7 m row plots due to higher LAI estimate from 0.35 m row plots at a high plant density.
The seasonal patterns of the LAI at the main stem and branches differed depending on cultivars, sowing dates, and plant densities ( Figure 5). The LAI of Fukuyutaka was larger than that of Enrei at R2 and R5, particularly for early sowing (June sowing in 2013), because the growth duration of Fukuyutaka was longer than that of Enrei (Table 1). In addition, although the LAI at branches was similar or higher than that at the main stem in 0.7 m row plots at R2, it was smaller in 0.35 m row plots. Our model did not calibrate the LAI at the main stem and branches separately, but it estimated the overall trend of the LAI differences between main stem and branches for these conditions. Moreover, this model estimated not only the leaf area increase before R2 but also the leaf senescence after R2 on the main stem.

Discussion
In this study, we aimed to develop a mechanistic model to estimate LAI, considering the increase of branches and branch nodes and the expansion of each leaflet. In addition, to avoid complicating the model with additional equations and parameters, we simplified the model structure assuming the regularity of branch increase and the stability of leaf shape in soybean.
A few models estimated LAI directly from the number of main stem leaves using empirical allometric equations   (Robertson et al., 2002;Sinclair, 1984). However, these LAI models are difficult to apply for various plant densities because they do not account for branch growth (Cera et al., 2017). Compared with the previous models, our model is advantageous because it evaluates the number of branch and branch nodes (Table 3 and Figure 4). Introducing a function that indicates suppression of branch and branch node increase due to an increase in FRI, our model quantified the branch growth in response to the variation in canopy light environment. This response enabled the model to estimate LAI at the plots of narrower row width than the plots used to estimate the model parameters (Table 3). In addition, the estimation of branch growth facilitated the prediction of the LAI separately for the main stem and branches ( Figure 5). These results indicated that our model is applicable to a wide range of plant densities because it considers the branch and branch node increase. Although the regression line slope for LAI differed from unity in 2013, it was mainly caused by the underestimation of LAI above the value of 9 (Figure 4). However, this underestimation is not expected to have a large effect on the estimation of light interception because a canopy closure in Japanese soybean cultivars occurs at LAI > 4 (Sameshima, 2000). Furthermore, the branch and branch node number estimations for the validation data were better than those for the calibration data that used measured FRI values (Table 3). These results indicate that the proposed model describes the process of the leaf area development reasonably well by integrating the increases in branches and branch nodes with the expansion of each leaflet.
However, there were two empirical functions determined experimentally from periodic length measurements: the value of IIL and the distribution of the fully expanded leaflet length along with the nodal position (Figures 1 and 2). For a more theoretical approach to determine the value of IIL, the plastochron concept can be used to estimate the leaf appearance rate, although additional measurement of the leaflet length is needed at the node number measurement (Hofstra, Hesketh & Myhre, 1977). The distribution of the fully expanded leaflet length along with the nodal position is associated with the difference of the formation of leaf primordia. The first five leaves of maize and sorghum, which are preformed in the seed, remain small (Tardieu, Granier & Muller, 1999), and soybean also forms five leaf primordia (third trifoliolate leaf) at the cotyledon stage (Tenorio et al., 2017). However, the distribution pattern depends on the difference of growth habit of soybean (Saitoh et al., 2004). More data are needed to clarify the relationship between the leaflet area and nodal position. Finally, in this study, the drought and phosphorus stresses were not considered for leaf area development because this study was primarily conducted in paddy fields where soil moisture is relatively high and phosphorus was applied sufficiently as basal fertilizer. However, these stresses are known to decrease the leaf area development in soybean, especially due to decreasing branches and individual leaf area (Chiera, Thomas & Rufty, 2002;Gutierrez-Boem & Thomas, 2001;Randall & Sinclair, 1988). Therefore, as the next step, the model should be improved with incorporation of these effects on branch growth and leaf elongation.
Although further studies are needed, the model presented in this study is the first LAI model for soybean that accounts for branch development and predicts the potential leaf area development over a wide range of plant densities.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was partly supported by the Council for Science, Technology, and Innovation (CSTI), Cross-ministerial Strategic Innovation Promotion Program (SIP), 'Technologies for creating next-generation agriculture, forestry and fisheries' (funding agency: Bio-oriented Technology Research Advancement Institution, NARO).