# L-GrassF: a functional-structural and phenological model of *Lolium perenne* integrating plant morphogenesis and reproductive development

Simon Rouet<sup>1,2,3,*</sup>, Jean-Louis Durand<sup>1</sup>, Denis Leclercq<sup>4</sup>, Marie-Hélène Bernicot<sup>5</sup>, Didier Combes<sup>1,●</sup>, Abraham Escobar-Gutiérrez<sup>1,●</sup> and Romain Barillot<sup>1,●</sup>

Handling Editor: Jochem Evers

Editor-in-Chief: Stephen P. Long

**Citation:** Rouet S, Durand J-L, Leclercq D, Bernicot M-H, Combes D, Escobar-Gutiérrez A, Barillot R. 2022. L-GrassF: a functional-structural and phenological model of *Lolium perenne* integrating plant morphogenesis and reproductive development. *In Silico Plants* 2022; diac012; doi: 10.1093/insilicoplants/diac012

## ABSTRACT

In the context of climate change and agrosystem complexification, process-based models of the reproductive phenology of perennial grasses are essential to optimize the agronomic and ecologic services provided by grasslands. We present a functional-structural model called L-GrassF, which integrates the vegetative and reproductive development of individual *Lolium perenne* plants. The vegetative development in L-GrassF was adapted from a previous model of perennial ryegrass where leaf elongation and tillering dynamics partially result from self-regulated processes. Significant improvements have been made to this vegetative module in order to deal with the whole growing cycle during which plants are exposed to contrasting temperatures. The reproductive module is a new functionality describing the floral induction of the individual tiller from daily temperature and photoperiod as well as its phenological state. From the interactions between the vegetative and reproductive developments, L-GrassF simulates the dynamics of plant architecture, the floral transition and heading date (HD) at tiller level. A sensitivity analysis was performed on L-GrassF and showed that parameters controlling the kinetics of leaf elongation and leaf appearance rate have a significant impact on HD. After calibration, L-GrassF was able to simulate the HD on seven *L. perenne* cultivars grown in a broad range of environmental conditions, as provided by an independent data set. We conclude that L-GrassF is a significant step towards better prediction of grassland phenology in contrasted conditions.

**KEYWORDS:** Floral transition; grassland; individual-based model; perennial grass; photoperiod; temperature.

## 1. INTRODUCTION

The reproductive phenology focuses on the series of periodical events occurring during the reproductive development of plants. In perennial grasses, it deeply affects plant morphogenesis, biomass production in terms of quantity and quality, as well as grassland perenniality

(Rouet et al. 2021). In spring, an increase of plant growth is commonly observed, which cannot only be explained by temperature conditions and was related to the phenological state of the plants (Davies 1971; Parsons and Robson 1980). The acceleration of forage biomass production goes along with a decrease of biomass quality mainly due to

<footnote n="1"><sup>1</sup>INRAE, URP3F, 86600 Lusignan, France</footnote>

<footnote n="2"><sup>2</sup>CIRAD, UPR GECO, F-97130 Capesterre-Belle-Eau, Guadeloupe, France</footnote>

<footnote n="3"><sup>3</sup>GECO, Univ. Montpellier, CIRAD, 34000 Montpellier, France</footnote>

<footnote n="4"><sup>4</sup>GEVES, 86600 Lusignan, France</footnote>

<footnote n="5"><sup>5</sup>GEVES, 21100 Bertenières, France</footnote>

<footnote>*Corresponding author's e-mail address: simon.rouet@cirad.fr</footnote>
an increased proportion of internodes in the biomass, which lowers digestibility for ruminants (Chapman et al. 2014). The tiller demography, defined as the number of tillers per unit area and their characteristics (age, phenological state...), is also affected by reproductive phenology as the reproductive tillers die at the end of the reproductive development, while vegetative ones survive and ensure plant perenniality (Matthew et al. 2000). The heading of reproductive tillers, i.e. the appearance of the spike tip outside the pseudo-stem formed by the sheaths (Gillet 1980; Thomas 1980), usually occurs in late spring in temperate regions. It is an easily observable phenological state used as a routine indicator for grassland management and ranking cultivar earliness.

Heading of reproductive tillers results from processes which start several months before. First, the transition of the tiller from a vegetative to reproductive development is subordinated to the floral induction of its terminal apex. In *Lolium perenne*, a major temperate perennial grass, the floral induction follows dual induction with two strictly successive phases: (i) the primary induction requiring several weeks of exposure to low temperatures and short photoperiods, (ii) the secondary induction requiring a few days of exposure to long photoperiods and high temperatures (Heide 1994). In all cases, the floral transition occurs at the end of the secondary induction and determines the start of reproductive development. After floral transition, the length of the apex increases as well as the rate of primordium production (Cooper 1950; Gonthier and Francis 1989; Kemp et al. 1989). The primordia accumulated at the apex after floral transition forms the last leaves of the tiller for the most basal ones, while the most distal ones will form spikelets (Gonthier and Francis 1989). Heading takes place only after all leaves have expanded and is the result of the combined growth of the spike and the peduncle below it.

It appears that the reproductive phenology of perennial grasses is the result of a complex interaction between the effect of environmental conditions on floral induction and growth processes, both exhibiting high genetic variability (Aamlid et al. 2000; Keep et al. 2020). Computer modelling is required to understand and predict perennial grass phenology, which will further help to optimize grassland management and the related ecosystem services (Durand et al. 2016; Kipling et al. 2016). Existing grassland models are highly empirical about phenology, which is a constraint for generic usage, as they should potentially be recalibrated for each genetic × environment combination (for a review, see Rouet et al. 2021). In these models, phenology is either dependent on thermal time (Bonesmo and Bélanger 2002; Jouven et al. 2006; Jégo et al. 2013) or is regulated according to a fixed calendar provided in inputs (Fiorelli et al. 2001; Barrett et al. 2005). In the current context of climate change and agrosystem complexification, phenological models of perennial grasses have to better consider the environmental drivers of reproductive development, the genetic variability of the response to these signals at the biologically most relevant scale.

The aim of this study is to propose an original and generic model of perennial grass phenology, which accounts for the interactions between the environmental factors, floral induction and tiller morphogenesis. We developed a functional-structural plant model (FSPM) named L-GrassF which represents the complete development of individual tillers from their initiation to heading. One of the main features

of the model is to consider the interactions between the reproductive and vegetative development. The model simulates the main reproductive stages (floral transition, heading date [HD]) and morphological variables (e.g. final leaf number and dimensions, number of spikelets per spike). Here, we analysed the behaviour of L-GrassF and its ability to simulate the HD of different cultivars of *L. perenne* grown under contrasted natural conditions. To this purpose, a sensitivity analysis was performed in order to determine the parameters with the highest effect on the HD. The most sensitive parameters were then estimated for seven cultivars on a subset of a data set provided by GEVES, the French Variety and Seed Study and Control Group. Finally, L-GrassF was evaluated for its ability to predict the HD on the second subset of data. The present analyses were performed on the main stems of plants only, as the available data did not allow us to assess the ability of the model to simulate the proportion of reproductive tillers.

# 2. MATERIALS AND METHODS

## 2.1 Model description

2.1.1 Model overview. L-GrassF is an FSPM of perennial grasses integrating the reproductive phenology as the result of the interaction between the vegetative development, the reproductive development and environmental factors. The visible expression of reproductive phenology, used to manage grasslands or describe the species variability is heading, which is the main phenological prediction for which we calibrated the model in the present study. L-GrassF therefore simulates, in relation to environmental conditions, the vegetative (appearance and elongation of leaves, tillering) and reproductive development (floral inductions, floral transition, internode and spike production, heading) of individual tillers from their initiation to heading. L-GrassF was developed from L-Grass, an existing FSPM which describes the aerial vegetative development in *L. perenne* (Verdenal et al. 2008; Verdenal 2009).

Briefly, the initial version of L-Grass describes the dynamics of leaf growth and tillering, as well as the effects of defoliation for each individual plant of the canopy. In L-Grass, the plant is seen as a set of individual tillers, each being represented as a succession of phytomers, i.e. an internode, an axillary bud and a leaf (made of a lamina and a sheath). The successive nested sheaths constitute a so-called pseudo-stem through which the leaves grow before they emerge. The above-ground morphogenesis was implemented as a self-regulated system: (i) the rhythm of leaf initiation is coordinated with leaf emergence, (ii) leaf elongation rate and final dimensions depend on the time spent growing in the pseudo-stem and (iii) the rhythm of tiller emergence is coordinated with that of the leaves on the mother tiller, while the probability of tiller emergence is negatively regulated by the leaf area of the plant and its direct neighbours. The model was calibrated for *L. perenne* and is able to produce contrasted phenotypes from turf to forage types. The model is based on the L-system formalism, which allows a 3D and botanic representation of plant dynamics (Prusinkiewicz 2004). It was implemented in L-Py, an L-system framework in Python distributed in the OpenAlea platform (Boudon et al. 2012).

L-GrassF was built upon the version of L-Grass implemented in the L-Py framework. The decomposition of the plant into individual tillers and organs was kept identical to that of L-Grass (Fig. 1). The
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Overview of L-GrassF functioning. L-GrassF accounts for tillering as well as leaf initiation, elongation and emergence for each tiller in relation with air temperature (1, 2). A phenological model was developed, it calculates the rate of primary and secondary induction according to temperature and photoperiod (3, 4). The progress of floral induction defines the rate of primordium initiation at the apex (5). Once the floral induction is completed, the floral transition is carried out and the final number of leaves is defined. The final number of spikelets is defined when the last leaf starts to elongate (6). Finally, after the elongation of all leaves is complete (7), the model simulates the heading of the spike from the ligulation date of the flag leaf of the tiller (8).</figcaption></figure>

main principles used to model the vegetative development in L-GrassF were derived from L-Grass, which however required several adaptions in order to extend L-GrassF to the reproductive stages. Briefly, we had to adapt two main aspects: (i) the time unit of the model was converted from thermal time to astronomic time and temperature-dependent processes were adapted by integrating a non-linear response function to temperature (Section 2.1.2; see **Supporting Information**—Fig. S1) and (ii) the kinetics of leaf growth prior to its emergence and the determination of the final leaf length were adapted to account for plants that remain vegetative for long periods, while the initial model was limited to 1500 °C day (Verdenal et al. 2008), which is not long enough to represent the complete reproductive development in natural temperature conditions (Section 2.1.3.2). In order to account for the transition from vegetative to reproductive stages, a new submodel of floral induction was developed at tiller scale and implemented in L-GrassF. Floral

induction is represented as two successive phases: primary induction depending on temperature (Section 2.1.4.1; Fig. 1.3) and secondary induction depending on photoperiod (Section 2.1.4.2; Fig. 1.4). Floral induction dynamically affects tiller morphogenesis (Section 2.1.5; Fig. 1.5, 6, 7) thereby determining the final morphology of the tiller (final number of leaves, number of long internodes and spikelets). Because the dynamics of internode and spike elongation are poorly documented for perennial grasses, we chose to approximate the HD as the date of ligulation of the flag leaf rather than from geometrical calculations (Fig. 1.8). L-GrassF was calibrated on *L. perenne* but the underlying processes implemented are generic to grasses. The code of the model is open source, available on GitHub (https://github.com/openalea-incubator/lgrass), and released under a CecILL-C license. The archive of the code used to generate the present results is available at https://doi.org/10.5281/zenodo.6873725.
### 2.1.2 Time unit of the model and temperature effects.
In L-GrassF, the time unit is astronomical time (expressed in hours), contrasting with the initial version of the model which was based on thermal time. Since the model is intended to be used in contrasting environments with variable temperature ranges and some processes like floral induction are directly dependent on temperature, it appears that thermal time based on the sum of daily temperatures is no longer a relevant time unit (Bonhomme 2000). Temperature-dependent processes were adapted by integrating a non-linear function of temperature (Yan and Hunt 1999; Zaka et al. 2017, pers. comm.).

Briefly, the kinetics of leaf elongation and senescence now depend on leaf physiological age (Age<sub>n</sub>; h), which accounts for a non-linear effect of temperature on leaf ageing. Leaf physiological age increases as a function of temperature as described in Equation (1) and **Supporting Information—Fig. S1**.

$$ \frac{dAge_n}{dt} = \begin{cases} 0 & \text{if } T_{\text{min}} \le T(t) \le T_{\text{max}} \\ \left(\frac{T(t)-T_{\text{min}}}{T_{\text{ref}}-T_{\text{min}}}\right)^q \times \left(\frac{T_{\text{max}}-T(t)}{T_{\text{max}}-T_{\text{ref}}}\right); & \text{else} \end{cases} \quad (1) $$

where T(t) is the mean temperature of the day (°C), T<sub>min</sub> and T<sub>max</sub> are the minimum and maximum temperatures below and above which leaf age does not increase, q is a shape parameter and T<sub>ref</sub> is the reference temperature (Table 1). According to the observations of Zaka (2016), we assumed that cardinal temperatures were identical for all genotypes, which is consistent with previous findings in a wide range of cultivated species (Parent and Tardieu 2012).

### 2.1.3 Vegetative morphogenesis.
This section describes tiller functioning during the vegetative stages. Although this part of the model is mainly based on the initial L-Grass model, we had to adapt the functions of leaf elongation in order to further extend the model to reproductive stages.

#### 2.1.3.1 Dynamics of apex functioning.
In L-GrassF, each tiller has a module called 'apex' declared in the L-system string which represents the terminal apical meristem. During the vegetative stages, the tiller apex initiates a new phytomer of rank n (leaf, internode and axillary bud) at the emergence of leaf n-1 (Fig. 2). Thus, the phyllochron and plastochron are equal but they depend upon the time taken by the leaf to exceed the length of the pseudo-stem, which usually corresponds to the sheath length of the previous leaf.

#### 2.1.3.2 Leaf elongation.
The elongation of each leaf is divided into two distinct phases: from leaf initiation to leaf emergence ('hidden phase') and from leaf emergence to the end of leaf elongation ('visible phase').

**Hidden phase**

The length of a leaf n ($Y_n$, mm) is initiated to 1 mm and increases following a 'dual' function depending on leaf age (Equation (2A)). While leaf age is lower than a threshold ($Age_{max}$, h), leaf length follows a beta function (Equation (2B)) which tends to $Y_m^0$ (mm), as proposed in the initial version of L-Grass (Verdenal et al. 2008). If leaf age exceeds $Age_{max}$, leaf length then increases linearly with leaf age (Equation (2A)). The slope of the linear function is equal to the beta function derived at $Age_{max}$.

$$ Y_n(t) = \begin{cases} \frac{\text{Beta}(t)}{\frac{d\text{Beta}(Age_{\text{max}})}{dt}} & \text{for } Age_n \le Age_{\text{max}} \\ (Age_n(t) - Age_{\text{max}}) + \text{Beta}(Age_{\text{max}}); & \text{for } Age_n > Age_{\text{max}} \end{cases} \quad (2A) $$

$$ \text{Beta}(t) = Y_m^0 \left( 1 + \frac{\text{LED}^0 - \text{Age}_n(t)}{\text{LED}^0 - \text{Age}_{\text{max}}} \right) \left( \frac{\text{Age}_n(t)}{\text{LED}^0} \right)^{\frac{\text{LED}^0}{\text{LED}^0 - \text{Age}_{\text{max}}} \quad (2B) $$

where $Y_m^0$ (mm) is the asymptote of the beta function which is reached at LED<sup>0</sup> (h) (Table 1).

### Leaf emergence and final length determination

The leaf emerged when it exceeds the length of the pseudo-stem. The length of the pseudo-stem is calculated as the length of the longest sheath of the tiller. At leaf tip emergence, the final length of the leaf (Ym) is defined as well as the duration of its elongation (LED<sub>0</sub>) (Verdenal et al. 2008).

First, final leaf length (Ym<sub>n</sub>, mm) is assumed to increase with the time spent growing in the pseudo-stem until a maximum is reached, if this duration exceeds $Age_{max}$ (Equation (3)).

$$ Y_{m_n} = C \times \left( T_{01} e^{k_1 (min(Age_n(t_{E_n}), Age_{max}) + T_{b_1})} - T_{02} e^{k_2 (min(Age_n(t_{E_n}), Age_{max}) + T_{b_2})} - L \right) \quad (3) $$

where C is a scaling factor, $Age_n(t_{E_n})$ (h) is leaf n age at emergence and $T_{01}$ (mm), $T_{02}$ (mm), $k_1$ (h<sup>-1</sup>), $k_2$ (h<sup>-1</sup>), $T_{b_1}$ (h), $T_{b_2}$ (h) and L (mm) are parameters fitted to data (Verdenal et al. 2008). Parameters were taken from Verdenal et al. (2008) except for time related ones ($k_1, k_2, T_{b_1}$ and $T_{b_2}$) which were converted from degree days to hours (Table 1). In contrast to the original version of L-Grass, we also assumed that final leaf length reaches a maximum if the duration of leaf elongation during the initial phase exceeds a threshold. With the exception of $Y_m^0$, all parameters were set to constant values between cultivars in the following simulations (Table 1).

Second, leaf elongation duration (LED<sub>0</sub>) is assumed to increase with the sheath: lamina length ratio, which is itself calculated as a linear function of leaf age (for more detail, see in Verdenal et al. 2008).

**Visible phases**

After leaf emergence, leaf length is calculated from the beta function defined in (Equation (2B)) using the parameters $Y_m$ and LED in place of $Y_m^0$ and LED<sup>0</sup>, respectively. Leaf elongation stops when leaf length reaches $Y_m$.

#### 2.1.3.3 Cutting.
One of the main features of grassland management is their frequent cutting. The cutting module implemented in L-GrassF is similar to that described in (Verdenal et al. 2008): every organ segment located above a horizontal plane, mimicking the cut height, is removed. If the cutting shortens the pseudo-stem of a given tiller, it will affect the emergence dates and final lengths of the non-emerged growing leaves of the tiller. If an apex is located higher than the cutting height, it will trigger the death of the tiller.

#### 2.1.4 Floral induction.
The floral induction proceeds in two independent phases, each explained by a single variable. This accounts for the observations on the floral induction in a sufficiently realistic way and
<table><caption>Table 1. Model parameters. *Values were studied by the sensitivity analysis.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Unit</th><th>Related equation</th><th>Reference</th></tr></thead><tbody><tr><td>Organ ageing</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>T<sub>min</sub></td><td>Minimum temperature</td><td>-2.88</td><td>°C</td><td>Equation (1)</td><td>Zaka et al. (2017), pers. comm.</td></tr><tr><td>T<sub>max</sub></td><td>Maximum temperature</td><td>38.9</td><td>°C</td><td>Equation (1)</td><td>Zaka et al. (2017), pers. comm.</td></tr><tr><td>T<sub>ref</sub></td><td>Reference temperature</td><td>23.97</td><td>°C</td><td>Equation (1)</td><td>Zaka et al. (2017), pers. comm.</td></tr><tr><td>q</td><td>Shape parameter</td><td>1.56</td><td>Dimensionless</td><td>Equation (1)</td><td>Zaka et al. (2017), pers. comm.</td></tr><tr><td>Leaf elongation</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Ω</td><td>Parameter to calculate LED<sup>a</sup> and LED</td><td>1.91</td><td>h mm<sup>-1</sup></td><td>—</td><td>Adapted from Verdenal et al. (2008)</td></tr><tr><td>a</td><td>Parameter to calculate LED</td><td>1.18 × 10<sup>-3</sup></td><td>h<sup>-1</sup></td><td>—</td><td>Adapted from Verdenal et al. (2008)</td></tr><tr><td>b</td><td>Parameter to calculate LED</td><td>-0.1</td><td>Dimensionless</td><td>—</td><td>Verdenal et al. (2008)</td></tr><tr><td>R<sub>ta</sub></td><td>Parameter to calculate LED</td><td>0.2</td><td>Dimensionless</td><td>—</td><td>Verdenal et al. (2008)</td></tr><tr><td>Hidden phase</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LED<sup>0</sup></td><td>Leaf age at which the beta function is maximal</td><td>Y<sub>m</sub><sup>0</sup> × Ω</td><td>h</td><td>Equation (2)</td><td>Verdenal et al. (2008)</td></tr><tr><td>Age<sub>max</sub></td><td>Leaf age threshold</td><td>0.76 × LED<sup>0</sup></td><td>h<sup>-1</sup></td><td>Equation (2)</td><td>Verdenal et al. (2008)</td></tr><tr><td>Visible phase</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>T<sub>01</sub></td><td>Parameter to calculate Y<sub>m</sub><sup>0</sup></td><td>2.718</td><td>mm</td><td>Equation (3)</td><td>Verdenal et al. (2008)</td></tr><tr><td>T<sub>02</sub></td><td>Parameter to calculate Y<sub>m</sub><sup>0</sup></td><td>2.718</td><td>mm</td><td>Equation (3)</td><td>Verdenal et al. (2008)</td></tr><tr><td>k<sub>1</sub></td><td>Parameter to calculate Y<sub>m</sub><sup>0</sup></td><td>2.249 × 10<sup>-2</sup></td><td>h<sup>-1</sup></td><td>Equation (3)</td><td>Adapted from Verdenal et al. (2008)</td></tr><tr><td>k<sub>2</sub></td><td>Parameter to calculate Y<sub>m</sub><sup>0</sup></td><td>3.375 × 10<sup>-2</sup></td><td>h<sup>-1</sup></td><td>Equation (3)</td><td>Adapted from Verdenal et al. (2008)</td></tr><tr><td>Tb<sub>1</sub></td><td>Parameter to calculate Y<sub>m</sub><sup>0</sup></td><td>144</td><td>h</td><td>Equation (3)</td><td>Adapted from Verdenal et al. (2008)</td></tr><tr><td>Tb<sub>2</sub></td><td>Parameter to calculate Y<sub>m</sub><sup>0</sup></td><td>27.98</td><td>h</td><td>Equation (3)</td><td>Adapted from Verdenal et al. (2008)</td></tr><tr><td>L</td><td>Parameter to calculate Y<sub>m</sub><sup>0</sup></td><td>30</td><td>mm</td><td>Equation (3)</td><td>Verdenal et al. (2008)</td></tr><tr><td>Floral induction</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Primary induction</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>T<sub>min</sub><sup>I<sup>R</sup></sup></td><td>Minimal temperature for primary induction</td><td>-3</td><td>°C</td><td>Equation (5)</td><td>Woodward et al. (2020)</td></tr><tr><td>T<sub>opt</sub><sup>I<sup>R</sup></sup></td><td>Optimal temperature for primary induction</td><td>4</td><td>°C</td><td>Equation (5)</td><td>Woodward et al. (2020)</td></tr><tr><td>T<sub>max</sub><sup>I<sup>R</sup></sup></td><td>Maximal temperature for primary induction</td><td>12</td><td>°C</td><td>Equation (5)</td><td>Woodward et al. (2020)</td></tr><tr><td>Δind<sub>max</sub><sup>I<sup>R</sup></sup></td><td>Maximal daily increment for primary induction</td><td>*</td><td>Dimensionless</td><td>Equation (5)</td><td>Calibrated on data</td></tr><tr><td>Secondary induction</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PP<sub>max</sub></td><td>Minimal photoperiod for secondary induction</td><td>*</td><td>h</td><td>Equation (7)</td><td>Calibrated on data</td></tr><tr><td>L/NM</td><td>Minimum number of leaves for secondary induction at maximum daily rate</td><td>*</td><td>Dimensionless</td><td>Equation (7)</td><td>Calibrated on data</td></tr><tr><td>PPRM</td><td>Maximum daily rate of the secondary induction</td><td>*</td><td>Dimensionless</td><td>Equation (7)</td><td>Calibrated on data</td></tr><tr><td>Reproductive morphogenesis</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>index<sub>n</sub><sup>R</sup></td><td>Threshold of the primary phase for leaf elongation increase</td><td>1</td><td>Dimensionless</td><td>—</td><td>Calibrated on data</td></tr><tr><td>index<sub>n</sub><sup>J<sup>R</sup></sup></td><td>Threshold of the secondary phase for leaf elongation increase</td><td>1</td><td>Dimensionless</td><td>—</td><td>Calibrated on data</td></tr><tr><td>kY<sub>m</sub><sup>0</sup></td><td>Multiplicative factor of Y<sub>m</sub><sup>0</sup></td><td>*</td><td>Dimensionless</td><td>—</td><td>Calibrated on data</td></tr><tr><td>ρ</td><td>Number of primordia initiated at leaf emergence</td><td>*</td><td>primordium per leaf</td><td>—</td><td>Calibrated on data</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 2. Dynamics of apex functioning during the floral induction process in L-GrassF. In this schematic representation of the apical zone, light green circles are leaf primordia, dark green structures are leaves growing in the pseudo-stem, light grey structures are mature leaves emerged from the pseudo-stem, yellow circles are spikelet primordia and yellow ovals are spikelets. Dotted red lines indicate newly formed structures. Before the secondary induction, phyllochron and plastochron are equal, i.e. one leaf primordium is produced at each leaf emergence. From the start of the secondary induction, plastochron decreases relatively to phyllochron, and then ρ primordia (ρ > 1) are produced per newly emerged leaf. All primordia produced before floral transition will become leaves, while primordia produced after floral transition will become spikelets. Primordia production ceased when the last leaf primordia is differentiated (for colour figure refer to online version).</figcaption></figure>

at the same time avoids ‘overparameterization’ of the model compared with our quantitative knowledge of the processes (Heide 1994). Each induction phase is represented by a variable ($ind_{apex}^{IR}$, $ind_{apex}^{II^R}$, respectively) ranging from 0, at apex initialization, to 1 when the induction requirements are completed. These variables are incremented daily according to the average temperature and photoperiod of the day, respectively. Both induction phases are irreversible and modify the apex status (Equation (4)). The consequences of the floral induction and consecutive floral transition on apex dynamics and leaf elongation are detailed below.

$$apex_n^{\text{status}} = \begin{cases} \text{vegetative} & \text{if } ind_{apex, n}^{IR} < 1 \\ \text{primary induced} & \text{if } ind_{apex, n}^{IR} = 1 \\ \text{reproductive} & \text{if } ind_{apex, n}^{IR} = 1 \end{cases} \quad (4)$$

#### 2.1.4.1 Primary induction.
Although the photoperiod could modulate the effect of low temperatures (Heide 1994), the primary induction is assumed to be dependent only on temperature in L-GrassF. In the literature, the primary induction is generally characterized by three cardinal temperatures: a minimum, optimum and maximum (He et al. 2012; Woodward et al. 2020).

The daily increment of the primary induction of an apex n ($Δind_{apex}^{IR}$, Equation (5)) linearly increases with mean air temperature from the minimum temperature ($T_{min}^{IR}$, °C) to a maximum ($Δind_{max}^{IR}$) reached for the optimal temperature $T_{opt}^{IR}$. Above that temperature, $Δind_{apex}^{IR}$

linearly declines down to zero at the maximum temperature ($T_{max}^{IR}$, °C, Fig. 1.3).

$$\Delta ind_n^{IR}(d) = \begin{cases} 0; & \text{for } T(d) < T_{min}^{IR} \text{ or } T_{max}^{IR} > T(d) \\ \Delta ind_{max}^{IR} \times \frac{T(d) - T_{min}^{IR}}{T_{opt}^{IR} - T_{min}^{IR}}; & \text{for } T_{min}^{IR} \le T(d) \le T_{opt}^{IR} \\ \Delta ind_{max}^{IR} \times \frac{T(d) - T_{max}^{IR}}{T_{opt}^{IR} - T_{max}^{IR}}; & \text{for } T_{opt}^{IR} \le T(d) \le T_{max}^{IR} \end{cases} \quad (5)$$

where $T(d)$ is the mean daily temperature (°C). The progress towards full primary induction ($ind_{apex}^{IR}$) is simulated as a time integral (Equation (6)), Day 1 being the start of the simulation.

$$ind_n^{IR} = \sum_{day=1}^{d} \Delta ind_n^{IR}, \text{ with } ind_n^{IR} \in [0, 1] \quad (6)$$

#### 2.1.4.2 Secondary induction.
The secondary induction only starts when the primary induction is complete and at least one leaf of the tiller has emerged. The latter condition reflects the need for photoperiod-sensitive organs to perceive the signal. The secondary induction increases at a constant rate above a minimal photoperiod (Equation (7)); Fig. 1.4). Also, we assumed that the maximal rate of secondary induction is only reached when a minimum number of leaves have fully emerged.

$$\Delta ind_n^{IR}(d) = \begin{cases} 0; & \text{for } PP(d) \le PP_{min} \\ \min(PPRM, PPRM \times \frac{LN_{x}}{LN_{M}}); & \text{for } PP(d) > PP_{min} \end{cases} \quad (7)$$
where PP(d) is the photoperiod of the day (h), PP_min (h) is the minimum photoperiod for secondary induction, PPRM is the maximal daily increment of the secondary induction, LN_n is the number of emerged leaves of the tiller n and LNM is the minimum number of leaves above which the secondary induction increment is maximal (PPRM).

The progress towards full secondary induction (ind_n^II^R) is simulated as a time integral Equation (8), day × being the first day after the end of the primary induction.

$$ind_n^{II^R} = \sum_{day=x}^{d} \Delta ind_n^{II^R}, \text{ with } ind_n^{II^R} \in [0, 1] \quad (8)$$

When the secondary induction is complete (ind_n^II^R = 1), the floral transition occurs, meaning that the apex switches from vegetative to reproductive development.

**2.1.5 Reproductive morphogenesis.** This section describes tiller functioning during the reproductive stages. After the floral transition of a tiller, L-GrassF assumes that (i) the rate of primordium initiation by the apex is increased, (ii) the primordia initiated by the apex progressively produce spikelets instead of leaves, (iii) leaf elongation is faster and (iv) long internodes start to elongate. The heading of the spike is approximated from the date of flag leaf ligulation and therefore results from both vegetative and reproductive developments.

**2.1.5.1 Dynamics of apex functioning.** Until the secondary induction begins, primordia are initiated at the same rate as leaf emergence and they will later develop into leaves (Fig. 2) (Hay and Kemp 1990). As soon as the secondary induction begins, the rate of primordium initiation is increased (Gillet 1980; Thomas 1980; Gonthier and Francis 1989): ρ-primordia are produced at each leaf emergence, with ρ > 1 (Table 1). This is in accordance with experimental studies which report accumulation of primordia at the apex level during the floral induction (Gonthier and Francis 1989).

All primordia initiated during the secondary induction will develop into leaves, while those initiated after the secondary induction, i.e. after the floral transition, will later develop into spikelets (Fig. 2). Consequently, the final number of leaves is determined at the time of floral transition, i.e. the end of the secondary induction. After floral transition, the initiation of primordia remains coordinated with leaf emergence and continues at the same pace, forming spikelets, until the last leaf starts to elongate. The cessation of primordium initiation marks the death of the terminal apex and no supplementary organ will then be produced.

**2.1.5.2 Leaf elongation.** There has been a debate in the literature about whether the increase in leaf elongation, frequently observed in spring in temperate climates, was directly caused by floral transition or whether these two processes were affected independently by common factors, leading to synchrony in most growth conditions (Davies 1971; Parsons and Robson 1980). Davies (1971) finally concluded that an acceleration of leaf elongation could be observed in non-induced plants, which leads to the conclusion that the acceleration of leaf elongation was not only caused by floral transition.

In the model, the increase of leaf elongation rate is triggered after the completion of a dual induction described by the same equations

as for the floral induction (Equations (5–8)). The thresholds to be reached for each phase completion are given by two parameters (index_n^II^R, index_n^II^R). For the sake of simplicity, these indices were set to 1 in the following analyses, assuming exact synchrony between leaf growth increase and floral transition in the explored environments. The increase in leaf elongation rate is implemented by using a multiplicative factor, kY_m^0, applied on Y_m^0 (Equation (2B)), which affects the rate of leaf elongation, the duration of the hidden phase and final leaf length ultimately leaf elongation rate. As an example, kY_m^0 < 1 increases the initial rate of leaf elongation and then decreases the duration of the hidden phase (earlier emergence). Consequently, kY_m^0 < 1 will decrease final leaf length and therefore increase the rate of leaf elongation outside of the pseudo-stem.

**2.1.5.3 Internode elongation.** The floral transition has been shown to trigger the elongation of internodes (Gillet 1980). However, the number of long internodes, the kinetics of their elongation and their final dimensions are poorly described in the literature for perennial grasses. Although internode extension dramatically increases plant height and apex position (Hazard et al. 2006), it does not affect the HD, at least geometrically, as the distance between the spike tip and the highest ligule is not affected by the extension of the internode located below the apex (Rouet et al. 2021). As a consequence, internode elongation is extremely simplified in L-GrassF, with the main objective to provide realistic 3D representations of reproductive tillers. Coordination rules between internode and leaf elongation were used as proposed for wheat and maize (Vidal and Andrieu 2009; Gauthier et al. 2020): internode n elongation requires the flower transition to be induced and starts at the ligulation of leaf n. Internode length increases linearly until its final length equals that of sheath n.

**2.1.5.4 Heading of the spike.** The heading of a tiller is defined as the emergence of the terminal spikelet from the pseudo-stem (Gillet 1980; Thomas 1980). In a majority of cases, the pseudo-stem at this time corresponds to the length of the last sheath of the tiller. Therefore, the date of the terminal spikelet emergence depends on the length of the flag leaf sheath but also on the elongation of the peduncle bearing the spike, the number of spikelets and the elongation of the short internodes which bear each of the spikelets. As for internodes, the literature on peduncle and spike elongation is very sparse. Instead of simulating the heading of the spike from geometrical calculations involving uncertain elongation rules of the different organs cited above, we preferred to use the flag leaf ligulation as a proxy for heading.

## 2.2 Experimental data set

The model was evaluated using a data set which describes the HD of several cultivars of L. perenne for several locations and years. The data set was obtained by GEVES, the French Variety and Seed Study and Control Group, which is responsible for the evaluation of candidate cultivars for their insertion in the French catalogue of cultivars. The assessment criteria of the cultivars include the HD, measured by GEVES as the first date when at least 10 spikes are visible per linear metre on row trials. The measurement consists of a visual assessment of the number of headed tillers every 3 days in six different sites in France. The HD of a cultivar corresponds to the average of three
replicates. The sward management follows the National protocol (GEVES 2006) repeated in all locations and years. Seeds are sown in spring in common gardens, in 6-m-long rows (5 g per linear metre) with three repetitions per cultivar. The rows are mowed during the first summer and the early winter. Plants are fertilized with N during their whole development to avoid nutritive stress.

In the present study, we selected seven forage cultivars with contrasting HD and used as control varieties by GEVES (Bronsyn, Indiana, Bargala, Lactal, Milca, Carillon and Escal). The detailed characteristics of these cultivars and indicative values for the variables of interest are given in Table 2. The data set covers a period of 17 years (2001–17) in six French locations (Table 3; see Supporting Information—Fig. S2). However, not all combinations of cultivar × location × year are available in the data set. Because the model does not account for drought effects, the data set was restricted to environmental conditions where water resources before heading were not limiting. Environments where the water shortage would have potentially affected plant growth were identified using a water balance approach as follow. The soil relative water content (SRWC) was considered maximum (value of 1) on 1 December 2000 for all sites and updated daily considering rainfalls and evapotranspiration (after Kunrath et al. 2015). In accordance with FAO (Allen et al. 1994), each day when SRWC was lower than 0.4 was considered detrimental for the development of *L. perenne*. Finally, a site-year combination was removed from the data set when the period between the start of regrowth in spring and the HD (Table 2) included at least 25 % of days with SRWC < 0.4.

## 2.3 Model evaluation

Given the available data, the calibration and evaluation of the model were performed for main stems only; the proportion of reproductive tillers was therefore not assessed in this study. As a first step, a sensitivity analysis was performed to identify the parameters affecting the HD of the main stems the most. Then, these parameters were calibrated for each of the seven cultivars using a subset of the data set described above. Finally, the model was evaluated for its ability to simulate the genetic diversity of the HD across all locations and years using a second subset of the data set (Tables 2 and 3; see Supporting Information—Figs S1 and S2).

### 2.3.1 Model inputs and initial conditions of simulations.
L-GrassF was used with an hourly time step to simulate the HD of seven *L. perenne* cultivars over 14 years and six locations (situations without water shortage). To date, very little is known on the HD variability between tillers of the same plant and this aspect of the model has not been assessed yet. Therefore, the 10th tiller to head in each experimental trial was simulated in L-GrassF. Virtual sowing was carried out on 15 September and a cut at 5 cm height was performed on 15 November in all locations. The virtual sowing was set later than the actual sowing to reduce the simulation time, assuming that this would not affect the simulation of HD. Mean daily temperatures from spring to 15 September were higher than the maximum temperature allowing for the primary induction (12 °C). So simulating plant development during that period would not have modified the floral status of the apices. Concerning the vegetative development, the 2 months between the sowing and the cutting date in our simulations allowed the tillers to reach a permanent

<table><caption>Table 2. Characteristics of the cultivars of *Lolium perenne* used for model evaluation. The trials were conducted by the GEVES institute and measurements were carried out on lines. Registration year is the year of the cultivar registration in the French catalogue. Start of regrowth in spring is the first date from 1 January when the height of stretched leaves exceeds 20 cm. Heading start is the first date from 1 January when at least 10 spikes are visible per linear metre.</caption><thead><tr><th>Cultivar</th><th>Breeder</th><th>Country</th><th>Registration year</th><th>Ploidy</th><th>Start of regrowth in spring (doy)</th><th>Heading start (doy)</th><th>Number of observations</th></tr></thead><tbody><tr><td>Bronsyn</td><td>Agriseeds Ltd</td><td>NZ</td><td>2001</td><td>2n</td><td>68–98</td><td>124</td><td>35</td></tr><tr><td>Indiana</td><td>Force Limagrain SA/DFL Trifolium</td><td>FR/DK</td><td>2004</td><td>2n</td><td>82–98</td><td>135</td><td>36</td></tr><tr><td>Bargala</td><td>Barenbrug Holland BV</td><td>NL</td><td>2000</td><td>4n</td><td>82–92</td><td>135</td><td>32</td></tr><tr><td>Lactal</td><td>RAGT/R 2n</td><td>FR</td><td>2004</td><td>4n</td><td>92–105</td><td>142</td><td>26</td></tr><tr><td>Milca</td><td>Carneau Frères SA</td><td>FR</td><td>1996</td><td>2n</td><td>92–105</td><td>145</td><td>27</td></tr><tr><td>Carillon</td><td>Carneau Frères SA</td><td>FR</td><td>2002</td><td>2n</td><td>92–112</td><td>159</td><td>26</td></tr><tr><td>Escal</td><td>R 2n</td><td>FR</td><td>2005</td><td>4n</td><td>—</td><td>—</td><td>23</td></tr></tbody></table>

DK: Denmark; FR, France; NL, Netherlands; NZ, New Zealand.
<table><caption>Table 3. Growth conditions of the seven *Lolium perenne* cultivars used for model evaluation. Temperatures and rainfall are means calculated for the period from 2001 to 2017.</caption><thead><tr><th>Location</th><th>Latitude</th><th>Longitude</th><th>Altitude (m)</th><th>Mean winter temperature (°C)</th><th>Mean spring temperature (°C)</th><th>Photoperiod at summer solstice (h)</th><th>Soil available capacity (cm)</th><th>Annual rainfall (mm)</th></tr></thead><tbody><tr><td>Exmes</td><td>48.759°</td><td>0.126°</td><td>205</td><td>4.5</td><td>9.8</td><td>16.2</td><td>50</td><td>755</td></tr><tr><td>Guyancourt</td><td>48.771°</td><td>2.097°</td><td>161</td><td>4.7</td><td>10.9</td><td>16.2</td><td>190</td><td>577</td></tr><tr><td>Lusignan</td><td>46.418°</td><td>0.119°</td><td>150</td><td>5.1</td><td>11.1</td><td>15.8</td><td>166</td><td>787</td></tr><tr><td>Ploudaniel</td><td>48.504°</td><td>-4.328°</td><td>85</td><td>7.2</td><td>10.4</td><td>16.1</td><td>161</td><td>1152</td></tr><tr><td>Theix</td><td>45.723°</td><td>3.018°</td><td>830</td><td>1.2</td><td>6.9</td><td>16.0</td><td>50</td><td>800</td></tr></tbody></table>

leaf production regime. Cutting the plants on 15 November as in real practice reduced the pseudo-stem length to 5 cm, leading to plant architecture close to that of the data set. Meteorological data used for the simulations were collected from the INRAE weather network (data available in the CLIMATIK database; https://intranet.inrae.fr/climatik_v2) or from the 8-km gridded Météo France 'SAFRAN' database when the test site was too far from an INRAE weather station. From these databases, we extracted mean daily temperatures (which affect plant development and primary induction) as well as daily rainfall and potential evapotranspiration to exclude drought situations from the data set. The daily photoperiod was calculated by using the Python package 'ephem' (https://pypi.org/project/ephem/). L-GrassF was run in a high-performance computing centre (MESO@LR, Université de Montpellier, France).

**2.3.2 Sensitivity analysis.** A sensitivity analysis was performed on seven parameters selected for their role in leaf elongation ($Y_m^0 kY_m^0$), apex dynamics ($\rho$) and floral induction ($\Delta ind_{max}^R PP_{min}^R$ PPRM, LNM). To optimize the number of simulations of the sensitivity analysis, we used a screening method based on a fractional factorial design (Morris 1991). Each of the seven parameters studied had six possible values around a value obtained by preliminary tests (not shown). The method uses a set of model runs called trajectories, each composed of successive simulations which differs by a unique parameter value (One at A Time—OAT), and each parameter varying once in a trajectory. In our case, a trajectory corresponded to 8 simulations (an initial random run and 7 simulations corresponding to the successive variations of the 7 parameters). Overall, 40 trajectories were performed, which made 320 simulations instead of 279 936 simulations ($6^7$) with a complete factorial design.

The effect of each parameter on the HD was estimated by computing the main effect indices ($\mu_i^*$) and the sensitivity indices for second-order interactions ($\sigma_i$), calculated using Equation (9).

$$ 
\begin{aligned}
\Delta_i^r &= \frac{HD(T_{i-1}^r) - HD(T_i^r)}{\delta} \\
\mu_i &= \frac{1}{R} \sum_{r=1}^{R} \Delta_i^r \\
\mu_i^* &= \frac{1}{R} \sum_{r=1}^{R} |\Delta_i^r| \\
\sigma_i &= \sqrt{\frac{1}{R-1} \sum_{r=1}^{R} (\Delta_i^r - \mu_i)^2}
\end{aligned}
\qquad (9) $$

where HD is the HD predicted by L-GrassF, $T^r$ is the trajectory r, R is the number of trajectories and $\delta$ is the step between two levels of parameter i.

The same sensitivity analysis was performed for three contrasted environmental conditions to test the independence between the results of the sensitivity analysis and the environmental conditions. Since the photoperiod was not significantly different between the sites of the data set, the three environments were chosen to account for temperature variability: a cold environment (Theix 2009–10), a temperate environment (Ploudaniel 2009–10) and a warmer environment
(Ploudaniel 2015–16). The daily temperature variations in these environments are given in **Supporting Information—Fig. S3**.

**2.3.3 Calibration and validation of the model.** The objective was to identify sets of parameters which would be able to simulate the HD of seven cultivars in any occurrence of the experimental data set. The data set was divided randomly into two independent subsets: 60 % (n = 120) of the data set was used for calibration and the remaining 40 % (n = 85) for model validation. The calibration was done only on the three most sensitive parameters revealed by the sensitivity analysis (Y<sub>m</sub><sup>0</sup> kY<sub>m</sub><sup>0</sup> and PPRM). The other parameters were set to the mean value of the range used in the sensitivity analysis (Table 4). A complete factorial plan with six levels for Y<sub>m</sub><sup>0</sup> kY<sub>m</sub><sup>0</sup> and PPRM was set up, leading to 216 scenarios (6<sup>3</sup>) which were run in each site-year combination. For each cultivar, the root-mean-square deviation (RMSE) was calculated between the observed and the simulated HD for each combination of parameters and environments. The final set of parameters selected for each cultivar was the one which minimized the RMSE on HD. By mapping RMSE values in three-dimensional space (Y<sub>m</sub><sup>0</sup> kY<sub>m</sub><sup>0</sup> and PPRM), we checked if this RMSE value corresponded to the global minimum of the space explored.

The final sets of parameters identified for each cultivar were used for model validation. We evaluated the ability of the model to simulate the observed HD by calculating the RMSE and the R² from the second subset of the data set.

**2.3.4 Statistical analyses.** Analyses of variance and non-linear regression were performed with R software (R Core Team 2021). Analyses of variance (ANOVAs) were performed following a one-factor linear model as Equation (10).

$$y_{ij} = \mu + C_i + \varepsilon_i \quad (10)$$

where y is any dependent variable, μ is the mean value of y, Cᵢ is the effect of cultivar i and ε is the random error.

Normal distributions of the residuals of ANOVAs were tested using the Shapiro–Wilk test. Homoscedasticity was checked by random distribution of the residuals.

The root-mean-square deviation (RMSE) was calculated as Equation (11).

$$RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (S_i - O_i)^2} \quad (11)$$

where n is the number of individuals, Sᵢ and Oᵢ are, respectively, the simulated and observed value of individual i for the HD.

### 3. RESULTS

The geometrical interpretation of the L-system provided 3D views of the plant architecture at different stages of development (Fig. 3A). L-GrassF also produced original insights in the dynamics of the floral induction and leaf ligation in relation to the environmental conditions (Fig. 3B). Simulating floral transition allows us to estimate the final number of leaves of the tiller, as well as tiller demography since reproductive tiller dies after flowering. Finally, the integration of the interactions between

<table><caption>Table 4. Model parameters used for the sensitivity analysis and calibrated for each cultivar.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Range of values for sensitivity analysis</th><th>Range of values for calibration</th><th>Unit</th></tr></thead><tbody><tr><td>Δind<sup>R</sup><sub>max</sub></td><td>Maximum daily rate of the primary induction</td><td>0.015–0.07</td><td>0.0425</td><td>Dimensionless</td></tr><tr><td>PP<sub>min</sub></td><td>Minimum photoperiod for secondary induction</td><td>10–12</td><td>11</td><td>h</td></tr><tr><td>LNM</td><td>Minimum number of leaves for secondary induction at maximum daily rate</td><td>15–25</td><td>20</td><td>Dimensionless</td></tr><tr><td>ρ</td><td>Number of primordia initiated at leaf emergence</td><td>2–7</td><td>4.5</td><td>Primordium per leaf</td></tr><tr><td>Y<sup>0</sup><sub>m</sub></td><td>Maximum leaf length before emergence</td><td>80–120</td><td>80–120</td><td>mm</td></tr><tr><td>kY<sup>0</sup><sub>m</sub></td><td>Multiplicative factor of Y<sup>0</sup><sub>m</sub></td><td>0.05–1.5</td><td>0.05–1.5</td><td>Dimensionless</td></tr><tr><td>PPRM</td><td>Maximum daily rate of the secondary induction</td><td>0.1–1</td><td>0.1–1</td><td>Dimensionless</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Figure 3. Dynamic 3D representations of plant architecture and numeric outputs from L-GrassF. The model was run in the environmental conditions of Lusignan, with sowing on 15 September 2010, cutting at 5 cm on 15 November 2010. In this simulation, tillering was activated for illustration only as the present study focused on the main tiller. (A) Snapshots of the 3D representation of four development stages of a virtual *Lolium perenne* (1, 2, 3 and 4). Green structures are leaves, red structures are internodes and purple structures are spikelets. (B) Daily numeric outputs of the model for the main tiller. Stages 1, 2, 3 and 4 are the same as in (A). Blue dashed line is the primary induction variable *ind*<sup>I1</sup> and yellow dashed line is the secondary induction variable *ind*<sup>I2</sup>. Light green line accounts for the cumulated number of ligulated leaves from tiller appearance, HD being consider as occurring at last leaf ligulation. Dark green line represents the cumulated length of all leaves (sheath and lamina) from tiller appearance, without senescence (for colour figure refer to online version).</figcaption></figure>
the vegetative and reproductive developments provided estimates of the
leaf area which is harvestable at each date according to the temperature,
photoperiod and stages of development (Fig. 3B).

<figure><img src="image_4.png" /><figcaption>Figure 4. Morris sensitivity indices for the HD. Mean effects (μ*) and standard deviation of the elementary effects (σ) of the parameters studied in the three sensitivity analyses, corresponding to three environmental conditions (blue dots: Theix 2009–10, orange triangles: Ploudaniel 2009–10, red diamonds: Ploudaniel 2015–16) (for colour figure refer to online version).</figcaption></figure>

**3.1 Sensitivity analysis: the predominant effect of leaf growth and photoperiodic-related parameters**

Three main groups of parameters were identified through the sensi-
tivity analysis (Fig. 4): $Y_m^0$ and $kY_m^0$ were highly sensitive parameters,
PPRM, $PP_{min}$ and $\rho$ had a moderate effect on HD while LNM and
$\Delta ind_{max}^{IR}$ little affected the HD. This pattern was similar across the
three tested environments.

The parameters $Y_m^0$ and $kY_m^0$, which are involved in the determina-
tion of final leaf length, presented a strong principal effect of, respec-
tively, 52 and 66 days on HD, whatever the environment. $kY_m^0$ positively
affected the HD (i.e. later heading) for all the values in the tested range
(0.05–1.5). $Y_m^0$ positively affected the HD for values ranging between
80 and 96 and showed a plateau beyond this range (Fig. 5). Second-
order interactions were also significant for the two parameters, reflect-
ing a crucial interaction between the parameters ($\sigma Y_m^0 \approx 34$ days
and $\sigma kY_m^0 \approx 21$ days). The maximal rate of the secondary induction
(PPRM) had a smaller mean effect than $Y_m^0$ and $kY_m^0$ on the HD but this
effect was still significant in the three environments, ranging from 23 to
32 days (Fig. 5). However, the effects of PPRM were mainly observed
for a narrow range of values (between 0.1 and 0.28; Fig. 5). In addi-
tion, the HD occurred earlier with increasing PPRM values since this
parameter controls the rate of the secondary induction. The parameter
$PP_{min}$ appeared to have a strong mean effect ($\mu^* \approx 25$ days) in all envi-
ronments. The mean primary effect of the maximum primary induc-
tion rate ($\Delta ind_{max}^{IR}$) on HD was low ($\approx 10$ days). The sensitivity of
$\Delta ind_{max}^{IR}$ increased with increasing mean temperature, meaning that
in warm environments, the rate of the primary induction has a stronger
impact on the HD than in cold environments, where the primary

<figure><img src="image_5.png" /><figcaption>Figure 5. Relations between HD and the three more sensitive parameters $Y_m^0$ and $PPRM$ HD values presented in the figure are outputs of the sensitivity analysis. For each of the three most sensitive parameters, all simulations of the sensitivity analysis are considered. Therefore, not all combinations of parameters have been studied. Tukey's style: -1.5 * interquartile range (IQR), first quartile, median, third quartile, +1.5 * IQR.</figcaption></figure>
induction is very soon regardless of $\Delta ind_{max}^{IR}$. In the three sensitivity analyses, the main effect of LNM was low, which can be explained by the fact that most of the plants had developed a large number of leaves at the completion of the primary induction, which led to a maximal rate of secondary induction whatever the values of LNM. This suggests a weak effect of this parameter in natural conditions where sowing is often carried out long before the flowering period, while this parameter could be useful to simulate experimental results obtained in artificial conditions. In the following, our study is limited to the three parameters that affected the HD the most: $Y_m^0$ and $kY_m^0$ and PPRM.

## 3.2 Comparison of simulated and observed HD for contrasted environments and cultivars

The analysis of the data set revealed that the cultivars presented a wide range of HD values between 2001 and 2017, with over 1 month between the earliest and the latest ones (Fig. 6). A considerable variability of the HD was also observed between sites and years, as shown by the boxplots (SD = 23 days for Escal to 36 days for Indiana). Despite the environmental variability, cultivar ranking according to their earliness was similar regardless of the environmental conditions (results not shown). Therefore, predicting the HD in *L. perenne* requires consideration of the genetic effect (supported by the cultivar), the environmental effect and their interactions (Keep *et al.* 2020).

The calibration of the model allowed us to identify a combination of parameters for each of the seven cultivars of the data set (Table 4). Using the optimal set of parameters for each cultivar, L-GrassF satisfactorily simulated the diversity of HD of the training set (Fig. 7, grey dots). Overall, the RMSE was lower than 8 days, but the coefficients of determination obtained from the linear regressions between simulated and observed HD were low (<0.60) (Table 5). This is due to the small amount of data available and the low temporal range of HD observed within cultivars. For cultivars Escal and Carillon, we found the same set of parameters, which is in accordance with their close observed HD (Fig. 7).

<figure><img src="image_6.png" /><figcaption>Figure 6. Diversity of the HD of seven cultivars of *Lolium perenne* as observed in GEVES network for six locations [see Supporting Information—Fig. S2] and 14 years. In this plot, HD corresponds to the definition used by GEVES, i.e. the first date when at least 10 spikes are visible per linear metre of trial lines. Number of observations were Bronsyn: 35, Indiana: 36, Bargala: 31, Lactal: 26, Milca: 27, Carillon: 26, Escal: 23. Cultivars were not represented in all environments. Anova with model $Y_i = \mu + C_i + \varepsilon_i; F: 63.918, P < 0.05$.</figcaption></figure>

The validation of L-GrassF on the independent data set showed variable results depending on the cultivar (Fig. 7). In Bronsyn, Indiana, Lactal, Carillon and Escal, the model provided good overall estimates of the HD (RMSE < 10 days) but the coefficients of determination were low. This could be explained by the narrow range of observed HD for these cultivars, especially for Carillon. For Milca, the elevated RMSE was partly explained by one environment (Lusignan 2014–15) for which the model clearly overestimated the HD (172 instead of 137) (Fig. 7E). Calculating the RMSE without taking this into account reduced it from 12.17 to 5.72 days (–53 %). For Bargala, the RMSE was higher (16.67 days) and the determination coefficient was very close to zero (0.08), meaning that the model was not able to give better prediction than the mean observed HD (Fig. 7D). Without taking into account the environment (Ploudaniel 2014–15) where the model clearly overestimated the HD (173 instead of 128) reduced it from 16.67 to 10.79 days (–35 %). We did not find any clear explanation for the model's behaviour in these situations which were not characterized by unusual environmental conditions. Considering all cultivars together highlighted an overall efficient prediction of the HD (RMSE ~ 10 days) and a good representation of cultivar earliness ($R^2$ = 0.48).

## 4. DISCUSSION

The environmental factors responsible for the floral induction, mainly temperature and photoperiod, have been well identified in perennial grasses (Heide 1994). Despite their importance for the success of reproductive development and the consequences on the functioning and management of grasslands, attempts to mechanistically predict their effects on reproductive phenology are sparse (Rouet *et al.* 2021). In this study, we developed a sequential two-phase model of floral induction, as proposed for *L. perenne* by Heide (1994). This kind of representation is similar to the propositions made by Woodward *et al.* (2020) for forage grasses and to the model Sirius developed for wheat (Brooking 1996; Jamieson *et al.* 1998; He *et al.* 2012; Brown *et al.* 2013). Unlike the latter models, the floral induction in L-GrassF occurs independently for each tiller, although this aspect was not evaluated in the present study. Based on the literature, we assumed that in *L. perenne*, the primary induction is mainly dependent on temperature and the secondary induction on photoperiod (Heide 1994). Nevertheless, some interactions between these two factors were observed during the induction process (Heide 1994; Aamlid *et al.* 2000). For primary induction, short photoperiods are known to increase the effect of low temperature. For secondary induction, high temperatures increase the effect of photoperiod. Some authors even suggested that specific secondary induction conditions could compensate for incomplete primary induction (Cooper 1960; Heide 1994). If the qualitative effects of the temperature, the photoperiod and their interactions during each phase are well known, there is a lack of quantitative assessments of floral induction response to temperature and photoperiod, in particular in natural conditions. Furthermore, in all environmental conditions of our data sets, periods of low temperature were always concomitant with periods of short photoperiod, respectively, long photoperiod and high temperatures. For primary induction, the cardinal temperatures (Equation (5)) were identical for all cultivars and only the maximum induction rate ($\Delta ind_{max}^{IR}$) was considered variable. This hypothesis was supported by studies on floral induction (Aamlid *et al.* 2000; Paina
<figure><img src="image_7.png" /><figcaption>Figure 7. Comparison of observed and predicted values of HD obtained with L-GrassF for seven cultivars. (A) Bronsyn, (B) Indiana, (C) Lactal, (D) Bargala, (E) Milca and (F) Carillon, (G) Escal. Empty and filled circle represent the sub-data sets used for model calibration and validation, respectively.</figcaption></figure>

<table><caption>Table 5. Values of parameters identified by calibration. Parameters $\Delta ind_{max}^{I^R}$, $PP_{min}$, LNM and $\rho$ were calibrated with the sensitivity analysis.</caption><thead><tr><th>Cultivar</th><th>Y<sub>m</sub><sup>0</sup></th><th>kY<sub>m</sub><sup>0</sup></th><th>PPRM</th><th>Δind<sub>max</sub><sup>I<sup>R</sup></sup></th><th>PP<sub>min</sub></th><th>LNM</th><th>ρ</th><th>RMSE</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>Bronsyn</td><td>112</td><td>0.05</td><td>0.1</td><td>0.0425</td><td>11</td><td>20</td><td>4.5</td><td>6.67</td><td>0.57</td></tr><tr><td>Indiana</td><td>104</td><td>0.34</td><td>0.1</td><td>0.0425</td><td>11</td><td>20</td><td>4.5</td><td>6.28</td><td>0.57</td></tr><tr><td>Lactal</td><td>96</td><td>0.63</td><td>0.1</td><td>0.0425</td><td>11</td><td>20</td><td>4.5</td><td>7.19</td><td>0.32</td></tr><tr><td>Bargala</td><td>120</td><td>0.34</td><td>0.82</td><td>0.0425</td><td>11</td><td>20</td><td>4.5</td><td>7.20</td><td>0.49</td></tr><tr><td>Milca</td><td>112</td><td>1.21</td><td>0.82</td><td>0.0425</td><td>11</td><td>20</td><td>4.5</td><td>5.75</td><td>0.45</td></tr><tr><td>Carillon</td><td>104</td><td>0.63</td><td>0.1</td><td>0.0425</td><td>11</td><td>20</td><td>4.5</td><td>6.20</td><td>0.24</td></tr><tr><td>Escal</td><td>112</td><td>0.63</td><td>0.1</td><td>0.0425</td><td>11</td><td>20</td><td>4.5</td><td>5.03</td><td>0.51</td></tr></tbody></table>

et al. 2014) which showed that cultivars of contrasting geographical origins could be induced by exposure to the same temperature but that the exposure duration was critical. For secondary induction, a threshold photoperiod is generally observed but with different duration requirements (Heide 1994; Aamlid et al. 2000). This justifies our choice to have secondary induction controlled by two parameters, a threshold ($PP_{min}$) and an induction rate ($PPRM$).

## 4.1 Importance of morphogenesis in HD prediction

The main originality of L-GrassF is to account for the interactions between floral induction processes and morphogenesis. The sensitivity analysis of the model revealed that two of the morphogenetic parameters ($Y_m^0$ and $kY_m^0$) had the greatest influence on the HD simulation. $Y_m^0$ and $kY_m^0$ are parameters affecting leaf elongation rate and the rhythm of leaf emergence. High values of $Y_m^0$ and $kY_m^0$ produce long leaves with long sheaths. Because of the self-regulation rules

implemented in the model, this increases the time spent by leaves inside the pseudo-stem, which delays the initiation of the next leaves. Therefore, for a given number of leaves to produce before heading, the expanding of all leaves will take longer with high values of $Y_m^0$ and $kY_m^0$ and HD will be delayed. Delayed HD for cultivars producing long leaves is in opposition with the observations made by Hazard (2006) on short and long leaf morphotypes. This may be due to different apex dynamics resulting in different final number of leaves. Further investigations with monitoring of the number of leaves would be required to establish the relationship between the final leaf length and HD. The parameter $Y_m^0$ played a central role in the first stages of leaf morphogenesis, similarly to that in the initial version of L-Grass, controlling leaf elongation rate before leaf emergence (Verdenal et al. 2008). With the objective to predict the HD of unknown cultivars, $Y_m^0$ and $kY_m^0$ are parameters of interest to be calibrated, as they can easily be estimated by measuring the leaf profiles (final lengths of laminas and sheaths)
in winter for $Y_m^0$ and in late spring for $kY_m^0$. Given the importance of morphogenesis on heading, it appears necessary to integrate the main factors controlling morphogenesis in L-GrassF. In particular, recent progress in the field of FSPMs has been made on the interactions between leaf morphogenesis and the trophic and hydraulic status of plants (Coussement et al. 2018; Albasha et al. 2019; Gauthier et al. 2020).

## 4.2 L-GrassF satisfactorily predicted HD despite many unknowns on the data set

L-GrassF was calibrated and validated in a range of environments which correspond to those of the North-West of France. Therefore, the range of latitudes explored by the data set did not provide any significant differences in photoperiod between the sites. The climatic contrast between the trials thus lays mainly in temperature range and distribution. In addition, our data set was limited to trials that had not been exposed to water stress or lack of mineral elements. Consequently, the model validity is only ensured for environmental conditions similar to the study. As shown by Blanco-Pastor et al. (2019), the genetic diversity used for the selection of new cultivars was mainly derived from natural ecotypes originating from North-West Europe only. The area fully encompasses the extension of our study area. Therefore, it would be of interest to verify that the model is also able to accurately simulate the behaviour of genotypes from geographic areas with very different climates and photoperiodic conditions like in the Northern Europe, where HD is delayed after snow melting, or the Mediterranean area, where HD occurs early to avoid drought (Heide 1994; Barre et al. 2017; Rouet et al. 2021).

Under the trial conditions used for this study, the HD was always observed, meaning that floral induction was always reached by at least 10 tillers per linear metre (see Materials and Methods). As the absence of floral induction was not observed, our data set did not allow us to discover what the behaviour of the model would be on the same cultivars grown under environmental conditions which would not allow floral transition, for instance with winter temperatures too high to complete primary induction. Owing to an explicit implementation of the floral induction, our model should be able to simulate plant behaviours in these critical situations. Nevertheless, it would be of interest to conduct trials in more contrasted sites.

Although all trials of our data set were theoretically managed according to a common protocol, there was significant uncertainty about the practical application of this protocol for each year and location. The sowing dates in spring may have been different according to the local climatic conditions. For instance, it may have been delayed if the soil was not dry enough to allow access to farm machinery. By initializing L-GrassF just before the first autumn (instead of the previous spring), we considered that the period preceding that date had no effect on the floral induction, in particular because during that period, temperatures were too high to induce the primary induction and because the development only affected the size of the vegetative organs. In order to homogenize organ sizes, a virtual cutting was made on the tillers on 15 November for all conditions. The cutting was analogous to the regularization cutting generally made by GEVES just before the winter (no precise date available). In practice, GEVES

cutting was however unevenly made between the trials depending on plant height and weather conditions.

Despite all these unknowns, L-GrassF was overall able to simulate the genetic and environmental variability in HD between cultivars, years and locations with a precision of circa 9 days (all locations, years and cultivars together). However, due to the lack of measurements at tiller scale, we could not compare all output variables produced by the model (final leaf number, leaf emergence date, spikelet number). With such variables available, an improvement of the prediction accuracy could be expected. In comparison with existing empiric models, the main gain of the calibration of L-GrassF is expected to cover a wide range of conditions and to be independent from the calibration data set.

## 4.3 Perspective of use

From this study, L-GrassF has shown interesting performances in predicting the HD of *L. perenne* cultivars. In the current climatic context, where interannual temperature variations are limited, the use of simple prediction models (function of thermal time or indicative dates provided by GEVES) already provides a sufficiently accurate prediction of the HD. The current version of L-GrassF is a heuristic tool which can be used for testing ecophysiological hypotheses (synchronization between tillers, implementation of response laws...) and designing original grassland agrosystems. As a first step, we presented here the ability of L-GrassF to simulate the HD for contrasted environments and cultivars which is a prerequisite before using more aspects of the model.

The individual representation of tillers in L-GrassF is particularly useful for studying their life cycles and the consequences on tiller demography. First, it would help to predict the variations observed in tiller demography (Matthew et al. 2000) and discriminate tillers dying due to their floral development from tillers dying due to competition with others tillers (Ong et al. 1978; Colvill and Marshall 1984; Matthew et al. 2000). It could also possibly be used to understand the determinisms of aftermath heading, i.e. a second or third period of tiller heading occurring later in the season (mid-summer to autumn) and after cutting of the first period of heading. Indeed, this phenomenon is generally associated with late heading cultivars but the determinants remain unknown (Sampoux et al. 2011; Arojju et al. 2016). In L-GrassF, the simulation of each individual tiller and its floral induction coupled to the possibility of simulating cutting give all the keys to test hypotheses about the tillers concerned by aftermath heading and the links to environmental conditions.

In comparison to other existing grassland models, L-GrassF simulates a detailed structure of perennial grasses, including the number of spikelets. The possibility of simulating the quantity and timing of spikelets and grain production (Byrne et al. 2009) as a function of the environmental conditions is a significant step towards coupling between ecophysiological and genetic models. Therefore, coupling an advanced version of L-GrassF (including flower, pollen and seed production) with a model of genetic crossing would allow simulation of the evolution of the genetic structure of the population in grasslands.
## SUPPORTING INFORMATION

The following additional information is available in the online version of this article—

**Figure S1.** Increment of the physiological age in relation to temperature.

**Figure S2.** Location of the GEVES study sites in France. (1) Erdre-en-Anjou, (2) Exmes, (3) Guyancourt, (4) Lusignan, (5) Ploudaniel and (6) Theix

**Figure S3.** Daily temperatures and mean photoperiod, between 15 September and 30 July, in the three environmental conditions used for the sensitivity analysis (blue: Theix 2009–10, orange: Ploudaniel 2009–10, red: Ploudaniel 2015–16, black: photoperiod in all sites).

## ACKNOWLEDGEMENTS

We would like to thank Jade Debourbe for the analyses on the climatic data sets and the preliminary analyses of the GEVES data set (internship founded by SOERE TEMPO). This work has been carried out with the support of MESO@ LR-Platform at the University of Montpellier.

## SOURCES OF FUNDING

S.R.'s PhD was funded by a convention between the Region Nouvelle-Aquitaine and INRAe (Convention n° 2017-1R20201 - 00013154 - Thèse).

## CONTRIBUTIONS BY THE AUTHORS

Conceptualization, S.R., R.B. and J.-L.D.; methodology, S.R., R.B., D.C. and A.E.-G.; software, S.R., R.B. and J.-L.D.; data collection, D.L. and M.-H.B., data analyses and visualization, S.R., R.B. and J.-L.D.; validation, S.R., R.B. and J.-L.D.; writing—original draft preparation, S.R., R.B. and J.-L.D.; writing—review and editing, S.R., R.B. and J.-L.D.; supervision, J.-L.D.; project administration, R.B. and J.-L.D.

## CONFLICT OF INTEREST

None declared.

## LITERATURE CITED

Aamlid TS, Heide OM, Boelt B. 2000. Primary and secondary induction requirements for flowering of contrasting European varieties of *Lolium perenne*. *Annals of Botany* **86**:1087–1095.

Albasha R, Fournier C, Pradal C, Chelle M, Prieto JA, Louarn G, Simonneau T, Lebon E. 2019. HydroShoot: a functional-structural plant model for simulating hydraulic structure, gas and energy exchange dynamics of complex plant canopies under water deficit—application to grapevine (*Vitis vinifera*). *in silico Plants* **1**:diz007; doi: **10.1093/insilicoplants/diz007**.

Allen R, Smith M, Perrier A, Pereira L. 1994. An Update For The Definition Of Reference Evapotranspiration. *ICID Bulletin of the International Commission on Irrigation and Drainage* **43**:1–35.

Arojju SK, Barth S, Milbourne D, Conaghan P, Velmurugan J, Hodkinson TR, Byrne SL. 2016. Markers associated with heading and aftermath heading in perennial ryegrass full-sib families. *BMC Plant Biology* **16**:160.

Barre P, Ruttink T, Muylle H, Lootens P, Sampoux J-P, Rohde A, Combes D, Roldán-Ruiz I. 2018. Natural diversity in vegetative

and reproductive investments of perennial ryegrass is shaped by the climate at the place of origin. *Grass and Forage Science* **73**:193–205.

Barrett PD, Laidlaw AS, Mayne CS. 2005. GrazeGro: a European herbage growth model to predict pasture production in perennial ryegrass swards for decision support. *European Journal of Agronomy* **23**:37–56.

Blanco-Pastor JL, Manel S, Barre P, Roschanski AM, Willner E, Dehmer KJ, Hegarty M, Muylle H, Ruttink T, Roldán-Ruiz I, Ledauphin T, Escobar-Gutiérrez A, Sampoux J-P. 2019. Pleistocene climate changes, and not agricultural spread, accounts for range expansion and admixture in the dominant grassland species *Lolium perenne* L. *Journal of Biogeography* **46**:1451–1465.

Bonesmo H, Bélanger G. 2002. Timothy yield and nutritive value by the CATIMO model: I. Growth and nitrogen. *Agronomy Journal* **94**:337–345.

Bonhomme R. 2000. Bases and limits to using 'degree.day' units. *European Journal of Agronomy* **13**:1–10.

Boudon F, Pradal C, Cokelaer T, Prusinkiewicz P, Godin C. 2012. L-Py: an L-system simulation framework for modeling plant architecture development based on a dynamic language. *Frontiers in Plant Science* **3**:76.

Brooking IR. 1996. Temperature response of vernalization in wheat: a developmental analysis. *Annals of Botany* **4**:507–512.

Brown HE, Jamieson PD, Brooking IR, Moot DJ, Huth NI. 2013. Integration of molecular and physiological models to explain time of anthesis in wheat. *Annals of Botany* **112**:1683–1703.

Byrne S, Guiney E, Barth S, Donnison I, Mur LAJ, Milbourne D. 2009. Identification of coincident QTL for days to heading, spike length and spikelets per spike in *Lolium perenne* L. *Euphytica* **166**:61–70.

Chapman DF, Lee JM, Waghorn GC. 2014. Interaction between plant physiology and pasture feeding value: a review. *Crop and Pasture Science* **65**:721–734.

Colvill KE, Marshall C. 1984. Tiller dynamics and assimilate partitioning in *Lolium perenne* with particular reference to flowering. *Annals of Applied Biology* **104**:543–557.

Cooper JP. 1950. Day-length and head formation in the ryegrasses*. *Grass and Forage Science* **5**:105–112.

Cooper JP. 1960. Short-day and Low-temperature Induction in *Lolium*. *Annals of Botany* **24**:232–246.

Coussement J, Henke M, Lootens P, Roldán-Ruiz I, Steppe K, De Swaeft T. 2018. Modelling leaf spectral properties in a soybean functional-structural plant model by integrating the prospect radiative transfer model. *Annals of Botany* **122**:669–676.

Davies A. 1971. Growth rates and crop morphology in vernalized and non-vernalized swards of perennial ryegrass in spring. *The Journal of Agricultural Science* **77**:273–282.

Durand J-L, Andrieu B, Barillot R, Barre P, Combes D, Enjalbert J, Escobar Gutiérrez A, Faverjon L, LeCarpentier C, Litrico I, Louarn G, Migault V, Sanchez Rodriguez L. 2016. Progresser dans la simulation mathématique des performances des mélanges de variétés fourragères pour composer et améliorer les prairies. *Fourrages* **225**:21–28.

Fiorelli CF, Woodward SJR, Wastney ME, Thom ER, Bahmani I. 2001. Modelling factors affecting reproductive development of perennial ryegrass in Waikato dairy pastures. *Proceedings of the Conference-New Zealand Grassland Association* **63**:165–170.

Gauthier M, Barillot R, Schneider A, Chambon C, Fournier C, Pradal C, Robert C, Andrieu B. 2020. A functional-structural model of grass development based on metabolic regulation and
coordination rules (R Hancock, Ed.). *Journal of Experimental Botany* **71**:eraa276.

GEVES. 2006. *Classement des variétés de plantes fourragères inscrites en liste A sur le Catalogue Officiel français*. Lusignan, France: CTPS - Comité Technique Permanent de la Sélection.

Gillet M. 1980. *Les graminées fourragères: description, fonctionnement, applications à la culture de l'herbe*. Paris, France: Gauthier-Villars.

Gonthier R, Francis D. 1989. Changes in the pattern of cell division in the shoot and root meristems of *Lolium perenne* during the transition from vegetative to floral growth. *Journal of Experimental Botany* **40**:285-292.

Hay RKM, Kemp DR. 1990. Primordium initiation at the stem apex as the primary event controlling plant development: preliminary evidence from wheat for the regulation of leaf development. *Plant, Cell & Environment* **13**:1005-1008.

Hazard L, Betin M, Molinari N. 2006. Correlated response in plant height and heading date to selection in perennial ryegrass populations. *Agronomy Journal* **98**:1384-1391.

He J, Le Gouis J, Stratonovitch P, Allard V, Gaju O, Heumez E, Orford S, Griffiths S, Snape JW, Foulkes MJ, Semenov MA, Martre P. 2012. Simulation of environmental and genotypic variations of final leaf number and anthesis date for wheat. *European Journal of Agronomy* **42**:22-33.

Heide OM. 1994. Control of flowering and reproduction in temperate grasses. *New Phytologist* **128**:347-362.

Jamieson PD, Semenov MA, Brooking IR, Francis GS. 1998. Sirius: a mechanistic model of wheat response to environmental variation. *European Journal of Agronomy* **8**:161-179.

Jégo G, Bélanger G, Tremblay GF, Jing Q, Baron VS. 2013. Calibration and performance evaluation of the STICS crop model for simulating timothy growth and nutritive value. *Field Crops Research* **151**:65-77.

Jouven M, Carrère P, Baumont R. 2006. Model predicting dynamics of biomass, structure and digestibility of herbage in managed permanent pastures. 1. Model description. *Grass and Forage Science* **61**:112-124.

Keep T, Sampoux J-P, Blanco-Pastor JL, Dehmer KJ, Hegarty MJ, LeDauphin T, Litrico I, Muylle H, Roldán-Ruiz I, Roschanski AM, Ruttink T, Surault F, Willner E, Barre P. 2020. High-throughput genome-wide genotyping to optimize the use of natural genetic resources in the grassland species perennial ryegrass (*Lolium perenne* L.). *G3 Genes* | *Genomes* | *Genetics* **10**:3347-3364.

Kemp DR, Eagles CF, Humphreys MO. 1989. Leaf growth and apex development of perennial ryegrass during winter and spring. *Annals of Botany* **63**:349-355.

Kipling RP, Virkajärvi P, Breitsameter L, Curnel Y, De Swaef T, Gustavsson A-M, Hennart S, Höglind M, Järvenranta K, Minet J, Nendel C, Persson T, Picon-Cochard C, Rolinski S, Sandas DL, Scollan ND, Sebek L, Seddaiu G, Topp CFE, Twardy S, Van Middelkoop J, Wu L, Bellocchi G. 2016. Key challenges and priorities for modelling European grasslands under climate change. *Science of the Total Environment* **56**:851-864.

Kunrath TR, de Berranger C, Charrier X, Gastal F, de Faccio Carvalho PC, Lemaire G, Emile J-C, Durand J-L. 2015. How much do sod-based rotations reduce nitrate leaching in a cereal cropping system? *Agricultural Water Management* **150**:46-56.

Matthew C, Assuero SG, Black CK, Sackville Hamilton NR. 2000. Tiller dynamics of grazed swards. In: Lemaire G, Hodgson J,

Moraes A, eds. *Grassland Ecophysiology and Grazing Ecology*. Wallingford, UK: CABI Publishing. 127-150.

Morris MD. 1991. Factorial sampling plans for preliminary computational experiments. *Technometrics* **33**:161-174.

Ong CK, Marshall C, Saoar GR. 1978. The physiology of tiller death in grasses. 2. Causes of tiller death in a grass sward. *Grass and Forage Science* **33**:205-211.

Paina C, Byrne SL, Domnisoru C, Asp T. 2014. Vernalization mediated changes in the *Lolium perenne* transcriptome. *PLoS One* **9:e107365**.

Parent B, Tardieu F. 2012. Temperature responses of developmental processes have not been affected by breeding in different ecological areas for 17 crop species. *New Phytologist* **194**:760-774.

Parsons AJ, Robson MJ. 1980. Seasonal changes in the physiology of S24 perennial ryegrass (*Lolium perenne* L.). 1. Response of leaf extension to temperature during the transition from vegetative to reproductive growth. *Annals of Botany* **46**:435-444.

Prusinkiewicz P. 2004. Modeling plant growth and development. *Current Opinion in Plant Biology* **7**:79-83.

R Core Team. 2022. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rouet S, Barillot R, Leclercq D, Bernicot MH, Combes D, Escobar-Gutiérrez A, Durand JL. 2021. Interactions between environment and genetic diversity in perennial grass phenology: a review of processes at plant scale and modeling. *Frontiers in Plant Science* **12**:672156-672156.

Sampoux J-P, Baudouin P, Bayle B, Béguéir V, Bourdon P, Chosson J-F, Deneufbourg F, Galbrun C, Ghesquière M, Noël D, Pietraszek W, Tharel B, Viguié A. 2011. Breeding perennial grasses for forage usage: an experimental assessment of trait changes in diploid perennial ryegrass (*Lolium perenne* L.) cultivars released in the last four decades. *Field Crops Research* **123**:117-129.

Thomas H. 1980. Terminology and definitions in studies of grassland plants. *Grass and Forage Science* **35**:13-23.

Verdenal A. 2009. De la simulation de la morphogénèse de l'appareil aérien du ray-grass anglais (*Lolium perenne* L.): exploration d'un schéma cybernétique inspiré du concept d'auto-organisation et applications. PhD Thesis, University of Poitiers, France.

Verdenal A, Combes D, Escobar-Gutiérrez AJ. 2008. A study of ryegrass architecture as a self-regulated system, using functional-structural plant modelling. *Functional Plant Biology* **35**:911-924.

Vidal T, Andrieu B. 2009. Contrasting phenotypes emerging from stable rules: A model based on self-regulated control loops captures the dynamics of shoot extension in contrasting maize phenotypes. *Annals of Botany* **126**:615-633.

Woodward SJR, Oijen MV, Griffiths WM, Beukes PC, Chapman DF. 2020. Identifying causes of low persistence of perennial ryegrass (*Lolium perenne*) dairy pasture using the Basic Grassland model (BASGRA). *Grass and Forage Science* **75**:45-63.

Yan W, Hunt LA. 1999. An equation for modelling the temperature response of plants using only the cardinal temperatures. *Annals of Botany* **84**:607-614.

Zaka S, Ahmed LQ, Escobar-Gutiérrez AJ, Gastal F, Julier B, Louarn G. 2017. How variable are non-linear developmental responses to temperature in two perennial forage species? *Agricultural and Forest Meteorology* **232**:433-442.