Are we focusing on the right parameters? Insights from Global Sensitivity Analysis of a Functional-Structural Plant Model

Abstract

Performing global sensitivity analysis on functional-structural plant models (FSP models) can greatly benefit both model development and analysis by identifying the relevance of parameters for specific model outputs. Setting unimportant parameters to a fixed value decreases dimensionality of the typically large model parameter space. Efforts can then be concentrated on accurately estimating the most important input parameters. In this work, we apply the Elementary Effects method for dimensional models with arbitrary input types, adapting the method to models with inherent randomness. Our FSP model simulated a maize stand for 160 days of growth, considering three outputs, namely yield, peak biomass and peak leaf area index (LAI). Of 52 input parameters, 12 were identified as important for yield and peak biomass and 14 for LAI. Over 70 % of parameters were deemed unimportant for the outputs under consideration, including most parameters relating to crop architecture. Parameters governing shade avoidance response and leaf appearance rate (phyllochron) were also unimportant; variations in these physiological and developmental parameters do lead to visible changes in plant architecture but not to significant changes in yield, biomass or LAI. Some inputs identified as unimportant due to their low sensitivity index have a relatively high standard deviation of effects, with high fluctuations around a low mean, which could indicate non-linearity or interaction effects. Consequently, parameters with low sensitivity index but high standard deviation should be investigated further. Our study demonstrates that global sensitivity analysis can reveal which parameter values have the most influence on key outputs, predicting specific parameter estimates that need to be carefully characterized.

1. Introduction

Models in the biological and environmental sciences typically have many free parameters (Razavi and Gupta 2015; Qian and Mahdi 2020; Puy et al. 2021). This is, in particular, true for functional-structural plant models (FSP models) (Evers and Bastiaans 2016; Henke et al. 2016; Zhu et al. 2017) which simulatpe growth and development of plants in 3D as a function of environmental factors such as light, temperature and nutrition. Calibration of these parameters often requires empirical data, which can be costly or simply impossible to obtain (see Qian and Mahdi (2020) and the references therein). However, often only a small subset of input parameters in a system has a significantinfluence on a specific system output (Box and Meyer 1986; Razavi et al. 2021). As such, it is beneficial for model development to identify unimportant parameters, so that they may be set to a fixed value and so that efforts can be concentrated on accurately estimating the most important input parameters. This can greatly decrease dimensionality of the model parameter space while increasing trust in the model. Sensitivity analysis (SA), the study of how variability in the model output can be attributed to the variability in the model inputs, is a common tool for identifying (un)important input parameters.

Local sensitivity analysis techniques in the form of one-at-a-time (OAT) methods involve changing one parameter at a time from a fixed base point and assessing the effect on the model output (Pianosi et al. 2016). This assessment may involve the use of (potentially discretized) derivatives or might simply involve visual inspection of the model outputs (Pianosi et al. 2016). Because of their simplicity and ease of implementation, these methods for local SA are commonly used, although local SA methods may only prove informative in very specific situations (e.g. inverse problems or approximating a model output in a small region of output space) (Saltelli and Annoni 2010). In general, local sensitivity analysis is not recommended for rigorous SA; global methods (GSA) should be used instead (Saltelli et al. 2008; Saltelli and Annoni 2010). Global sensitivity analysis methods address the limitations of local OAT techniques by considering the entire space of variability of the input parameters (Pianosi et al. 2016). A variety of approaches have been proposed for GSA (Fig. 3 in Pianosi et al. (2016) compares some of them), based on different philosophies and theories, leading to different notions of (global) sensitivity (Saltelli and Tarantola 2002; Razavi and Gupta 2015). While there are certainly examples of GSA in FSP modelling (Cournède et al. 2013; Mathieu et al. 2016; Streit et al. 2016; Sainte-Marie and Cournède 2019; Zhu et al. 2021), in other cases, SA has not yet been applied or analyses may consist solely of OAT or visual inspection of model outputs (Vos et al. 2007; Buck-Sorlin et al. 2011; Henke et al. 2016; Coussement et al. 2018; de Vries et al. 2020; Bailey and Kent 2021; Gauthier et al. 2021; Pao et al. 2021; van der Meer et al. 2021; Pao et al. 2021). This does not invalidate the modelling exercise, but does highlight that GSA is not yet typically incorporated as standard into the study of FSP models.

The objective of this study is to show the benefit of performing GSA in FSP models. In particular, we apply the elementary effects method (EE) (also known as Morris’ method Morris (1991)) to an FSP model simulating maize to identify (un)important parameters, assess the potential for model improvement and simplification and identify which biological processes might be more or less important for three model outputs, namely leaf area index (LAI), aboveground biomass and yield. The elementary effects method was chosen because it is a simple and commonly used GSA method designed to screen for (un)important parameters and is capable of identifying non-linear or interaction effects. EE has been used to analyse crop models (e.g. Richter et al. 2010; Vanuytrecht et al. 2014; Casadebaig et al. 2016; Caubel et al. 2017; Silvestro et al. 2017), but applications to FSP models are rare. We use the recently published version with scaled effects here, as described in Rutjens et al. (2023), which makes the method applicable to dimensional models with inputs of arbitrary range and inputs of integer or Boolean type and prevents erroneous rankings. Several other common approaches, including those that take the dependency structure of the inputs into account (Owen 2014) are described in Saltelli et al. (2008), Razavi and Gupta (2015), Qian and Mahdi (2020) and the references therein.

2. Elementary Effects Method

The idea behind EE is to characterize the space of model outputs by a relatively low number of strategically placed simulation points. From these simulation points, we are then able to calculate finite differences (called elementary effects) as measures of how the output changes when one input changes. Finally, by aggregating these effects for each combination of input and output in a certain way, we obtain measures of (global) sensitivity of the outputs for the inputs. In this work, we define the sensitivity of output $Yj$ to input parameter $Xi$ to be the relative contribution of the variability in the input parameter to the variability of the output. The asserted (finite) range of an input parameter is used here as input variability and variance in the output as output variability (which is a common choice (Razavi and Gupta 2015)), but other notions of input and output variability can in principle be used as well, such as mean, standard deviation or interquartile range. This section follows (Rutjens et al. 2023). We refer to that work and the references therein for further detail.

2.1. Definition of elementary effects

Let $Xi$⁠, $i=1,…,k$ be dimensional input parameters with units $[Xi]$⁠, taking values in $[mini,maxi]$ uniformly. If the parameter can only take integer values, it takes values in the set ${mini,mini+1,…,maxi}$⁠. The same holds for Boolean parameters, but then $mini=0$ and $maxi=1$⁠, where 0 encodes false and 1 stands for true. $xi$ denotes the dimensionless equivalent scaled to the unit interval, that is,

henceforth referred to as scaled dimensionless parameters. The dimensional outputs of interest are denoted by $Yj$⁠, $j=1,…,q$ with corresponding unit $[Yj]$⁠. An elementary effect for output $Yj$ and input parameter $Xi$ at base point $Xn=(X1n,…,Xkn)$ is then defined as the finite difference

Here, the superscript $n$ is simply an index to distinguish different $X$⁠, to emphasize that the EE can be calculated at numerous points in the parameter space. $Δin$ represents the step in the actual parameter space, which depends on the scaled dimensionless step $δin∈[0,1]$ in the following way:

The value of $δin$ follows from the chosen trajectory generation method and is calculated a posteriori (see Section 2.2 and Supporting Information—Supplementary Material S1). The effect $E⁢Ei⁢jn$ given by Equation (2) is dimensional with units $[E⁢Ei⁢jn]=[Yj]/[Xi]$⁠.

While here we assume uniformly distributed inputs, this can be relaxed to include arbitrary distributions. Scaling the sampled parameter values from $[0,1]$ to the actual interval should then be done using the corresponding inverse cumulative density function (CDF).

2.2. Trajectory generation

It is generally not feasible (nor desirable) to calculate every possible effect; even with discrete-valued inputs, the number of required simulations quickly becomes prohibitive. Instead, the goal is to generate a small set of $Q=r⁢(k+1)$ simulation points, typically $Q∼1000$⁠, which still provides good coverage of the parameter space. Here, $r$ is the number of effects for each input. The sensitivities are then characterized by some aggregation measure and standard deviation over $r$ effects. In this work, the $Q$ simulation points are clustered into $r$ star-shaped trajectories (Fig. 1; also called a radial design), where the base point (⁠$x1$ in Fig. 1) is used repeatedly in combination with the perturbed points (⁠$x2$–$x4$⁠) for the calculation of the elementary effects. To ensure a uniform distribution of the $r$ base points in the parameter space, a quasi-random (QR) or low-discrepancy sequence is typically used (Morokoff and Caflisch 1994; Hickernell 2000; Saltelli et al. 2010; Campolongo et al. 2011). QR sequences are designed to produce point sets that cover a space both efficiently (i.e. with a low number of points) and evenly (i.e. approximating a uniform distribution). In this work, we consider the recently presented $Rd$ sequence (Roberts 2021), which has been shown to perform well in EE (Rutjens et al. 2023). The QR sequence is used to generate a set of $r$ base points (gathered in a $r×k$ matrix $A$⁠) and a set of $r$ perturbation vectors (⁠$r×k$ matrix $B$⁠). Each pair of base point and perturbation vector then leads to a star-shaped trajectory in parameter space [see Supporting Information—Supplementary Material S1].

For Boolean or integer input parameters, an adjustment is required to avoid non-integer/Boolean parameter values in the actual parameter space for these parameters (Rutjens et al. 2023). We pin the base point coordinates corresponding to integer/Boolean inputs to a discrete value and then use a fixed step size. In other words, given a base coordinate $xi$ in $[0,1]$ generated by a QR sequence, we transform the coordinate to a discrete value $x~i$ by:

where $pi$ (the number of discrete values input $i$ can take) satisfies

for some $m∈ℕ$⁠. While $x~i$ is not necessarily integer/Boolean, the coordinate in the actual parameter space is given by:

for some $m∈ℕ$⁠, which is an integer/Boolean. Subsequently, the perturbed point is found by stepping with fixed step size

where $n∈{1,2,…,pi−1}$ is to be chosen, such that the perturbed point still lies in the parameter space. Whenever possible, $pi$ should be chosen to be even and $n=pi/2$ to ensure equal sampling probabilities (Morris 1991).

2.3. Scaled effects and sensitivity measure

Scaling effects in the input-direction and using a dimensionless sensitivity measure are essential when input parameters are dimensional and take values on non-unit intervals, in order to prevent erroneous rankings and to obtain results consistent with the notion of global sensitivity (Rutjens et al. 2023). This is not standard practice yet (see, e.g. Menberg et al. 2016; Pao et al. 2021et al.; Uys et al. 2021) and not all third-party GSA implementations include both of these additional steps (in particular the use of a dimensionless sensitivity measure) Iooss et al. (2021); Pianosi et al. (2015). Following Rutjens et al. (2023), we scale effects by the input range

The scaled effect $E⁢Ei⁢jn⋅cxi$ has units $[Yj]$⁠. Subsequently, we define the dimensionless normalized sensitivity measure

where $χi⁢j$ is the median of the absolute effects ${|E⁢Ei⁢jn|}n=1,…,r$⁠. The $Sχ⁢(i,j)$ take values in $[0,1]$ and sum to unity (over $i$ for fixed output $j$⁠). This allows for a standardized way of identifying the (un)important parameters (Wu 2020). For a given output $Yj$⁠, the $Sχ⁢(i,j)$’s are sorted in ascending order leading to a sequence

The $q$ unimportant parameters are then those for which

where $h$ is a predefined percentage (10 % in Wu (2020), although they note a low $h$-value may lead to incorrectly identifying parameters as important). Important parameters for a given output are those with a value of the sensitivity index $Sχ$ above a pre-determined threshold $S0⁢j$⁠. This quantity may be a function of the $Sχ⁢(i,j)$’s, making it a function of the threshold $h$ through Equation (10) (essentially $q=q⁢(h)$⁠); here

is used (as in Wu (2020) and Rutjens et al. (2023)), where $μ^0⁢j$ and $σ^0⁢j$ are the sample mean and standard deviation of the $Sχ⁢(i,j)$’s corresponding to unimportant variables, that is,

Finally, to gain insight into the non-linear response of outputs to inputs and detect potential interaction effects, we consider the relative standard deviation (RSD) defined by

where $σi⁢j$ is the standard deviation of the (non-absolute) effects ${E⁢Ei⁢jn}n=1,…⁢r$ and $μi⁢j⋆$ is the mean of the absolute effects ${|E⁢Ei⁢jn|}n=1,…⁢r$⁠. In other works (Ruano et al. 2012; Yang and Becerik-Gerber 2015; Menberg et al. 2016), ratios of $μ⋆$ and $σ$ are used to classify levels of non-linearity, interaction effects or general importance, but in this work, RSD is purely used as a tool to inform which parameters might deserve a more detailed follow-up (i.e. those with a small sensitivity index $Sχ$ but a high RSD).

3. Model description

The model considered here is a general modular FSP model using the modelling platform GroIMP (Hemmerling et al. 2008). It simulates aboveground plant growth and architectural development, driven by competition between plants for light and nutrients. It is an evolution of the model presented in Evers and Bastiaans (2016).

In this study, we simulate maize plants at a daily time step based on the principles of assimilate supply and demand. Plant organs serve either as a pure sink for assimilates (the generative organs and the roots) or both as a source and a sink depending on their age and surface area (the leaves and the stem segments).

3.1. Environment

To accurately represent a realistic light field, a mix of direct and diffuse light sources are placed in the scene. Seventy-two diffuse light sources are configured in 6 rings at different heights and 24 direct light sources describe an arc in the sky to emulate the sun’s path (Evers and Bastiaans 2016). The angle and orientation of the arc depend on the day of the year and latitude; the power of the light sources depends on the day of the year, the location in the sky, latitude and transmissivity. The relative contribution of the diffuse and direct light sources is constant over a simulation, with the fraction of diffuse light of the total incoming radiation set to 0.8. Further details of the diffuse/direct radiation calculation can be found in Goudriaan and Van Laar (1994), Spitters et al. (1986) and Spitters (1986). At each time step, rays of photosynthetically active radiation (PAR) cast by the light sources are absorbed, reflected or transmitted by leaves and stems according to organ-independent optical coefficients, using the stochastic ray tracer capabilities of GroIMP (Hemmerling et al. 2008).

The average daily temperature follows a sinusoidal pattern over the year, reaching its minimum value around day 20, mean at day 111, and its maximum around day 202 (Evers and Bastiaans 2016). It is characterised by a mean value (tav_a) and the amplitude of the deviation around that mean (tav_b). Other atmospheric parameters such as CO$2$ level, vapour pressure deficit and O$2$ level are kept constant throughout a simulation.

3.2. Physiology

At the organ level, the absorbed PAR is used to calculate the photosynthesis rate as a function of organ nitrogen level (Evers and Bastiaans 2016). Assimilated CO$2$ is converted into growth substrates, and maintenance costs are deducted. This leads to a daily pool of substrates available for organ growth.

Photosynthesis is described by a C4-equivalent of the Farquhar–von Caemmerer–Berry (FvCB) model as in Yin and Struik (2009). The FvCB model implemented here predicts the net photosynthesis rate as the minimum of the Rubisco-limited and the electron transport-limited rates of CO$2$ assimilation. The FvCB model is the standard in relating photosynthetic carbon assimilation to the concentration of intercellular CO$2$ and absorbed photosynthetically active radiation (Dubois et al. 2007). While parameter values for the original C3 model have been calibrated fairly precisely, there are a number of parameters in the C4-equivalent with more uncertainty (Yin and Struik 2009). We take nine of the FvCB parameters into account in the sensitivity analysis (Table 1, parameters 35–43).

Potential organ growth rate is defined as the organ demand for growth substrates (its sink strength), implemented using the first derivative of a beta growth function (Yin et al. 2003; Evers and Bastiaans 2016). To determine actual growth rate, the relative sink strength concept is used (Heuvelink 1996), in which the sink strength of an organ is expressed as a fraction of total plant sink strength. Depending on substrate availability and relative sink strength, organs grow at or below their potential rate. Any excess growth substrates stored from all organs are made available for growth in the next time step.

Finally, organ size is updated based on the substrates received by each organ, using parameters for leaf mass per unit of leaf area (LMA) for the leaves, and specific internode length (SIL) for the internodes. Specifically for internodes, additional extension due to shade avoidance (Huber et al. 2021) is implemented by making SIL dependent on the level of competition experienced by the plant (Evers and Bastiaans 2016).

3.3. Plant development

Plant development is temperature driven, with a daily thermal time increment depending on the average daily temperature (see Section 3.1) and the species base temperature. We suppose plants germinate on the same day. Subsequent leaves are initiated and appear at constant thermal time intervals (plastochron and phyllochron, respectively). The total number of phytomers produced during a plant’s lifespan is taken as constant and equal for all plants. If a leaf reaches a certain age or receives less light over a day than a given threshold, it will be shed.

3.4. Plant architecture

To mimic maize architecture, the first given number of internodes do not elongate (parameter 18: nrShortInternodes), resulting in the corresponding leaves emerging near the soil level. Leaves are represented by narrow oblong surfaces (see Evers and Bastiaans (2016) and Figs. 7 and 8). Consecutive leaves appear along the stem at a constant phyllotactic angle. The root system architecture is not explicitly modelled.

4. Simulations and Analysis

For each of the selected simulation points (as described in Section 2), we simulate a field of 100 maize plants until harvest at day 159. The field is cloned 10 times in each direction to eliminate border effects. Three outputs are considered: yield at the final simulation day (grain biomass per m$2$⁠), peak above ground biomass and peak LAI (in m$2$ of leaf area per m$2$ of ground area).

4.1. Model parameters

Fifty-two parameters are taken into account in this analysis. Their input ranges and brief descriptions are included in Table 1. Inputs are grouped into four categories, namely architectural, developmental, environmental and physiological (Table 1).

For each parameter, 40 effects (hence 40 trajectories) are considered in a radial design as described in Supporting Information—Supplementary Material S1, leading to a total of $40×(52+1)=2120$ simulations per replicate (Section 4.2). For determining what parameters are (un)important, an $h$-level (Equation (10)) of $30$ % is used.

Furthermore, a number of model parameters are fixed throughout the analysis as they have known and fixed values (e.g. latitude) or concern management practices (e.g. row distance). Their values and descriptions are listed in Table 2 in Supporting Information—Supplementary Material S2.

4.2. Inherent model randomness

EE is designed for deterministic models. Our model, however, contains inherent randomness, not caused by the input parameters, but, for example, by small random perturbations in internode orientation or leaf angle (to capture the natural variation of a species population) or by inherent randomness in the light submodel. To address this, we repeat each of the simulations three times, and average the outputs of the three replicates at each simulation point. We then calculate effects and sensitivity indices of these averaged values. Since the relative standard devation (RSD) of the output values of the different replicates at each simulation point—the standard deviation of three replicates divided by the mean—is low (Fig. 2), we argue three replicates is sufficient in this case.

4.3. Effect outliers

Because we use a quasi-random (QR) sequence Roberts (2021) to generate the simulation points, in a small number of cases the distance from base point to perturbed point is relatively very small (Fig. 3). This is not a problem if the model is deterministic, but since our model contains inherent randomness it can cause effect outliers, which could cause changes in the sensitivity indices. This can be seen as follows. We split the output $Y$ in a deterministic part (i.e. the value if the model was deterministic) and a random part, which leads to the (unscaled) effect for input $xi$ and trajectory $n=1,…,r$ given by:

A problem arises if there are a few $n$ for which $|δin|≪1$⁠. The deterministic part $E⁢Edet,in/δin$ should be the same order of magnitude for all $n$⁠, but for the $n$ with a very small (relative) step size the stochastic part $E⁢Ernd,in/δin$ may blow up and become dominant (relative to the other $n$ with larger step sizes).

In preliminary experiments (not shown here), we found up to a thousandfold increase in the elementary effect compared to the mean for that parameter. To remedy this issue, we remove effect outliers. Specifically, assuming coordinates generated by a QR sequence are identical and independent uniformly distributed (so $xin∼U⁢([0,1])$⁠, the scaled dimensionless step size $|δin|$ follows a triangle distribution with CDF $F⁢(|δin|)=2⁢|δin|−|δin|2$ for $0≤|δin|≤1$⁠.

Outliers are then those effects for which the step size is an extreme value of this distribution, here defined as those for which $F⁢(|δin|)<0.005$⁠. It follows that all effects for which $|δin|<1−398/20≈0.0025$ are classified as outliers, and are thus removed from the analysis. In the results presented below, we removed 13 effects distributed over 12 inputs out of a total of $40⋅52=2080$ elementary effects (Fig. 3).

5. Results

With an unimportance threshold of $h$=30 % (see Equation (10)) and the importance threshold given in Equation (11), 12 out of 52 parameters were classified as important (39 as unimportant) for the outputs yield and peak biomass, and 14 parameters were found to be important (37 unimportant) for peak LAI (Fig. 4 and Table 1). This leaves one parameter for each output that is neither unimportant norimportant.

5.1. Important parameters

A number of parameters in the model currently take a value from the same input range for all species and/or organs, although it is known these are actually species- or even organ-specific. The fact that some of these—for example, growth respiration (parameter 44: rg) or the conversion factor for CO$2$ to biomass (parameter 47: fCO$2$⁠)—are identified as important suggests it might be appropriate to make these parameters species- or organ-dependent, for example, give them different input ranges for different species or organs. Likewise, average temperature parameters being identified as important for yield and peak biomass indicates that i) it is important to have accurate temperature data or predictions; and ii) it might be beneficial to add a more detailed description (in time and/or space) of temperature.

This sensitivity analysis also suggests the description of some species-dependent parameters is too simplistic. As an example, leaf mass per unit area (parameter 6: LMA) is currently modelled as homogeneous in space and constant in time. In reality, this quantity is likely variable in both space and time (Zhou et al. 2020). This input is among the most important ones for all outputs, which implies it may be useful to include more detail in this aspect of the model.

Interestingly, some of the important parameters cannot generally be measured in an experimental setting, or at best can only be roughly estimated. Among these are theoretical maximum biomass of the root system, flower and leaf (parameters 19: wmaxRoot, 20: wmaxFlower and 22: wmaxLeaf).

5.2. Unimportant parameters

While two architectural parameters are identified as important for the output peak LAI (leaf curvature (parameter 12: leafCurve) and insertion angle of upper leaves (parameter 11: leafAngleUpper)), the vast majority are deemed unimportant for all outputs (Fig. 4). This suggests that, for the current maize model, parameters such as internode length per unit biomass (parameter 5: specificInternodeLength), petiole length per unit biomass (parameter 14: specificPetioleLength), angle between consecutive leaves (parameter 16: phyllotaxis) and a number of leaf shape parameters (parameters 7: lwRatio, 8: maxWidth, 9: shapeCoeff), which might be time-consuming or costly to measure, can be fixed or estimated more coarsely.

Furthermore, the parameters governing the strength of the shade avoidance syndrome (SAS; parameters 31 and 32) (in this case shade-induced internode extension) are identified as unimportant, even though the range for these parameters was quite wide. This does not mean that the SAS mechanism itself is irrelevant, as plant architecture changes significantly depending on the amount of competition Huber et al. (2021), but instead our analysis suggests that variation in the strength of the SAS response does not lead to significantly different LAI, yield or biomass.

Surprisingly, phyllochron (parameter 50) and plastochron (parameter 49) are also classified as unimportant. Simulations with phyllochron set to the minimum and maximum values of its input range show a noticeable effect on the architecture (Fig. 8), but, like SAS, our analysis suggests that this does not translate to a significant impact on field-level outputs yield, peak biomass or peak LAI (Fig. 6).

5.3. Relative standard deviation

A number of inputs have a low sensitivity index but high RSD Equation (14), such as PAR reflectance (parameter 33: reflectancePAR) and transmittance (parameter 34: transmittancePAR) by leaves and stem, and a photosynthesis parameter used in the calculation of the maximum rate of Rubisco activity-limited carboxylation (parameter 37: Vcmax25_a) (Fig. S6 in Supporting Information—Supplementary Material S4). While the mean response of an output to changes in such an input is low, there are high fluctuations around this mean, which would indicate non-linearity or interaction effects. This makes it unclear whether these parameters are actually unimportant. We found that there can be a variety of reasons for a high RSD (see Fig. S6 in Supporting Information—Supplementary Material S4). On one end of the spectrum are inputs with a single (absolute) effect that is much larger than the others, solely causing a high standard deviation of effects. This might indicate very local non-linearity or interaction effects, or suggest important parts of parameter space have not been sufficiently covered by the chosen number of trajectories. On the other end are the inputs for which effects are evenly distributed around the mean, indicating a more general non-linear trend or more global interaction effects. There does not seem to be a clear correlation between large effects and small step sizes for such effects (Fig. 3 and Supporting Information—Supplementary Fig. S6).

5.4. Illustrative (local) OAT simulations

We performed an illustrative (local) OAT study for three parameters that were identified as important (nitrogen content of fully lit leaf (parameter 27: nitro), maize base temperature (parameter 48: tb), conversion factor of CO$2$ to biomass (parameter 47: fCO$2$⁠)) and three inputs found to be unimportant for all outputs (leaf shape coefficient (parameter 9: shapeCoeff), phyllochron (parameter 50: phyllochron), PAR reflectance by stem and leaves (parameter 33: reflectancePAR)) (Fig. 6). The baseline is taken to be the mean of each input’s range, and the six inputs are then varied over their ranges uniformly. Consistent with the results from the GSA, the OAT method shows that the three important parameters have a significant effect on the outputs, while the unimportant parameters have negligible effects (Fig. 6). This is purely an illustrative example; OAT is not a replacement for rigorous GSA, and the output response could be different at another baseline in parameter space. Figures 7 and 8 show the effect on the plant architecture of varying leaf nitrogen content (parameter 27: nitro) and phyllochron (parameter 50: phyllochron), respectively.

5.5. Further observations

While maximum biomass is typically attained at the last simulation day, there is a large variation in the time when maximum LAI is achieved, as the precise pattern of leaf senescence hardly plays a role for biomass accumulation, while it does for LAI (Fig. 9). The simulation day at which maximum LAI is achieved approximately follows a normal distribution with mean 103 (days) and standard deviation 17.74 (days), noting that all the simulations where the maximum was at (or would be after) the final simulation day are grouped in the binfor day 159.

Plotting the sensitivity indices for yield against those for peak biomass (Fig. 10; $R2=0.8946$⁠) shows that for only 6 out of 52 parameters, the sensitivity indices lead to different classifications (e.g. important for yield, neither for peak biomass), and in only three cases are the classifications opposite (e.g. important for yield but unimportant for peak biomass), indicating biomass and yield are generally correlated. This agrees with previous findings from field trials (Sinclair et al. 1990).

Finally, we considered the evolution of the sensitivity indices for LAI, aboveground biomass and yield over time (Supporting Information–Supplementary Material S3), and used collinearity analysis (following De Swaef et al. (2019); Coudron et al. (2021)) to assess identifiability of the model. The collinearity index used here was proposed by Brun et al. (2001) as a measure of linear dependencies between model parameters. Our results show that the model is generally identifiable (Figs. S4 and S5), with all collinearity values for LAI and biomass (and most for yield) below the commonly used threshold of 20, and the majority of values for all three outputs below the stricter threshold of 15.

6. Discussion

Over recent decades, many modellers have developed functional-structural plant models to gain understanding of plant development. However, as these models become more detailed, they typically involve huge parameter sets, making it challenging to explore the parameter space and infer how the choice of parameters influences model predictions. Furthermore, the size of the parameter space makes it impractical to identify optimal parameter regions that would be predicted to maximise outputs like light capture, biomass production or grain production. Modelling studies are therefore typically limited to using a defined parameter set and exploring the influence of a handful of plant traits. A notable exception is found in Rangarajan et al. (2022), who performed multi-objective optimization on an FSP model to identify optimal root phenotypes.

Whenever models use parameter values estimated from the experimental literature, it is often unclear whether any errors in measurement or differences between species would influence model predictions and where more detailed measurements would be beneficial.

To address these challenges, we have explored the use of global sensitivity analysis (GSA), which identifies important (and unimportant) parameters for a given model output. Our study has demonstrated that GSA provides useful insights both for our biological understanding of plant development and for guiding further model development andparameterization.

Modelling insights   In terms of model development, GSA enables us to identify unimportant parameters, so that they may be set to a fixed value, thereby decreasing dimensionality of the typically large model parameter space. Efforts can then be concentrated on accurately estimating the most important input parameters. In addition, identifying the most important parameters can significantly ease the task of finding the optimal set of plant characteristics for a given trait.

Out of 52 input parameters, only 12 were identified as important for yield and peak biomass (14 for LAI), while over 70 % of inputs were deemed unimportant (Fig. 4). This highlights the benefit of incorporating GSA in the modelling process, as it suggests a rough estimate suffices for the majority of input parameter values. Our analysis revealed that the important parameters included leaf nitrogen content (parameter 27: nitro), base temperature (parameter 48: tb), conversion factor of CO$2$ to biomass (parameter 47: fCO$2$⁠) and leaf mass per unit area (parameter 6: LMA). Interestingly, some of the important parameters cannot generally be measured in an experimental setting, or at best can only be roughly estimated. Among these are theoretical maximum biomass of the root system, flower and leaf (parameters 19: wmaxRoot, 20: wmaxFlower and 22: wmaxLeaf).

The FvCB photosynthesis model is calibrated fairly precisely for C3 species, but as Yin and Struik (2009) note, several parameter estimates in the C4-equivalent contain more uncertainty. None of the FvCB parameters taken into account (parameters 35–43) are classified as important in our analysis, which suggests the C4 model can be reliably used with the values proposed in Yin and Struik (2009).

Some inputs identified as unimportant due to their low sensitivity index have a relatively high standard deviation of effects (Fig. 5, e.g., parameters 33: reflectancePAR; 34: transmittancePAR; 37: Vcmax25_a), with high fluctuations around a low mean, which could indicate non-linearity or interaction effects. To quantify the importance of these fluctuations there is no one rule that fits all, but instead one should investigate those parameters with low sensitivity indices but high RSD further on an ad hoc basis as described in Section Relative standard deviation. Expert knowledge of the biological processes or model equations related to such inputs can help identify likely interaction effects and quantify their importance.

Most parameters related to crop architecture were classified as unimportant. This does not mean that architecture itself is not relevant. If a model does not address aspects of architecture, questions that rely on architecture cannot be addressed. As changes in architectural parameters (e.g. leaf length/width ratio or phyllotaxis) do lead to noticeable changes in plant architecture, we hypothesize this is due to the choice of outputs (yield, aboveground biomass and LAI), which are crop-level outputs. As a result, many parameters that are routinely deemed to be important in conventional crop models are also identified as important in our analysis. Considering specific structural outputs (e.g. average leaf area, number of leaves) might result in more architectural parameters being classified as important. Additionally, architectural parameters might play a bigger role in polycultures, whereas here a monoculture was simulated.

More generally, sensitivity results may change under different scenarios (e.g. management or environmental) (Richter et al. 2010; Coudron et al. 2021). Architectural parameters might for example become more important with increased competition (higher plant densities) or with increased variability in plant traits (e.g. intercropping several maize varieties). Sensitivity indices may also change over time. In our analysis, this was the case for the output peak aboveground biomass (see Fig. S3), where parameters related to initial development (e.g. parameter 30: seed mass) were initially important, but became less relevant at later points in time, while, for example, parameter 27: leaf nitrogen content became more important over time. One should, therefore, carefully consider the applicability and generalizability of sensitivity results.

Biological insights   We expect several conclusions to extend to different crop species, crop designs or environmental conditions within this model. Parameters that are part of major conserved mechanisms, such as photosynthesis and growth, are likely to be equally insensitive independent of the species or conditions (within biologically sensible ranges). For stress events like drought or heat, GSA specific for such conditions would need to be done. Interestingly, we found that several architectural parameters were unimportant for model output, in contrast to some model-aided trait analyses for other species that do show the relevance of architectural parameters (Sarlikioti et al. 2011; Barillot et al. 2014; Zhu et al. 2015; Li et al. 2021). This discrepancy might have been caused by the difference in the output under consideration. We focused on LAI, biomass growth and yield, whereas the studies cited analysed the effect of architectural traits on light capture and photosynthesis without feedback on plant growth.

Parameters governing shade avoidance response were identified as unimportant. This is not surprising, since the relatively uniform maize canopy we simulated leads to approximately the same response in all plants. This would likely be different in more heterogeneous mixed-species stands where performance depends on plastic plant responses to local conditions (Zhu et al. 2015).

Leaf appearance rate (phyllochron) is known to have a significant impact on plant architecture, is typically deemed as an important parameter in crop models, and has received significant interest as a potential breeding target to bring forward the flowering date (Wilhelm and McMaster 1995; Birch et al. 1998; Padilla and Otegui 2005; Clerget et al. 2008; dos Santos et al. 2022). Interestingly, phyllochron was classified as unimportant in our study. Variations in phyllochron—as expected—did lead to visible changes in plant architecture (Fig. 8). However, these architectural changes did not lead to significant changes in field-level outputs yield, biomass or LAI. We hypothesize that the expected increase in total leaf area as a result of lower phyllochron is counteracted by increased competition between the higher number of leaves. This leads to a near-constant peak LAI, which subsequently does not lead to significant changes in peak biomass and yield.

To conclude, this work shows that including global sensitivity analysis in the modelling routine for FSP models can lead to a variety of insights about both the model and the biological processes it describes. GSA is applicable to any model, and should therefore become a standard consideration for any FSP modeller. However, as sensitivity results may change under different scenarios or models, the applicability and generalizability of sensitivity results should be inferred with caution.

Supporting Information

S1. Trajectory generation

S2. Input parameters with fixed values

S3. Exploratory data analysis

S4. Effects for three unimportant inputs with high RSD

They contain the following:

S1: Detailed information on the used trajectory generation method for our Elementary Effects analysis

S2: Table of parameter values used in EE analysis

S3: Additional figures showing results from EE analysis

S4: Figure of effects for three unimportant inputs with high relative standard deviation

Acknowledgments

This work was supported by the University of Nottingham Future Food Beacon of Excellence. We gratefully acknowledge the constructive feedback provided by two anonymous reviewers.

Conflict of Interest Statement

None declared.

Data Availability

The numerical results were obtained using code in the RGG language Kniemeyer (2008) (an extension of Java), compiled on the freely available modelling environment GroIMPHemmerling et al. (2008). Included in the code is a stand-alone implementation of our formulation of EE. The data are available in the form of Microsoft Excel spreadsheets. All code and data that support the findings of this study are openly available at https://github.com/pmxrr3/EE_applied_2023 and as Supporting Information—supplementary materials to this paper.

LITERATURE CITED