- © 2018 American Society of Plant Biologists. All Rights Reserved.

## Abstract

Process-based crop growth models are popular tools with which to analyze and understand the impact of crop management, genotype-by-environment interactions, or climate change. The ability to predict leaf area development is critical to predict crop growth, particularly under conditions of limited resources. Here, we aimed at deciphering growth coordination rules between wheat (*Triticum aestivum*) plant organs (i.e. between leaves within a stem, between laminae and sheaths, and between the mainstem and axillary tillers) to model the dynamics of canopy development. We found a unique relationship between laminae area and leaf rank for the mainstem and its tillers, which was robust across a range of sowing dates and plant densities. Robust relationships between laminae and sheath areas also were found, highlighting the tight control of organ growth within and between phytomers. These relationships identified at the phytomer scale were used to develop a simulation model of leaf area dynamics at the canopy level that was integrated in the wheat model SiriusQuality. The model was then evaluated using several independent experiments. The model accurately predicts leaf area dynamics under different scenarios of nitrogen and water limitations. It accounted for 85%, 64%, and 73% of the variability of the surface area of leaf cohorts, total leaf area index, and total green area index, respectively. The process-based model of the dynamics of leaf area described here is a key element to quantify the value of candidate traits for use in plant breeding and to project the impact of climate change on wheat growth.

The development of leaf area and its structural components are major determinants of light, carbon, nitrogen, and water capture by plants. It also largely defines plant interactions with pests and diseases (Andrivon et al., 2013). Leaf expansive growth also is one of most sensitive processes to temperature (Parent et al., 2010) and to a wide range of abiotic stresses (Pantin et al., 2012). Leaf growth and leaf area index (LAI; the surface area of green leaf laminae per unit of ground area) is a key trait with which to understand and model plant responses to environmental changes (Ewert, 2004).

Modeling approaches used to model canopy development assumed either that leaf growth is temperature driven (i.e. sink limited; Hammer et al., 1993; Jamieson et al., 1998b), carbon limited (i.e. source limited), where daily leaf mass production is converted to leaf area (Stöckle et al., 2003; Yin and van Laar, 2005), or a combination of the two (Olesen et al., 2002) . Experimental evidence shows that expansion growth is primarily sink limited (Tardieu et al., 2014), and in a comparison of different wheat (*Triticum aestivum*) growth models, Jamieson et al. (1998a) found that sink-limited leaf area models were able to simulate the response of LAI to water and N limitations more accurately than source-limited models. The choice between these modeling approaches also has critical implications on the ability to link model parameters with genetic information (Parent and Tardieu, 2014).

The level at which the canopy is considered in crop growth models depends on the objectives of the models. The conventional approach has been to consider the canopy as homogenous and to use descriptive allometric relationships. Although this approach has been very effective in simulating biomass production, it has several limitations for modeling plant N economy and the effects of abiotic and biotic stresses (Birch et al., 2003) or in linking model parameters to genetic information (Chenu et al., 2008). Therefore, individual leaf (scale) models have been developed, especially for maize (*Zea mays*; Birch et al., 1998; Lizaso et al., 2003; Chenu et al., 2008) and wheat (Porter, 1984; Lawless et al., 2005). To overcome the complexity of the dynamics of tillers, wheat leaf cohort models (on mainstem and axillary tillers) have been developed (Porter, 1984; Lawless et al., 2005), but these models have only been evaluated for their ability to predict total LAI dynamics. Three-dimensional functional-structural plant models (FSPMs) also have been developed (Fournier et al., 2003; Evers et al., 2005; Verdenal et al., 2008). However, current FSPMs do not simulate several important processes, such as leaf senescence, carbon and nitrogen remobilization, or the effect of belowground resource availability. Moreover, FSPMs remain difficult to parametrize for new genotypes and to evaluate in field conditions, and most of them do not simulate the entire crop growth period.

In cereal crops, the development of the canopy takes place through the sequential production of leaves on the mainstem and its axillary tillers. There is a strong coordination between the mainstem and its tillers and between the different leaves depending on their positions along the stem, and leaf development is dependent on the development of the previous leaf, either on the same tiller or on the mainstem (Bos and Neuteboom, 1998; Tivet et al., 2001; Zhu et al., 2015). Each leaf is composed of a lamina and a sheath enclosing an internode, which with an axillary bud constitute a phytomer. The emergence of the collar (the junction of sheath and lamina) is coordinated with the development of the lamina and internode (Andrieu et al., 2006). Leaf growth also is tightly related to changes in the apical meristem and, thus, to flowering time (Brooking and Jamieson, 2002). Therefore, vegetative and reproductive developments are not independent, since floral initiation occurs during leaf development and is strongly related to the final number of leaves on the mainstem through a reduction of the duration of the phase of leaf primordia production (Brown et al., 2013) and has a large impact on final leaf size (Borrill, 1959). The precise coordination between the mainstem and its tillers and between successive leaves on an axis have to be considered to model leaf expansive growth at the whole-plant level.

In this study, we aimed at identifying environmentally robust relationships describing the coordination between wheat shoot organ sizes (i.e. between the mainstem and its tillers and also between laminae and sheaths). Using these simple relationships, we developed and evaluated a new leaf expansion model implemented in the wheat crop model SiriusQuality (Martre et al., 2006, 2015).

## RESULTS

### Coordination of Final Lamina Area between the Leaves of the Mainstem and Axillary Tillers and between Leaves along an Axis

We used two experiments (experiments A and B in Table I) to decipher the coordination between the final lamina area of the leaves of the mainstem and axillary tillers. In these experiments, the winter wheat cv Soissons was sown at different dates and densities in two locations in France. For both experiments, the final number of leaves on successive tillers decreased. However, the final area of the penultimate leaf was similar for the different tillers (Fig. 1). Moreover, the slope of the relationship between the final lamina area and the leaf number for the last six leaves was similar for the different tillers (Fig. 1, A–E) and experimental conditions (Fig. 1F). This suggests that a unique relationship may exist between the lamina area and leaf rank for all the tillers, whatever the sowing date or density.

When the leaves were counted acropetally (from the base of the plant), the relationship between the final lamina area of the mainstem and that of an axillary tiller depended on the tiller number, in particular for tiller order greater than 3 (Supplemental Fig. S1). To verify if a unique pattern of leaf area along the different axes can be identified, the lamina areas of the different leaves and tillers were normalized by the area of the penultimate leaf lamina of the mainstem and were plotted against the leaf rank counted basipetally (from the top to the base of the plant). Moreover, in order to account for the increasing delay between the emergence of successive tillers (Fournier et al., 2003), the phytomer numbers of axillary tillers were shifted by a decimal number (Evers et al., 2005). After shifting the second tillers along the *x* axis by +0.35 phytomer, the third, fourth, and fifth tillers by +0.75 phytomer, and the sixth and seventh tillers by +1 phytomer, a unique relationship was observed between the normalized lamina area and the leaf rank counted from the flag leaf (Fig. 2).

The normalized lamina area of the first mainstem leaves (in this article, the first leaf always refers to the first leaf counted from the base of the plant) was similar, whereas it increased almost linearly for the last six leaves (Fig. 2). The normalized lamina area of the flag leaf was much more variable. In wheat, leaf growth is coordinated with the expansion of the internodes (Giunta et al., 2001). This coordination is most likely trigged early in development by floral initiation and regulated by GA accumulation in the shoot apex (Pharis and King, 1985; King and Evans, 2003). Thus, the number of large leaves depends on the number of leaves initiated after the floral initiation. For winter wheat sown in autumn, the number of leaves initiated after floral initiation may vary between four and six (Brown et al., 2013), which is close to the intercept of the relationship between the normalized lamina area and the leaf rank counted basipetally (Fig. 2F).

The relationship between the normalized lamina area and the phytomer number counted basipetally was used to model the potential final lamina area (i.e. in the absence of biotic or abiotic limitations) for the different phytomer numbers and tillers (). We distinguished two classes of leaves, the first small (juvenile) leaves and the last four to six large leaves. We assumed that of the juvenile leaves is constant (; cm^{2}). Then, of the different phytomers and tillers was modeled as:where *i* is the phytomer number counted from the base of the plant, *j* is the tiller order, (cm^{2}) is the potential final lamina area of the mainstem penultimate leaf, (dimensionless) is the ratio of the mainstem flag-to-penultimate laminae area, *N* (leaf) is the final leaf number, η (leaf) is the *x* intercept of the relationship between and the phytomer number for the mainstem large leaves, and *s* is the number of phytomers by which the phytomer number of the different tiller order is shifted (for parameter definitions and units, see Supplemental Table S1).

### Coordination of Final Organ Sizes between the Lamina and Its Associated Sheath

In experiment A (Table I), sheath surface area was measured to analyze the coordination between the final area of laminae and sheaths. The area of the sheath increased with leaf rank (Fig. 3, A and B). After normalizing the sheath area by and shifting the phytomer number of the tillers as for the laminae, a unique pattern of sheath area was observed for the different tiller orders and the two sowing densities (Fig. 3E). As a consequence, a unique relationship was found between the areas of the lamina and sheath, except for the flag leaf, whose sheath and lamina did not follow the coordination rule observed for the other leaves (Supplemental Fig. S1).

To model light interception by leaf sheaths, one has to model the area of sheaths between successive leaf collars, which corresponds to the exposed area of the pseudostem (; cm^{2}). This requires modeling the lengths of the sheaths and internodes for the last leaves with expanded internodes. Here, we found that also increased with the phytomer number (Fig. 3, C and D), and a unique relationship was found between normalized by and *s*, as for the total sheath area (Fig. 3F). This relationship provides a direct way to model. In wheat, the beginning of stem extension occurs after floral initiation, approximately at the terminal spikelet stage of development (Hay, 1978; Hay and Kirby, 1991), and the final number of elongated internodes is closely related to the number of leaves emerged after the terminal spikelet stage (Giunta et al., 2001). This explains why the sheath area starts to increase approximately one phytomer after that of the laminae. The relationship in Figure 3F was modeled as:where (cm^{2}) is the potential final sheath area of the first leaf.

As summarized in Figure 4, Equations 1 and 2 allow calculating the final area of leaf laminae and light-exposed sheaths of the mainstem and axillary tillers with only five parameters. In general, the equations predicted well the final area of the mainstem and axillary tiller leaves. However, the area of mainstem leaves 4 to 6 was underestimated, because we assumed that mainstem juvenile leaves have a constant size.

### Prediction of LAI under Water- and Nitrogen-Unlimited and -Limited Conditions

Equations 1 and 2 were used to develop a model of the dynamics of leaf area. The timing of leaf expansion and the responses to water and nitrogen limitations were modeled using five and six additional parameters, respectively (see “Materials and Methods”; Supplemental Table S1). The model was implemented in the wheat model SiriusQuality. In SiriusQuality, the variations of the time to anthesis associated with vernalization requirement and daylength sensitivity are described in terms of primordium initiation, leaf production, and final mainstem leaf number (He et al., 2012). Thus, the model of leaf area described here links the fate of the mainstem apex with the size of the organs produced by the meristem.

The model was evaluated using three independent field data sets (experiments C–E in Table I). First, we used a data set where the winter wheat cv Caphorn was grown under nonlimiting soil nitrogen and water supply (experiment C in Table I). At the canopy level, the model simulated well the dynamics of LAI and green (laminae + light-exposed sheaths) area index (GAI; Fig. 5, A and B). The root mean squared error (RMSE) was lower for LAI than for GAI (Supplemental Table S2). For both variables, the lack of correlation was the highest component of the mean squared error (MSE). For GAI, the nonunity slope was the second highest component and the slope of the relationship between measured and simulated GAI was greater than 1 (Fig. 5B). The measurement error for GAI at the end of the growing period was high, which could explain the large contribution of the lack of correlation component to MSE. The model was able to predict reasonably well the LAI of each leaf cohort (Fig. 5, C and D; Supplemental Table S2; *r*^{2} = 0.87).

A second independent experiment was used to assess LAI predictions under different scenarios of water deficit (experiment D in Table I). The model captured the effects of early and late droughts of varying durations on the dynamics of LAI from crop emergence to maturity (Fig. 6, A–C). Overall, the model explained 85% of the observed variation of LAI when considering all the scenarios (Fig. 6D; Supplemental Table S2). For all treatments, the lack of correlation was the highest component of the MSE. This was largely due to the noise in the measured LAI, especially for the well-watered condition (Fig. 6A). For all the water deficit scenarios, the time when the predicted LAI was maximum was delayed compared with the measurements.

The last experiment (experiment E in Table I) used to assess the model prediction skills consisted of two N treatments (low and high) and two different sites (Clermont-Ferrand, France, and Sutton Bonington, United Kingdom) during two consecutive growing seasons (2006-2007 and 2007-2008). In Clermont-Ferrand in 2006-2007, a water deficit during the stem extension period led to a severe limitation of the maximum LAI (Fig. 7A). The cv Soissons was used in this experiment. The model was calibrated using the two data sets used for the development of the leaf area model (experiments A and B) to predict LAI and GAI for experiment E. Our model explained 67% and 69% of total LAI and GAI variability, respectively (Fig. 7; Supplemental Table S2). The best predictions were for LAI of the leaf cohorts (*r*^{2} = 0.84).

## DISCUSSION

In this study, we investigated the patterns of leaf lamina and sheath final areas in relation to their positions on the mainstem and on axillary tillers that emerge from coordination events between successive leaves. We found a robust relationship between the final area of leaf laminae and sheaths along an axis and between axes. Most of the previous studies identified relationships between organ length and width and model leaf area for wheat (Bos and Neuteboom, 1998; Fournier et al., 2005; Evers et al., 2007), barley (*Hordeum vulgare*; Kirby et al., 1982), rice (*Oryza sativa*; Tivet et al., 2001), or maize (Andrieu et al., 2006). Here, we showed that patterns similar to those reported for organ length and width can be found for their areas. This allowed us to develop a parsimonious and robust model of leaf and sheath area expansive growth and its response to water and nitrogen supply.

We confirmed previous observations, according to which the sizes of leaf laminae and sheaths on the mainstem and axillary tillers follow the same developmental patterns when the phytomer number is corrected to account for the difference in final leaf number of each axis. The reason for the similarity in mature leaf size between the mainstem and axillary tillers is likely the quasi-synchrony of floral transition between the different axes of a plant (Hay and Kirby, 1991). The concept of phytomer shift and relative phytomer number was introduced by Fournier et al. (2003) and is used in the architectural model ADEL-wheat (Evers et al., 2005). The value of the phytomer shift (*s*) is related to the final leaf number of the tillers, as the integer of the final leaf number of a tiller multiplied by *s* + 1 is equal to the final leaf number of the mainstem. This concept allows a simple description of final organ size along the mainstem and the axillary tillers, which are not modeled independently. As the number of juvenile leaves is variable in wheat (Dornbusch et al., 2011a), in our model, the number of phytomers is counted basipetally to reduce errors arising from the uncertainty regarding the number of juvenile leaves. Here, we used the area of the penultimate leaf of the mainstem to normalize the development pattern of leaf size. This is a convenient way to model the observed coordination of mature leaf size along the mainstem, and as mentioned above, it reduces the model error due to variations in the number of juvenile leaves. However, it is important to recognize that it is an emerging property of coordination rules in the sequence of succession of leaf emergence events, as shown for maize (Zhu et al., 2014).

To evaluate the overall deviation of the model simulations, we decomposed the MSE into its three components. The errors of the model were due mainly to a lack of correlation. At early developmental stages, errors could be explained by the overestimation of the leaf area of the juvenile leaves. Indeed, the final area of the first leaves has been shown to depend strongly on the temperature and photoperiod during their growth (Gallagher et al., 1979; Kirby et al., 1982). For wheat, it usually increases for spring sowings and decreases for winter sowings (Borrill, 1959; Dornbusch et al., 2011a). Under our growing conditions with winter sowing, the assumption of a constant potential final size of the juvenile leaves was consistent with our field measurements, and there were no significant carryover effects. For spring wheat sown in the spring, the floral transition occurs much earlier in the development than for winter wheat, and leaf size increases linearly from the first or second leaf (Evers et al., 2005; Moeller et al., 2014). In the current version of the model, this difference between spring and winter wheat is modeled via the genotypic parameter η. However, since vernalization saturation (i.e. when vernalization requirements are satisfied), which is modeled in SiriusQuality (He et al., 2012), occurs either at or a few days before floral initiation (Brown et al., 2013), the value of η can be predicted as the number of emerged leaves at vernalization saturation. Thus, the differences in the pattern of leaf size between spring and winter wheat can be modeled without the need of a parameter. However, changes in the potential final size of the juvenile leaves will require more attention to use our model to study the effect of early vigor traits on spring wheat performance, especially for winter sowings, where floral initiation may be delayed because of the short daylength.

For later developmental stages, model errors for LAI were explained mainly by the observed variations of the flag leaf area and of the ratio of flag-to-penultimate leaf area. During the period of flag leaf expansion, the plant demand for carbon is high, since at that time the growth rates of the stem and ear are maximum. Thus, the changes of the ratio of penultimate to flag leaf area might be due to carbon competition. Therefore, our model probably could be improved by considering the carbon supply-to-demand ratio during leaf growth, as was done for modeling tillering in rice (Luquet et al., 2006) or sorghum (*Sorghum bicolor*; Alam et al., 2014).

The model tended to underestimate LAI under water deficit conditions before anthesis for all watering scenarios. This was due mainly to a delay of LAI increase during early developmental stages. This delay shows the strong interplays between phenological development and LAI establishment in SiriusQuality. Although leaf expansion is very sensitive to water deficit, the duration of leaf expansion also might be sensitive to drought (Simane et al., 1993; Estrada-Campuzano et al., 2008). In the model, we considered the effect of water deficit on leaf expansion and leaf senescence. However, we did not consider any environmental effect on the phyllochronic duration of leaf expansion, except an independent effect related to the increase of the canopy temperature under drought.

We modeled the sheath area exposed to light rather than the total area of the sheaths. This simplification gives an accurate estimation of the sheath area to calculate light interception and carbon assimilation by the stems. However, this may lead to higher errors in modeling carbon and nitrogen contents of the sheaths (in the model, the carbon and nitrogen of the hidden sheath parts are allocated to the stem). The nitrogen mass per unit of sheath area is constant along individual sheaths, but it decreases strongly between phytomers from the top to the bottom of the canopy following the vertical light gradient (Bertheloot et al., 2012). Therefore, modeling the area of light-exposed sheath parts may result in an underestimation of plant nitrogen content, which can explain parts of the errors for LAI prediction under the different nitrogen supplies.

The model was able to simulate reasonably well the dynamics of LAI for fertile shoot numbers ranging from 390 (for cv Soissons at Clermont-Ferrand under low nitrogen) to 669 (for cv Soissons at Sutton Bonington under high nitrogen). Therefore, the simulation results support our modeling approach to model LAI.

Here, we did not attempt to model the growth of leaves before their emergence above the collar of the preceding leaf. This means that the effect of stresses during the early growth period of the leaf was not considered. The growth of the leaves before they emerge could be modeled by considering the coordination of events between successive leaves. In wheat, the blade of leaf *n* and the sheath of leaf *n* + 1 start to grow when the tip of leaf *n* + 1 appears (Fournier et al., 2005).

Several studies have focused on deciphering the coordination rules between organs to predict the balance between biomass accumulation and organ size (Song et al., 2015), the timing of plant growth (Fournier et al., 2005; Verdenal et al., 2008; Zhu et al., 2014), and the final size of organs (Evers et al., 2005). These relationships are parameterized mainly in FSPMs, which are descriptive models that need to be calibrated in each new condition. The individual leaf growth model developed here allows considering the effect of water and nitrogen limitations on dry matter accumulation and actual organ expansion at a daily time step. It can be used to analyze the response of the plant to the environment and to identify architectural trait combinations (i.e. ideotypes) to optimize light capture by the canopy and biomass production in targeted environments.

## MATERIALS AND METHODS

### Field Experiments Used to Analyze the Patterns of Final Leaf Laminae and Sheath Surface Area

We used two field experiments to decipher the patterns of final leaf laminae and sheath surface area. The first experiment (experiment A in Table I) was conducted with the winter bread wheat (*Triticum aestivum*) ‘Soissons’ at Grignon, France, on a deep loamy soil (Ljutovac, 2002; Fournier et al., 2005). The experiment consisted of two sowing density treatments, a normal density that yielded a canopy with 250 plants m^{−2} and a low-density treatment (70 plants m^{−2}) to allow for the development of additional tillers. The two sowing density treatments were conducted in adjacent areas as randomized block designs with three replicates. Subplot size was 30 × 8 m with an interrow spacing of 0.175 m. The second experiment (experiment B in Table I) was conducted at Villiers le Bâcle, France, where cv Soissons was sown with an optimal and a late sowing date. For the optimal sowing date, two sowing densities were used. The two sowing density treatments were conducted in adjacent areas as randomized block designs with three replicates. Subplot size was 10 × 1.5 m with an interrow spacing of 0.125 m.

In both experiments, crop inputs were applied at levels to prevent nutrients or pests, diseases, and weeds from limiting canopy expansion, and no growth regulator was used. In experiment A, a drip irrigation system was installed, but the crops were irrigated only once near the end of the measurement period. Sixty plants were collected every 2 to 3 d in each individual plot, and the median values of leaf stage and lamina length of the last ligulated mainstem leaf were calculated and used to identify 10 median plants, which were dissected to measure the lengths of individual laminae of the mainstem and each tiller. The areas of the individual leaf laminae were then calculated as:where *l* (cm) is the laminae length from the ligule to the tip, *w*_{max} (cm) is the maximum lamina width, and *f*_{t} (dimensionless) is an overall shape parameter. *f*_{t} varies slightly with leaf rank, and we used the values determined on an independent data set for cv Soissons (Dornbusch et al., 2011b). In experiment A, final lengths and maximum widths of the leaf sheaths also were measured and multiplied to obtain the final areas of the sheaths.

### Field Experiments Used for Model Evaluation

Three independent experiments were used to evaluate the robustness of the leaf area model. The first one (experiment C in Table I) was conducted at Clermont-Ferrand, France, with the winter wheat cv Caphorn. Seeds were sown on December 1, 2005, at a density of 300 seeds m^{−2} in a clay loam soil. Nitrogen fertilizer was applied as ammonium nitrate granules in four splits to give nonlimiting growth conditions. The data used here are from a larger experiment with two N treatments and four winter wheat cultivars, and a randomized complete block design was used in which N treatments were randomized on main plots, cultivars were randomized on the subplots, and each treatment was replicated three times. Subplots size was 25 × 2 m with an interrow spacing of 0.175 m. The crops were rainfed. Starting at the emergence of leaf 4 until maturity harvest, in each subplot three central rows were collected on 50 cm (approximately 60 plants per sample) every 7 d. The total fresh mass of the samples was determined, and a 25% subsample (by fresh mass) was selected. Thirty plants from the subsample were randomly selected and dissected into individual leaf laminae, stem (including leaf sheath), and ear. The projected surface areas of green laminae, stems, and ears were determined using a Li-3100 area meter (LI-COR). Actual stem and ear surface areas were calculated as their projected surface area multiplied by π/2 (Lang, 1991). Canopy GAI (m^{2} green tissue m^{−2} ground) was calculated as the sum of lamina, stem, and ear surface areas of the main and secondary shoots.

The second experimental data set (experiment D in Table I) used to evaluate the model was from a rainshelter experiment carried out at Lincoln, New Zealand, with the spring wheat cv Batten (Jamieson et al., 1995b). The mobile rainshelter (55 × 12 m) automatically covered the experimental crop during rainfall but was otherwise positioned 50 m away. The working section of the rainshelter was divided into 24 3.6- × 5-m subplots. Two strips 10 m wide immediately adjacent to the rainshelter were planted in crop to minimize edge effects. The soil was a deep (greater than 1.6 m) Templeton sandy loam (*Udic Ustochrept*; U.S. Department of Agriculture soil taxonomy) with an available water-holding capacity of 190 mm m^{−1} depth (Martin et al., 1992). Nitrogen fertilizer was applied through the irrigation system as urea dissolved in the irrigation water. N fertilizer rate was determined with the objective to obtain a grain yield of 10 t DM (dry mass) ha^{−1}. The experiment consisted of a randomized complete block design with two replicates. Six irrigation treatments were designed to subject the crops to drought of varying durations at different times during the growth of the crops. A control treatment was fully irrigated, and another one received no water so that the crop grew completely on stored water. Details of the irrigation treatments are given by Jamieson et al. (1995b). Plant material in a 0.2-m^{2} area per subplot was sampled by cutting at ground level at 2-week intervals until anthesis and thereafter at 5-d intervals. The leaf area and dry mass of a 10-tiller subsample were measured, and LAI was calculated from the leaf area ratio and the total sample dry mass.

The last experiment was conducted at Clermont-Ferrand, France, and at Sutton Bonington, United Kingdom, during two consecutive winter cropping cycles (experiment E in Table I; Gaju et al., 2011; Moreau et al., 2012). At both site, crops were sown at the recommended dates. Two rates of N fertilizer were used. The high-N treatment was intended to replicate commercial practice, and the rates of N fertilization were determined using the balance-sheet method to optimize grain yield (Rémy and Hébert, 1977). N was applied as ammonium nitrate granules in three (at Sutton Bonington) to four (at Clermont-Ferrand) splits. The amount of N applied in the low-N treatment was adjusted in each site-season according to the soil mineral N measured at the end of the winter period, with the aim of providing 100 kg N ha^{−1} from the combined soil mineral N and fertilizer N, corresponding to a moderate to severe N limitation sufficient to reduce grain yield by approximately 30% compared with the high-N treatment. All the crops were rainfed. At Sutton Bonington, plant growth regulator was applied as chlormequat at the onset of stem extension. The data used here are for cv Soissons, but the experiment comprised 16 genotypes, and at both sites, a split-plot design was used in which N treatment was randomized on main plots, genotypes were randomized on the subplots, and each treatment was replicated three times. Subplot size was 24 × 1.65 m at Sutton Bonington and 7 × 1.5 m at Clermont-Ferrand. Interrow spacing was 0.175 and 0.125 m at Clermont-Ferrand and Sutton Bonington, respectively. Total LAI and GAI were determined between one and five times around anthesis (Moreau et al., 2012) following the procedure described above for experiment C. LAI of the culm leaf cohorts (i.e. the last five leaves) also was determined at Clermont-Ferrand in both growing seasons and at Sutton Bonington in the second year.

### Model Description

The equations describing the final potential sizes of laminae and sheaths are given in “Results” (Eqs. 1 and 2). The equations and assumptions concerning the rate and duration of leaf expansive growth and senescence and their responses to water and N deficit are presented below. Values of all parameters of the leaf area model are given in Supplemental Table S1.

The leaf area model was developed as an independent executable component in the BioMA software framework (http://www.biomamodelling.org) and was integrated in the wheat model SiriusQuality. The source code and the binaries of SiriusQuality can be freely downloaded at http://www1.clermont.inra.fr/siriusquality/. The source code of the leaf area component can be obtained upon request from the corresponding author.

#### Duration and Rate of Leaf Expansion and Senescence

The length of wheat leaf laminae increases almost linearly after emergence of the tip of the leaf above the whorl of subtending leaves until the appearance of the leaf collar, after which it decreases rapidly (Gallagher, 1979; Fournier et al., 2005). When normalized by their mature length and plotted against phyllochronic time centered on the time of leaf tip emergence, a unique kinetic of leaf expansion is observed for all the leaf and tiller ranks (Fournier et al., 2005). According to Gallagher (1979), the phyllochronic duration of lamina expansion is close to 1.1 and that of the sheath is close to 0.3. Based on these results, a unique phyllochronic duration of lamina expansion is considered to model the expansive growth of the laminae of the different leaf ranks and tiller orders. Thus, in our model, the potential daily increase in lamina surface area () is calculated as:where (°Cdays) is the daily integral of thermal time, *P* (°Cdays) is the phyllochron (i.e. the thermal time requirement for the appearance of successive leaf tips), and *t*^{exp} (dimensionless) is the number of phyllochrons between the appearance of the tip and the collar of leaf *n* above the collar of leaf *n* − 1.

Assuming that the relative rate of expansion with thermal time is the same for the laminae and the sheaths, the duration of sheath expansion above the collar of the preceding leaf is calculated as:Daily minimum () and maximum () canopy temperatures are calculated using an energy balance assuming a neutral atmospheric stability approach (Jamieson et al., 1995a), and is calculated as the sum of eight contributions each day of a cosinusoidal variation between and modified from Weir et al. (1984):whereandwhere *T*_{min} and *T*_{opt} (°C) are the minimum (base) and optimum temperatures for leaf development and expansive growth, *f*(*T*) (dimensionless) is the temperature response function for leaf development and expansive growth, (°C) is the calculated 3-hourly canopy temperature contribution to estimated daily mean canopy temperature, *f _{r}* (dimensionless) is the fraction that each 3-h period during the day contributes to the thermal time for that day, and

*r*is the array index of the item.

Recent studies showed that all developmental and expansive growth processes follow a common curvilinear response to temperature after normalization by a common reference temperature (Parent et al., 2010; Parent and Tardieu, 2012; Wang et al., 2017). To model the temperature response of leaf growth, we use the nonlinear temperature function proposed by Wang and Engel (1998):wherewhere *T*_{max} (°C) is the maximum canopy temperature for leaf development and expansive growth. Equations 9 and 10 simulate the effect [0-1] of temperature between *T*_{min} and *T*_{max} and is used with the same cardinal temperature values to model the duration of expansion and the rate of leaf appearance.

As for the duration of the expansive growth period, the potential durations of the period during which the green area of a leaf is constant (mature period) and of the period during which its green area regresses due to senescence (senescence period) are assumed to be constant in phyllochronic time. The duration (in phyllochronic time) of the mature and senescence periods is longer for mature leaves than for juvenile leaves (see below). We assumed that the lamina and the sheath senesce sequentially, as is usually observed in the field.

#### Tillering

The delay between the emergence of the first leaf of a tiller and the emergence of the first leaf of its parent tiller is approximately constant to three phyllochrons (Kirby et al., 1985). Once the first tiller has emerged, the next tillers appear with a high probability (Whaley et al., 2000) and follow a fixed pattern according to the Fibonacci series (Kirby et al., 1985; Evers et al., 2007). In our model, tiller outgrowth follows the Fibonacci series, and based on the results presented here, leaves of the different tillers competed acropetally (i.e. leaf cohorts) and have the same characteristics and constitute homogenous layers (after shifting the phytomers by a constant number).

In the field, as the canopy develops, both a decrease of the red/far-red ratio (Evers et al., 2006; Sparkes et al., 2006) and the beginning of internode extension (Miralles and Richards, 2000) are able to trigger the end of tillering. Soon after the end of the period of tiller production, a proportion of tillers regresses until ear emergence (Alzueta et al., 2012). This proportion is highly variable and depends on the availability of carbohydrates (Dreccer et al., 2013), water (Davidson and Chevalier, 1990), and nutrients (Rodríguez et al., 1998; Alzueta et al., 2012). It is usually observed that the youngest tillers start to regress first (Davidson and Chevalier, 1990). Within a plot, plant-to-plant competition and local heterogeneity in plant density and resources often result in large variations of final tiller number between individual plants (Masle-Meynard and Sebillotte, 1981). However, at the canopy level, the number of fertile tillers is largely independent of the sowing density for sowing density greater than 100 seeds m^{−2} (Spink et al., 2000; Whaley et al., 2000), because of a compensation between tiller production and tiller survival (Alzueta et al., 2012). Therefore, here, we did not attempt to explicitly model the cessation of tillering and tiller death. Instead, we assumed that canopies have a fixed number of fertile tillers of 550 tillers m^{−2}. The number of tillers produced is calculated according to the plant density (PD; plants m^{−2}) to reach a fixed number of fertile tillers. For instance, if PD is equal to 200 plants m^{−2} and the final number of fertile shoots is 550 tillers m^{−2}, then the modeled crop will produce 200 mainstems and first-order tillers per square meter and 150 second-order tillers per square meter.

#### Water- and N-Limited Leaf Expansive Growth

The potential daily increase in green area index of the canopy (; m^{2} m^{−2}) is calculated by summing the potential daily increase in lamina () and sheath () surface area of all growing leaves multiplied by PD:The actual daily increase in green area index () is calculated from according, first, to soil water availability and then to N availability. For most crop species, leaf expansion decreases linearly as the fraction of available soil water (FPAW; dimensionless) decreases below a threshold value (Meyer and Green, 1981; Robertson and Giunta, 1994; Sadras and Milroy, 1996). Therefore, we modeled the effect of soil water deficit on leaf expansion using a bilinear stress function:where (m^{2} m^{−2}) is the water-limited daily change in canopy green area index, (dimensionless) is the threshold FPAW value at which the rate of leaf expansion starts to decrease, and (dimensionless) is the FPAW value at which the rate of leaf expansion equals zero.

Expanding grass leaves require a minimum nitrogen content and do not store nitrogen during their growth period (Gastal and Lemaire, 2002), and it was shown that the area-based leaf nitrogen mass of the light-exposed part of expanding leaves is constant and independent of the N status of the crop (van Oosterom et al., 2010). Based on these results, in our model, leaf expansion is reduced if there is not enough nitrogen available in the plant to maintain a critical area-based leaf nitrogen mass (; g N cm^{−2} leaf). Thus, the water- and nitrogen-limited daily increase in green area index (; m^{2} leaf m^{−2} ground) is calculated as:where (g N m^{−2} ground) is the mass of labile (remobilizable) N in the crop. The mass of labile N is calculated as the total mass of N in the crop minus the mass of structural N and plus N taken from the soil. The concentration of structural N is constant and is different for internodes (3 × 10^{−3} g N g^{−1} DM) and leaves (5 × 10^{−3} g N g^{−1} DM; Martre et al., 2006).

Finally, the water- and nitrogen-limited daily rates of expansion of each growing leaf lamina () and sheath () are calculated according to their potential rate of expansion and to the water- and nitrogen-limited rate of expansion of the whole canopy:and

#### Leaf Senescence

The durations of the mature and senescence periods are calculated using a constant phyllochronic duration modified by a drought stress factor (Semenov et al., 2009). In wheat, heat stress significantly accelerates the rates of leaf aging and senescence (Wollenweber et al., 2003; Tewolde et al., 2006; Zhao et al., 2007; Eyshi Rezaei et al., 2015), and in our model, this effect was considered by accelerating the thermal time as proposed by Asseng et al. (2011). Leaf aging during the mature period is modeled as a time integral (*D*_{lag}; dimensionless). A leaf’s senescence starts when *D*_{lag} has reached a value of 1.wherewhere (°Cdays) is the daily integral of thermal time for leaf aging and senescence, and (dimensionless) are the potential phyllochronic durations of the mature period for the juvenile and mature leaves, respectively, DSF (dimensionless) is a drought stress factor that accelerates leaf aging and senescence, is the threshold FPAW value at which the rate of leaf senescence is accelerated, is the FPAW value at which the maximum value of DSF is reached, and DSF_{max} is the maximum value of DSF.

Similar to the mature period, during the senescence period, the daily change in leaf area is accelerated by water deficit and heat stress:where and (dimensionless) are the potential phyllochronic durations of the leaf senescence period for the juvenile and mature leaves, respectively. As in the previous version of SiriusQuality, during the grain-filling period, N dynamics accelerate leaf senescence in low-N crops (Martre et al., 2006).

To account for the shortening of the mature and senescence phases caused by heat, the 3-hourly canopy temperatures used to calculate the duration of these two leaf ontogenic phases are multiplied by an accelerated leaf senescence factor (; dimensionless), which increases linearly from 1 when exceeds a threshold temperature (*T*^{L}; °C; Stratonovitch and Semenov, 2015):wherewhere *S*^{L} (°C^{−1}) is the slope of the senescence acceleration per unit of canopy temperature above *T*^{L}. As in the original version of SiriusQuality, grain filling is stopped prematurely if the canopy has fully senesced.

#### Model Calibration and Evaluation

The four parameters describing the potential size of the mature leaf laminae and sheaths were defined as genotypic parameters. All the parameters of the leaf area model and the three genotypic parameters of the phenology model (He et al., 2012) used to evaluate the model are given in Supplemental Table S1. For experiments C and D, measured anthesis date, final leaf number, and LAI measured in the high-N or well-watered treatments were used to estimate the genotypic parameters. For experiment E, we used the genotypic parameters of the leaf area model estimated independently with experiments A and B for cv Soissons and the three genotypic parameters of the phenology model estimated by He et al. (2012) for this cultivar. The other parameters of the leaf area model were obtained from the literature and were set constant for all cultivars (Supplemental Table S1).

To assess the quality of the model, measured () and simulated () values were compared using ordinary least square regression and the RMSE. The RMSE was calculated as the square root of the MSE:The MSE was further decomposed into its three components: the nonunity slope (NU), the squared bias (SB), and the lack of correlation (LC; Gauch et al., 2003):where *b* is the slope of the regression of on and *r*^{2} is the coefficient of correlation.

The three components of the MSE, which add up to give MSE, represent different aspects of the overall deviation of the model simulations and have simple geometrical interpretations. NU reflects the rotation, SB the translation, and LC the scattering (random error) around the 1:1 line. This analysis was used as a complement of the classical least square linear regression.

### Supplemental Data

The following supplemental materials are available.

**Supplemental Figure S1.**Relationships between final sheaths and laminae area.**Supplemental Table S1.**Values of the parameters used for the model evaluation.**Supplemental Table S2.**Statistics of the model evaluation.

## Acknowledgments

We thank Philippe Gate (ARVALIS Institut-du-Végétal) for sharing unpublished data and Bruno Andrieu (Unité Mixte de Recherche EcoSys, Institut National de la Recherche Agronomique) and Christian Fournier (Unité Mixte de Recherche Laboratoire d'Ecophysiologie des Plantes sous Stress Environnementaux, Institut National de la Recherche Agronomique, Montpellier SupAgro) for sharing data and for useful discussion and comments on the model.

## Footnotes

The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (www.plantphysiol.org) is: Pierre Martre (pierre.martre{at}inra.fr).

P.M. conceived the research, analyzed the experimental data, and developed the model; A.D. performed the simulations; P.M. and A.D. wrote the article.

↵1 This work was supported by the European Union’s Seventh Framework Program (FP7/2007-2013; grant no. FP7-613556).

↵2 These authors contributed equally to the article.

- Received July 18, 2017.
- Accepted November 12, 2017.
- Published November 15, 2017.