# Assessing inter-individual genetic variability in peach sugar metabolism through reliable parameter estimation of a kinetic model

Hussein Kanso, Valentina Baldazzi, Bénédicte Quilot, Charlotte Baey, Mohamed-Mahmoud Memah

► To cite this version:

Hussein Kanso, Valentina Baldazzi, Bénédicte Quilot, Charlotte Baey, Mohamed-Mahmoud Memah. Assessing inter-individual genetic variability in peach sugar metabolism through reliable parameter estimation of a kinetic model. European Journal of Agronomy, 2025, 168, pp.127641. 10.1016/j.eja.2025.127641. hal-05036459

HAL Id: hal-05036459

https://hal.science/hal-05036459v1

Submitted on 16 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Assessing inter-individual genetic variability in peach sugar metabolism through reliable parameter estimation of a kinetic model

Hussein Kanso<sup>a,b</sup>, Valentina Baldazzi<sup>c,d</sup>, Bénédicte Quilot<sup>a</sup>, Charlotte Baey<sup>e</sup>,
Mohamed-Mahmoud Memah<sup>b,*</sup>

## ARTICLE INFO

Keywords:

Peach fruit

Kinetic model

Inter-individual variability

Genetic algorithm

Nonlinear mixed model

MCMC-SAEM

NSGAII

## ABSTRACT

Integrating genetic information into dynamical models is crucial for understanding genotype variation and improving crop performance in the face of current agronomic and ecological constraints. Calibration over a large number of genotypes is a key step in building gene-to-phenotype models but remains challenging. We compared two strategies for calibrating an Ordinary Differential Equations (ODE) kinetic model simulating sugar accumulation during peach fruit development. In the first strategy, a Genotype-Based (GB) approach, the model was calibrated independently for each genotype, using either a Single-Objective Optimization (GBS) or a Multi-Objective Optimization (GBM) formulation. In the second strategy, a Population-Based (PB) approach, the model was calibrated for all genotypes simultaneously. The two strategies were applied to simulated data and to a real dataset from 106 peach genotypes. Results showed that the GB strategy allowed for a high goodness of fit for most genotypes, especially in the GBS formulation. However, the estimated parameters suffered from a lack of practical identifiability as independent repetitions of the estimation algorithm did not always converge to the same value for most genotypes. The PB calibration strategy overcame this issue showing a good identifiability of the population parameter values, a goodness of fit comparable to the one obtained with the GB strategy, and a good characterization of parameter variations within the progeny, which is key to assess the inter-individual genetic variability. This study demonstrates the value of the PB calibration strategy in capturing inter-individual genetic variability, a critical step toward developing robust gene-to-phenotype models.

## 1. Introduction

Predicting genotype-to-phenotype relationships under contrasting environments is a big challenge for plant biology and breeding. Process-based ecophysiological models are considered as powerful tools to deal with such a challenge (Martre et al., 2015) as they are able to decompose plant phenotypic traits into developmental, architectural and environmental factors. Within this framework, differences among genotypic responses are usually encoded into a set of parameters values, ideally considered as the fingerprint of the underlying genetic information (Martre et al., 2015; Bertin et al., 2010). Ultimately, the aim of gene-to-phenotype modelling is to make this link explicit, by expressing the value of genetic parameters as a function of specific genetic loci

(QTL) or allelic combinations of genes. The existence of pure genotypic parameters is a strong modelling hypothesis which is largely accepted within the community although it should be carefully checked. Indeed, some papers showed that some supposed genotypic parameters may be affected by the environment, implying that they are in fact emergent properties of the models (as GxE are emergent properties) rather than true genetic fingerprints. Interested readers are referred to (Bertin et al., 2010; Lamsal et al., 2017,2018) for further explanations and illustrative examples.

From a practical perspective, the construction of gene-to-phenotype models is challenging as it requires their calibration on a large number of genotypes, usually with few data available. Moreover, process-based models are generally nonlinear and involve a large number of

<footnote><sup>a</sup> GAFL, INRAE, Montfavet, France</footnote>

<footnote><sup>b</sup> PSH, INRAE, Avignon, France</footnote>

<footnote><sup>c</sup> ISA, Université Côte d'Azur, INRAE, CNRS, Sophia-Antipolis, France</footnote>

<footnote><sup>d</sup> MACBES, Université Côte d'Azur, Inria, INRAE, CNRS, Sophia-Antipolis, France</footnote>

<footnote><sup>e</sup> Laboratoire Paul Painlevé, Université Lille, CNRS, Lille, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: memah@inrae.fr (M.-M. Memah).</footnote>
parameters among which correlations may exist (Chou and Voit, 2009, Gutenkunst et al., 2007). In such cases, non-identifiability or equifinality issues may arise, when parameter values cannot be unambiguously inferred from available data and multiple parameter sets provide equally good fits of the data (Lamsal et al., 2018). Although in many biological applications, useful insights may be derived from model predictions rather than parameter values, uncertainty in parameter estimates poses a problem for genetic analysis, in which the parameters values are expected to be the genetic fingerprint of the system (Hass et al., 2017, Becker et al., 2010). This last expectation requires that the estimated values of genotypic parameters are stable over repeated calibrations and environments (Lamsal et al. 2018). Moreover, spurious correlations among parameter estimates along with sloppy directions in the parameter space, may prevent the identification of true genetic correlations due to epistastic effects or common regulatory mechanisms. In general, a reliable estimate of the parameters depends on the structure of the model, the quality of the observations and the choice of the calibration strategy. Although model reduction approaches can be effectively used to improve model identifiability (Maiwald et al., 2016, Chou and Voit, 2009, Snowden et al., 2017), we focus here on the last point, namely the formalization and implementation of the calibration process.

In this work, we compared two different strategies for the calibration of an ODE kinetic model simulating the accumulation of different sugars during peach fruit development (Desnoues et al., 2018, Kanso et al., 2020). The first strategy corresponded to the common approach in plant modeling community: a set of genotypic parameter is defined a priori, based on preliminary studies, sensitivity analysis or simply by biological understanding, and their estimation is addressed for each genotype independently (Genotype-based, GB). Examples of this approach can be found in the study of many plant processes, including leaf elongation (Reymond et al., 2003), plant development (Yin et al., 2000, 2005, Messina et al., 2006), phenology (Nakagawa et al., 2005, Ravi Kumar et al., 2009), nitrogen adaptation (Laperche et al., 2006), fruit growth (Constantinescu et al., 2016) and quality (Quilot et al., 2005, Prudent et al., 2011). This strategy has been largely employed as it has the advantage to be straightforward: it can be performed with any kind of optimization algorithm, applied to any model and can be easily extended to newly acquired genotypes without recalibrating all the dataset. However, as estimations are performed for each genotype independently without sharing any information, the estimated parameters using GB approach could show large variations for the whole population and even for the same genotype. In the second strategy, we calibrated the model using a Population-based (PB) approach based on mixed-effects models, which assumes that individual parameters, defined for each genotype, are independent and identically distributed random variables following a common probability distribution (Davidian and Giltinan, 2003, Baey et al., 2016). Accordingly, the genetic population is treated as a whole and the model is calibrated for all genotypes simultaneously. The target of the estimation shifts from the individual profile to the average behaviour, in terms of parameter distribution across the population. Individual parameter values can be estimated a posteriori once the population model has been fitted. PB approach has been largely used in pharmacokinetics to characterize the inter-individual variability of observed drug-response dynamics (Donnet and Samson, 2013, Owen and Fiedler-Kelly, 2014). It has been recently applied to the calibration and improvement of functional-structural plant growth models, in relation to plant genetic diversity (Perez et al., 2016, Baey et al., 2016, 2018).

Although both strategies have proven successful in a number of cases, to our knowledge, they have never been compared. The aim of this paper is to test their capacities to assess the inter-individual variability of parameter values in the context of genetic studies, using the model of sugar metabolism in peach fruit as a case study. To this aim, the two strategies were first applied on a set of simulated data to evaluate the impact of each calibration strategy on the estimation of individual genotypic parameters, in terms of accuracy and uniqueness of the

solution. Then, we applied each approach on a real dataset derived from an interspecific population of 106 peach genotypes. For the GB strategy, two formulations of the calibration problem were evaluated, either as a Single-Objective Optimization (GBS) problem or as a Multi-Objective Optimization (GBM). In the first case i.e. GBS formulation, the estimation process results in a unique value for the objective function. Hopefully, the value of this objective function corresponds to a unique set of parameter values for each genotype of the progeny. On the other hand, using the GBM formulation, an entire set of trade-offs between contrasting objective functions (Pareto front) is identified for each genotype, among which the desired parameter set has to be selected based on some other criteria. This choice could be done using graphical tools, decision-making methods, or some particular heuristics designed for specific classification problems (Coello et al., 2007, Akdemir et al., 2019). This allows for a deeper analysis of possible trade-offs among physiological traits, in the case of complex phenotypes.

The paper is organized as follows. Section 2 details the materials and methods including mathematical model, the problem formulation and the calibration methods, the simulated and experimental datasets, and the indicators used for assessing the quality of the two calibration strategies. Results are reported in Section 3 from the point of view of gene-to-phenotype modelling. The different approaches and the resulting estimation are first compared based on the quality of their phenotypic predictions i.e. the agreement to the data. The accuracy and robustness of parameter estimates is then investigated, using both a simulated and a real dataset. Finally, a discussion concludes the paper in Section 4.

## 2. Materials and methods

### 2.1. Mathematical model

The model used in this paper is the one proposed by Kanso et al. (2020) as a reduced version of the model previously developed by Desnoues et al. (2018). The model describes the accumulation of four different sugars (sucrose, glucose, fructose, and sorbitol) in peach fruit during its development, as a system of ordinary differential equations (see section S7). Carbon enters the fruit from the plant sap which is transformed by a metabolic network, including enzymatic reactions and transport mechanisms between the cytosol and the vacuole (see Fig. 1). A total of 17 reactions are included in the model and represented by linear flows. The model makes explicit use of experimental data to describe the temporal evolution of fruit carbon supply and fruit volume. To this aim, measured fruit dry and fresh masses are represented by genotype-specific temporal functions and provided as input to the model (see section S1, Figure S1, and Table S1 of the supplementary information). The model contains 30 parameters, among which 21 are considered as known (measured or taken from literature) and identical for all genotypes. This work focuses on the estimation of the remaining 9 parameters (Table 1), supposed genotype-dependent.

### 2.2. Problem formulations and calibration strategies

The aim of this work was to estimate the values of the 9 genotype-dependent parameters using data on sugar concentrations, at six developmental stages, from a large number of genotypes. Four main sugars are considered: sucrose, glucose, fructose, and sorbitol. Two calibration strategies were compared i) parameters were estimated independently for each genotype (GB strategy) and ii) parameters were simultaneously estimated for the whole progeny (PB strategy). Fig. 2 gives a schematic representation of the calibration steps, detailed below.

We consider a common framework for the two calibration strategies. Let us denote by $y_{ijk}, i = 1, ..., N_{jk}, j = 1, ..., 4, k = 1, ..., G$ the i-th observation of the j-th sugar concentration for genotype k, and by $M_{ij}$ the corresponding model prediction. $N_{jk}$ represents the total number of measurements available for genotype k and sugar j. It may include
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Schematic network of the reduced model of sugar accumulation in the peach fruit. Arrows represent carbon flows that are annotated by the related parameters being estimated (see Table 1). The corresponding kinetic equations are presented in section S7 of the supplementary information.</figcaption></figure>

<table><caption>**Table 1**

Description of the 9 unknown genotype-dependent parameters of the sugar model.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Range</th><th>Unit</th></tr></thead><tbody><tr><td>LHx</td><td>sucrose proportion hydrolized in the apoplasm</td><td>[0, 1]</td><td>-</td></tr><tr><td>RAI</td><td>coefficient of the transfer function between sucrose and (glucose +fructose) under action of acid invertase (AI) enzyme</td><td>[0, 0.1]</td><td>day<sup>-1</sup></td></tr><tr><td>KFk</td><td>fructokinase (FK) affinity</td><td>[0, 100]</td><td>mggFW<sup>-1</sup></td></tr><tr><td>KHk</td><td>hexokinase (HK) affinity</td><td>[0, 500]</td><td>mggFW<sup>-1</sup></td></tr><tr><td>OCp</td><td>coefficient of the transfer function between hexoses phosphate and other compounds</td><td>[250,1500]</td><td>day<sup>-1</sup></td></tr><tr><td>RSS</td><td>coefficient of the transfer function between hexoses phosphate and sucrose</td><td>[0, 1000]</td><td>day<sup>-1</sup></td></tr><tr><td>TAS</td><td>coefficient of sucrose transport (active import) from cytosol to vacuole</td><td>[0, 400]</td><td>mggFW<sup>-1</sup>day<sup>-1</sup></td></tr><tr><td>TPF</td><td>coefficient of fructose passive transport between cytosol and vacuole and in the opposite direction</td><td>[0, 150]</td><td>mggFW<sup>-1</sup>day<sup>-1</sup></td></tr><tr><td>RSO</td><td>coefficient of the transfer function between sorbitol and glucose under action of sorbitol oxydase (SO) enzyme</td><td>[0, 10]</td><td>day<sup>-1</sup></td></tr></tbody></table>

measurements at different time-points as well as biological replicates, if available. The index *i* stands for time point and replicates.

We then denote by $\mathbf{y}_k := (y_{11k}, \dots, y_{N_{jk}1k}, \dots, y_{14k}, \dots, y_{N_{jk}4k})^T$ the vector of all observations for genotype *k*, and by $\mathbf{M}_k$ the associated vector of all model predictions $M_{ij}$ for genotype *k*. We assume the following model:

$$ \mathbf{y}_k = \mathbf{M}_k(\boldsymbol{\phi}_k) + \epsilon_k, \epsilon_k \sim N(\mathbf{0}, \boldsymbol{\Sigma}_k) \quad (1) $$

where $\boldsymbol{\phi}_k$ is the vector of genotypic parameters to be estimated, $\epsilon_k$ is the residual error assumed to follow a Gaussian distribution with mean $\mathbf{0}$ and covariance matrix $\boldsymbol{\Sigma}_k$ of the observations $y_k$. Three parsimonious models were driven from the above one based on some hypotheses detailed in the following subsections.

### 2.2.1. Model calibration using genotype-based strategy

In the framework of this first strategy, we used two alternative formulations of the model calibration as a single-objective optimization (GBS), and as a multi-objective optimization (GBM) problem.

#### 2.2.1.1. Calibration as Single-Objective Optimisation (GBS).
The single-objective optimisation formulation combines observations of the four sugars together into a single objective defined as the sum of the squared errors of each sugar.

Using Eq. (1), we assume that the variance of the error term is independent of the sugar type i.e. $\boldsymbol{\Sigma}^{(k)} = (\sigma_{GBS_k})^2 \mathbf{I}_{N_k}$, where $N_k = \sum_{j=1}^4 N_{jk}$ is the total number of measurements for genotype *k* and $\mathbf{I}_p$ is the identity

matrix of size $p \times p$, i.e. We denote by $\boldsymbol{\theta}_{GBS,k} = (\boldsymbol{\phi}_{GBS_k}, \sigma_{GBS_k})^T$ the vector of parameters for genotype *k*. The log-likelihood function for genotype *k* is denoted by $l_{GBS,k}$ and assuming independence between data within a given genotype, is defined as follows:

$$ l_{GBS,k}(\boldsymbol{\theta}_{GBS,k}) = \frac{-N_k}{2}\log(2\pi) - \frac{N_k}{2}\log((\sigma_{GBS,k})^2) - \frac{1}{2(\sigma_{GBS,k})^2}\sum_{j=1}^4 \sum_{i=1}^{N_{jk}} (y_{ijk} - M_{ij}(\boldsymbol{\phi}_{GBS,k}))^2 \quad (2) $$

A maximum log-likelihood estimator of $\boldsymbol{\theta}_{GBS,k}$ is defined as a solution to the maximization problem:

$$ \hat{\boldsymbol{\theta}}_{GBS,k} = \underset{\boldsymbol{\phi}_{GBS,k}}{\operatorname{argmax}} (l_{GBS,k}(\boldsymbol{\theta}_{GBS,k})) \quad (3) $$

In this Gaussian case, the maximum log-likelihood estimator is equivalent to the ordinary least-square estimator:

$$ \hat{\boldsymbol{\phi}}_{GBS,k} = \underset{\boldsymbol{\phi}_{GBS,k}}{\operatorname{argmax}} \left\{ - \sum_{j=1}^{4} \sum_{i=1}^{N_{jk}} (y_{ijk} - M_{ij}(\boldsymbol{\phi}_{GBS,k}))^2 \right\} \quad (4) $$

$$ (\hat{\sigma}_{GBS,k})^2 = \frac{1}{N_k} \sum_{j=1}^{4} \sum_{i=1}^{N_{jk}} (y_{ijk} - M_{ij}(\hat{\boldsymbol{\phi}}_{GBS,k}))^2 \quad (5) $$

#### 2.2.1.2. Calibration as Multi-Objective Optimisation (GBM).
In the multi-objective optimisation formulation, the optimisation is performed on a vector of objective functions. Contrary to GBS, each sugar is considered independently. Thus, the problem can be decomposed into four
<figure><img src="image_4.png" /><figcaption>Fig. 2. Graphical representation of the proposed steps to calibrate the sugar model. Two strategies were used to calibrate the model: 1) Genotype-Based strategy where two formulations were considered: Single-Objective Optimisation (GBS) and Multi-Objective Optimisation (GBM), 2) Population-Based (PB) strategy. To estimate the unknown 9 genotype-dependent parameters (See Tables 1), 3 algorithms were used: GA for GBS, NSGA-II for GBM and MCMC-SAEM for PB.</figcaption></figure>

objective functions, each of which addressing one sugar. These functions are considered as separate objectives to be optimised independently. The idea is to estimate the parameters by finding a solution that balances between objective functions. In this context, and using Eq. (1), we have $\Sigma_k = diag(\Sigma_{1k}, ..., \Sigma_{4k})$, where $\Sigma_{jk} = (\sigma_{jk})^2 I_{N_{jk}}$, i.e. we assume that the variance of the error term depends on the sugar type.

We denote by $\theta_{GBM,jk} = (\phi_{GBM,jk}, \sigma_{GBM,jk})^T$ the vector of parameters for genotype k and sugar j and by $\theta_{GBM,k} = (\phi_{GBM,k}, (\sigma_{GBM,jk})_{j=1,...,4})^T$ the vector of parameters for genotype k. The log-likelihood function for genotype k and sugar j is denoted by $l_{GBM,jk}$, and under the assumption of independence between the observations, it is given by:

$$l_{GBM,jk}(\theta_{GBM,jk}) = \frac{-N_{jk}}{2}\log(2\pi) - \frac{N_{jk}}{2}\log((\sigma_{GBM,jk})^2) - \frac{1}{2(\sigma_{GBM,jk})^2}\sum_{i=1}^{N_{jk}}(y_{ijk} - M_{ij}(\phi_{GBM,k}))^2 \quad (6)$$

We denote by $l_{GBM,k}(\theta_{GBM,k}) = (l_{GBM,1,k}(\theta_{GBM,1,k}), ..., l_{GBM,4,k}(\theta_{GBM,4,k}))^T$ the vector of log-likelihoods for genotype k, and the

maximum likelihood estimator for genotype k is given by:

$$\hat{\theta}_{GBM,k} = \underset{\phi_{GBM,k}}{\operatorname{argmax}}(l_{GBM,k}(\theta_{GBM,k})) \quad (7)$$

Maximising each log-likelihood $l_{GBM,jk}$ lead to the same explicit formulation for $\phi_{GBM,jk}$ and $\sigma_{GBM,jk}$ as in Eqs. (4) and (5). Note that here $\sigma_{GBM,jk}$ depends on j as we assumed that the error term depends on the sugar type.

### 2.2.1.3. Model calibration using population-based strategy.
We used the statistical approach developed by Baey et al. (2018) for a simultaneous parameter estimation of all the genotypes in the population. In the PB strategy, the statistical model is defined as a two-stage hierarchical approach. The first stage deals with the intra-individual variability. The observations of each genotype are modeled as in Eq. (1), with $\Sigma_k = \sigma_{PB}^2 I_{N_k}$.

The second stage deals with the inter-individual variability. More precisely, the individual parameters $\phi_k$ are considered as random variables with population mean $\beta$ and covariance matrix $\Gamma$. In our case, log transformations were applied to these parameters as they are all strictly
positive by definition:

$$ \varphi_k := \log(\varphi_k) = \beta + \xi_k, \xi_k \sim N(0, \Gamma) \quad (8) $$

In the sequel, the term random effects will refer equivalently to $\varphi_k$ or $\varphi_k$.

We denote by $\theta_{PB} = (\beta, \Gamma, \sigma_{PB})^T$ the vector of parameters for the PB model. We assume that the random effects $\varphi_k$ are mutually independent, and that the within-genotype observations are conditionally independent, given the random effects. Thus, the log-likelihood of the PB model is defined as:

$$ l_{PB}(\theta_{PB}) = \log \left( \int \prod_{k=1}^{G} \prod_{j=1}^{4} \prod_{i=1}^{N_{jk}} f(y_{ijk} | \varphi_k; \theta_{PB}) \prod_{k=1}^{G} f(\varphi_k; \theta_{PB}) d\varphi \right) \quad (9) $$

where $f(y_{ijk}|\varphi_k; \theta_{PB})$ is the conditional density of the observations $y_{ijk}$ given the individual parameters, and $f(\varphi_k; \theta_{PB})$ is the marginal density of $\varphi_k$.

Due to the non-linearity of the mathematical model M in Eq. (1), an analytical computation of the integral (9) is impossible in practice. However, mixed models can be seen as a special case of incomplete data models, with $(y_{ijk})$ being the observed data, and the random effects $\varphi_k$ being the missing data. In this case, maximum likelihood estimation can be done using an EM-type algorithm (Dempster et al., 1977).

Confidence intervals can be obtained for each component of $\theta_{PB}$ using a parametric Bootstrap approach once the EM-type algorithm has converged (Efron and Tibshirani, 1993).

## 2.3. Optimisation algorithms and their settings

For the GB strategy, we used two genetic algorithms for the numerical optimisation of the objective functions: GA (Goldberg, 1989, Conn et al., 1991, 1997) for the GBS formulation and NSGA-II (Deb et al., 2002) for GBM formulation.

For the PB strategy, we used a MCMC-SAEM algorithm (Kuhn and Lavielle, 2004, 2005), coupling a Markov chain Monte Carlo algorithm (Hastings, 1970) with a Stochastic Approximation Expectation-Maximization (Robbins and Monro, 1951). This method allows for an estimation of $\theta_{PB}$, but individual estimates for $\varphi_k$ (or equivalently $\varphi_k$) can also be obtained as a by-product of the algorithm. More details on the implementation of this algorithm in the context of plant growth modelling can be found in Baey et al. (2018).

Given the stochastic nature of the algorithms, the estimation procedure was repeated multiple times to ensure the good exploration of the parameter space and test the robustness of the approach. For GA and NSGA-II, the estimation procedure was repeated 10 (real data) to 20 (simulated data) times for each genotype k, whereas for MCMC-SAEM, the estimation procedure was repeated 5 (simulated data) to 10 (real data) times, for all genotypes simultaneously. In the case of GBS, the results of all repetitions were combined together. Then, the solution having the maximum likelihood was selected as the reference estimation, for each genotype k. In addition, the solutions having likelihoods highest than 95th percentile of all likelihood values were kept to analyse the intra-individual variability of the estimated parameter values. Unlike in GBS, GBM does not yield a single optimal solution but identifies a Pareto-optimal set i.e. solutions corresponding to the best reached trade-offs between our four likelihoods i.e. Pareto front. The issue is then how to select a single reference solution among all these options (Constantinescu et al., 2016). Many methods can be used for this purpose including graphical tools and multicriteria decision making methods. Here, all solutions resulting from the independent repetitions of the calibration process were first pooled together and then filtered in order to identify the final Pareto-optimal set. We then summed the four individual likelihoods for each solution in this Pareto-set and choose solution(s) with highest value of this sum. In the case of PB, the

estimation having the highest likelihood value was selected as the reference solution and used for the subsequent analyses.

## 2.4. Simulated and experimental data

### 2.4.1. Simulation study

Simulated datasets were generated in order to assess the performance of the different calibration approaches.

Given the intrinsic difference between the GB and PB strategies, distinct procedures were used and are detailed in the following sections.

#### 2.4.1.1. Datasets for the genotype-based strategy.
In the case of the GB strategy (GBS and GBM formulations), we used the model to build S (S=100) simulated datasets corresponding to repeated noisy measurements of the same genotype. Accordingly, a single set of parameters (see Table 2) and the same model input (see section S1, Figure S1, and Table S1 of the supplementary information) were used to simulate the dynamics of the four sugars by means of our ODE model. In order to be close to the experimental limitations, only 6 dates were selected along the dynamics and used as observations for the subsequent calibration. In order to obtain S different measurements, a random observation noise with variance $\sigma^2 = 0.05$ was added to the simulated values, at the considered time points.

#### 2.4.1.2. Datasets for the population-based strategy.
For the population-based approach, S populations (S=100) were generated, each one including G = 100 different genotypes. All the nine effects were considered as random and $\Gamma$ was assumed to be diagonal. The population mean $\beta$ was defined as the logarithm of the values in Table 2 and the variances were equal to 10 % of the mean, i.e. diag($\Gamma$) = 0.10 $\beta$.

For each simulated population, a set of G individual parameters was generated according to Eq. 8. Observations were then obtained by running the model according to Eq. 1, with $\sigma^2 = 0.05$. As in the previous case, six dates were selected, evenly distributed along the fruit developmental period in order to be close to experimental conditions. For each genotype, model inputs (fruit dry and fresh weights) were randomly assigned (see section S1, Figure S1, and Table S1 of the supplementary information).

### 2.4.2. Experimental data

A progeny of 106 genotypes was used in this study. As described by Quilot et al. (2004), these genotypes come from an inter-specific progeny obtained by two subsequent back-crosses between Prunus davidiana (Carr.) P1908 and Prunus persica (L.) Batsch 'Summergrand' and then 'Zephyr'. All trees are planted in the orchard of the National Institute for Agriculture, Food and Environment (INRAE) in Avignon (southern France). In 2012, the concentrations of different metabolites, namely sucrose, glucose, fructose, sorbitol, and hexoses phosphates, the fruit flesh fresh mass and the flesh dry matter content, as well as enzymatic capacities (maximal activity) of twelve enzymes were measured at six time points during fruit development, for all genotypes, as described in

<table><caption>**Table 2**

Parameters used to simulate the dynamics of sugar concentrations, in the GB strategy.</caption><thead><tr><th>LHx</th><th>0.03</th></tr></thead><tbody><tr><td>RAI</td><td>0.01</td></tr><tr><td>KFk</td><td>12.18</td></tr><tr><td>KHK</td><td>20.08</td></tr><tr><td>OCp</td><td>665.14</td></tr><tr><td>RSS</td><td>33.11</td></tr><tr><td>TAS</td><td>20.08</td></tr><tr><td>TPF</td><td>12.18</td></tr><tr><td>RSO</td><td>1.01</td></tr><tr><td>σ<sup>2</sup><sub>k</sub></td><td>0.05</td></tr></tbody></table>
Desnoues et al. (2014). At each time point, available data consists in 3 biological measurements for 10 genotypes and 1 biological measurement for the remaining 96 genotypes. Genotypes in the studied progeny can be sorted into two classes: those having a 'standard phenotype', namely a balanced fructose-to-glucose ratio at maturity, and those with a 'low fructose phenotype' i.e. a lower proportion of fructose than glucose at maturity (Desnoues et al., 2018).

## 2.5. Quality indicators

The different estimation methods were evaluated using several criteria dealing with i) the quality of phenotypic predictions, and ii) the accuracy and variability of the parameter estimates.

### 2.5.1. Quality of phenotypic predictions

#### 2.5.1.1. Correlation between observations and predictions.
Pearson's linear correlation coefficient Corr<sub>j</sub> was used to quantify the relationship between predictions and observations over the whole peach population, for each sugar j:

$$ \text{Corr}_j = \frac{\sum_{k=1}^{G} \sum_{i=1}^{N_{jk}} (y_{ijk} - \bar{y}_{jk})(M_{ij}(\phi_k) - \bar{M}_j(\phi_k))}{\sqrt{\sum_{k=1}^{G} \sum_{i=1}^{N_{jk}} (y_{ijk} - \bar{y}_{jk})^2 \sum_{k=1}^{G} \sum_{i=1}^{N_{jk}} (M_{ij}(\phi_k) - \bar{M}_j(\phi_k))^2} \quad (10) $$

Where $\bar{y}_{jk} = \frac{1}{N_{jk}} \sum_{i=1}^{N_{jk}} y_{ijk}$ and $\bar{M}_j(\phi_k) = \frac{1}{N_{jk}} \sum_{i=1}^{N_{jk}} M_{ij}(\phi_k)$. Values of the correlation coefficient can range from -1 to +1. A value of -1 indicates perfect negative correlation, while a value of +1 indicates perfect positive correlation. A value of 0 indicates no correlation between the variables. Although Pearson's linear correlation coefficient is widely used to assess linear relationship between two variables, it shows some limitations (assumes linearity and normality, sensitive to outliers and sample size, etc.). For example, regression dilution can occur when errors in the observed values are not ac-counted for, which can lead to an underestimation of the correlation coefficient between observed and predicted (Piñeiro et al., 2008; Correndo et al., 2021). This can make the model appear less accurate than it actually is. One possibility to overcome this limit is to combine this criterion with other goodness-of-fit metrics to complement our fitting evaluation.

#### 2.5.1.2. Normalized Root Mean Square Error of model predictions.
The Normalized Root Mean Square Error (NRMSE) provides a measure of the agreement between model predictions and data all along the dynamics of sugar accumulation. The NRMSE for genotype k and sugar j was defined as:

$$ NRMSE_{jk} = \frac{\sqrt{\frac{1}{N_{jk}} \sum_{i=1}^{N_{jk}} (y_{ijk} - M_{ij}(\phi_k))^2}}{\max(y_{ijk}) - \min(y_{ijk})} \quad (11) $$

In order to get a comprehensive view of the quality of model predictions, a weighted NRMSE (aNRMSE) was constructed to evaluate the corresponding sugar composition of the fruit as

$$ aNRMSE_k = \sum_{j=1}^{4} \alpha_j NRMSE_{jk} \quad (12) $$

where NRMSE<sub>jk</sub> is the Normalized Root Mean Square Error computed for sugar j and α<sub>i</sub> the relative proportion of sugar i in the fruit. α = (0.72, 0.13, 0.09, 0.05)<sup>T</sup> for sucrose, glucose, fructose, and sorbitol, respectively.

### 2.5.2. Quality of parameter estimates

The definition and assessment of the quality of parameter estimates vary depending on the type of dataset under consideration. We propose here three different indicators, each with a specific scope and

application domain. In the case of simulated datasets, the real values of parameters are known, allowing the definition of distance measures between the estimated and the real parameter values. This is the case of the Root Mean Squared error and the Expected Error. In contrast, when dealing with experimental datasets, the real parameter value is unknown. In this case, the quality of parameter estimates is related to their reproducibility as measured by the Relative Standard Deviation.

#### 2.5.2.1. Root mean squared error.
The root mean squared error (RMSE) of an estimator $\hat{\theta}$ is a measure of its accuracy. It quantifies the risk of using this estimator as a proxy for the real parameter value $\theta$, in the framework of parameter estimation. The lower the RMSE value, the more the estimator is reliable.

The RMSE is defined as:

$$ \text{RMSE}(\hat{\theta}) = \sqrt{\text{Bias}^2(\hat{\theta}) + \text{Var}(\hat{\theta})} \quad (13) $$

In practice, the bias is estimated as the difference between the average value of our reference estimates (as a proxy of the expected value) and the correct value which we are trying to predict:

$$ \text{Bias}(\hat{\theta}) = \frac{1}{S} \sum_{s=1}^{S} (\hat{\theta}_s - \theta) \quad (14) $$

where $\phi$ is the value used to generate the simulated data. The variance is estimated by

$$ \text{Var}(\hat{\theta}) = \frac{1}{S} \sum_{s=1}^{S} (\hat{\theta}_s - \bar{\theta})^2 \quad (15) $$

where $\bar{\theta}_S = \frac{1}{S} \sum_{s=1}^{S} \hat{\theta}_s$ is the mean value of our reference estimates.

GB and PB approaches differ in the set of parameters $\theta$ that they aim to estimate. Accordingly, the RMSE can be computed for the nine model parameters $\phi_k$ and the observation noises in the GB strategy, whereas for the PB strategy, the RMSE applies to the parameters associated with the population parameters but also the G individual parameters).

#### 2.5.2.2. Expected error on individual parameters.
In addition to the estimation of population parameters $\theta$, PB approach can also be used to derive estimates of the individual parameters $\phi$.

In the context of our simulation study, the quality of individual parameters estimation (by GB or PB approaches) can be assessed by computing the relative distance between the true value $\phi_{r,k}$ and the estimated value for the r-th parameter and the k-th genotype. Accordingly, the Expected Error $E_{r,k}$ (in %) was defined as:

$$ E_{r,k} = \left| \frac{\hat{\phi}_{r,k} - \hat{\phi}_{rk}}{\phi_{rk}} \right| \times 100 \quad (16) $$

Overall, for each of the 9 parameters, our simulation study yielded S values of $E_{r,k}$ in the case of GB approach and S × G values for the PB approach.

## 2.6. Relative standard deviation

The relative standard deviation (RSD) was used to quantify the estimate variability of the r-th parameter of genotype k, over R repetitions of the calibration algorithm. RSD is defined as:

$$ RSD_{rk} = \frac{\sqrt{\frac{1}{R-1} \sum_{h=1}^{R} (\hat{\phi}_{r,k}^{(h)} - \hat{\phi}_{rk})^2}}{\hat{\phi}_{rk}} \times 100 \quad (17) $$

where $\hat{\phi}_{rk} = \frac{1}{R} \sum_{h=1}^{R} \hat{\phi}_{r,k}^{(h)}$ is the average estimate of the considered parameter. Contrary to the RMSE and the Expected error, RSD can be computed both on real and simulated datasets.
# 3. Results

## 3.1. The Genotype- and Population-based strategies provide comparable fit quality on real data

GB and PB strategies were applied to a population of 106 peach genotypes for which data on fruit sugar composition were available (Section 2.4.2). As a widely used approach, we first plotted observed versus predicted values to visually assess the quality of phenotypic prediction. Fig. 3 shows the predictions - observations plots obtained for the four sugars using GB and PB calibration approaches. In general, both estimation strategies provided acceptable predictions for the four sugars, with particular emphasis on sucrose, the dominant sugar in peach fruit. Predictions were less accurate for sorbitol compared to the other sugars, regardless of the calibration method. This can be partially explained by the quality of experimental measurements, which can suffer from a significant error due to the low concentration of sorbitol in the fruit.

Scatter plots are useful tool for evaluating model performance. However, this evaluation is mostly visual since some issues can arise when computing OP (observed vs predicted) or PO (predicted vs

observed) regression (Piñeiro et al., 2008; Correndo et al., 2021). One possibility to overcome this limit is to combine different goodness-of-fit metrics to perform our fitting evaluation.

The visual results were confirmed by computing correlation values between predictions and observations for each problem (see Fig. 3). GBS and PB methods gave similar if not identical results for the prediction of sugar concentrations. GBM performed well on sucrose but was clearly outdistanced on minor sugars. Despite the multi-objective definition of the corresponding likelihood function, the fitting quality of GBM appeared to be driven by the dominant sugar. Moreover, differences in the fitting quality could be observed among individual genotypes.

When looking at the prediction of the sugar accumulation dynamics, most aggregate NRSME values ranged between 0.15 and 0.25, for both GB and PB methods, but outliers existed with large prediction errors, above 0.4 (see Fig. 4 and Figure S2 for more details). In particular, a handful of data points at low glucose concentration were consistently overestimated (Fig. 3), independently of the calibration algorithm.

These points suggest the presence of atypical genotypes in the population, for which the underlying model is not well-suited.

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Observations vs. predictions using the individual estimates over the peach population. Top: results from (GBS/GA). Middle: results from (GBM/NSGA-II). Bottom: results from (PB/MCMC-SAEM). Corr: Correlations computed between predictions and observations for each sugar and each calibration formulation/algorithm. The solid line corresponds to the bisector of the first and third quadrants $y = x$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Aggregate Normalized Root Mean Square Errors (aNRMSE) between model predictions and experimental data using GB and PB calibration approaches. The boxplots show the variability of aNRMSE over the peach population.</figcaption></figure>

## 3.2. Population-based strategy provides higher statistical accuracy

In the previous section, we were concerned by the overall goodness-of-fit of each calibration strategy. We now investigate the link between the phenotype and the values of the genetic parameters. Indeed, in order to be used for genetic studies, parameter estimations should be robust across repetitions, providing a unique mapping between the individual phenotype (sugar concentrations) and the underlying parameter values. In this section, we used data from our simulation study to evaluate the performance of each algorithm to estimate the model parameters, based on root mean squared error (RMSE, Section 2.5.2.1) of the associated statistical estimator. Note that given the design of our simulation study (high number of repetitions, large population size), the measured performances in terms of RMSE have to be considered as an upper limit of the statistical accuracy of the PB and GB methods. Results are reported in Table 3.

Both GB strategies resulted in larger RMSE values compared to PB. A large variability emerged among the RMSE values for individual

parameters, with values ranging from 0.003 to 1.2. When compared to the estimations obtained by means of the single objective formulation (GBS/GA), the GBM/NSGA-II was generally less accurate, yielding larger or equivalent RMSE values. In particular, smaller RMSEs were obtained for parameters LHx, RAI, KFk, when GBS/GA was used. One parameter, TPF, however, stood out as poorly predicted by the GBS/GA approach compared to the GBM/NSGA-II approach. It is interesting to notice that TPF corresponded to a very low-sensitive parameter of the GBS/GA objective function and this could explain why the algorithm wandered around a wide range of values without any impact on the prediction quality (see Figure S5 of the supplementary information). On the other hand, this parameter ranked third for the objective function of fructose and is involved in the interactions for the remaining sugars (see Figure S4 of the supplementary information). Therefore, the GBM/ NSGA-II better accounted for these multiple involvements and thus this could partly explain the difference between the GBS and GBM results in favour of the last one.

PB approach provided a reliable estimation of distribution parameters $\theta_{PB}$ as proved by the low RMSE values ranging from $1e^{-4}$ to $6e^{-2}$ (see Table 3). Good performances for the estimation of $\theta_{PB}$ also extended to the estimation of individual genotypic parameters. Fig. 5 shows a comparison between estimated and true individual parameter values obtained for the three considered algorithms (GBS/GA, GBM/NSGA-II, PB/SAEM), out of the 100 simulated datasets. Results clearly point to a greater accurateness of the PB approach, with expected errors below 5 % for all parameters. Among the GB approaches, the GBS formulation gives satisfactory results (expected error less than 10 %) for most parameters but completely fails the estimation of parameters TAS, TPF and RSO, consistently with their high RMSEs values. Results obtained with the GBM formulation are more homogeneous but of lower quality, with expected error close to 20 %, for most parameters.

## 3.3. Population-based approach provides robust estimates on real data

Results on simulated data were confirmed by estimations obtained over the peach progeny. When comparing results obtained over different repetitions of the calibration algorithm, GB strategies suffered from a substantial lack of reproducibility.

Fig. 6 shows the relative standard deviation of KHK estimates obtained for the 106 genotypes, with the three considered approaches. The distribution is much wider for GB approaches compared to the PB

<table><caption>**Table 3**

Comparison between parameter estimates and exact values: mean estimate and root mean square error based on 100 simulated datasets.</caption><thead><tr><th rowspan="3">Param.</th><th colspan="5">Genotype-based</th><th colspan="4">Population-based</th></tr><tr><th colspan="3">GBS</th><th colspan="2">GBM</th><th rowspan="2">Param.</th><th rowspan="2">Exact</th><th rowspan="2">Mean estimate</th><th rowspan="2">RMSE</th></tr><tr><th>Exact</th><th>Mean estimate</th><th>RMSE</th><th>Mean estimate</th><th>RMSE</th></tr></thead><tbody><tr><td>LHx</td><td>0.0302</td><td>0.0322</td><td>0.1250</td><td>0.0316</td><td>0.1774</td><td>β<sub>LHx</sub></td><td>0.0302</td><td>0.0300</td><td>0.0071</td></tr><tr><td>RAI</td><td>0.0111</td><td>0.0110</td><td>0.0684</td><td>0.0102</td><td>0.1689</td><td>β<sub>RAI</sub></td><td>0.0111</td><td>0.0111</td><td>0.0071</td></tr><tr><td>KFk</td><td>12.182</td><td>12.219</td><td>0.0270</td><td>11.579</td><td>0.0924</td><td>β<sub>KFk</sub></td><td>12.182</td><td>12.121</td><td>0.0071</td></tr><tr><td>KHK</td><td>20.085</td><td>20.264</td><td>0.0941</td><td>20.594</td><td>0.0876</td><td>β<sub>KHK</sub></td><td>20.085</td><td>19.985</td><td>0.0071</td></tr><tr><td>OCp</td><td>665.14</td><td>735.00</td><td>0.2535</td><td>614.85</td><td>0.1225</td><td>β<sub>OCp</sub></td><td>665.14</td><td>661.81</td><td>0.0071</td></tr><tr><td>RSS</td><td>33.115</td><td>34.212</td><td>0.2004</td><td>30.249</td><td>0.1654</td><td>β<sub>RSS</sub></td><td>33.115</td><td>32.950</td><td>0.0071</td></tr><tr><td>TAS</td><td>20.085</td><td>25.948</td><td>0.4684</td><td>20.244</td><td>0.1661</td><td>β<sub>TAS</sub></td><td>20.085</td><td>19.985</td><td>0.0071</td></tr><tr><td>TPF</td><td>12.182</td><td>14.310</td><td>1.2457</td><td>11.945</td><td>0.1484</td><td>β<sub>TPF</sub></td><td>12.182</td><td>12.121</td><td>0.0071</td></tr><tr><td>RSO</td><td>1.0101</td><td>0.9067</td><td>0.5032</td><td>0.8604</td><td>0.3126</td><td>β<sub>RSO</sub></td><td>1.0101</td><td>1.0092</td><td>0.0011</td></tr><tr><td>σ<sup>2</sup></td><td>0.05</td><td>0.0466</td><td>0.0034</td><td>-</td><td>-</td><td>Γ<sub>LHx</sub></td><td>0.0030</td><td>0.0030</td><td>0.0002</td></tr><tr><td>σ<sup>2</sup><sub>1</sub></td><td>0.05</td><td>-</td><td>-</td><td>0.0979</td><td>0.1383</td><td>Γ<sub>RAI</sub></td><td>0.0011</td><td>0.0011</td><td>0.0016</td></tr><tr><td>σ<sup>2</sup><sub>2</sub></td><td>0.05</td><td>-</td><td>-</td><td>0.1229</td><td>0.1726</td><td>Γ<sub>KFk</sub></td><td>1.2182</td><td>1.2288</td><td>0.0012</td></tr><tr><td>σ<sup>2</sup><sub>3</sub></td><td>0.05</td><td>-</td><td>-</td><td>0.1054</td><td>0.1488</td><td>Γ<sub>KHK</sub></td><td>2.0086</td><td>2.0158</td><td>0.0005</td></tr><tr><td>σ<sup>2</sup><sub>4</sub></td><td>0.05</td><td>-</td><td>-</td><td>0.1015</td><td>0.1434</td><td>Γ<sub>OCp</sub></td><td>66.514</td><td>64.461</td><td>0.0044</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>Γ<sub>RSS</sub></td><td>3.3115</td><td>3.3064</td><td>0.0002</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>Γ<sub>TAS</sub></td><td>2.0086</td><td>2.0156</td><td>0.0005</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>Γ<sub>TPF</sub></td><td>1.2182</td><td>1.2287</td><td>0.0012</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>Γ<sub>RSO</sub></td><td>0.1010</td><td>0.1011</td><td>0.0001</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>σ<sup>2</sup></td><td>0.0500</td><td>0.0697</td><td>0.0640</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Expected error (%) on the individual estimates of GBS/GA, GBM/NSGA-II and PB/SAEM approach over 100 simulated datasets. Simulated datasets consisted of 100 genotypes for GB approach (20 repetitions) and 100 populations of 100 genotypes for BP approach (5 repetitions). For PB approach, values are averaged for each genotype over all populations in order to obtain the same dataset size as in GB approach.</figcaption></figure>

strategy. The robustness of GBM estimations, in particular, strongly depended on the considered genotype, with standard deviation values ranging from zero to 70 %. More generally, large variations in the estimated values were observed for several parameters (see Figure S3 of the supplementary information), although these parameters did not systematically correspond to high RMSE values when estimated over a simulated dataset (see Table 3: simulated data). This suggest a lack of biological information (practical identifiability) for some of the genotypes, which adds to the estimation accuracy of GB algorithms. In agreement with the simulation study, the PB approach proved much more stable across repetitions, independently of the genotype. Table 4 reports the confidence intervals of the vector of parameters $\theta_{PB}$, characterizing the inter-individual distribution. For all estimated parameters, the confidence intervals were remarkably narrow, confirming the ability of PB approach to robustly estimate both mean and individual parameter values of the population.

### 3.4. *Candidate role of Kfk in low-fructose phenotype is confirmed over the whole peach population*

Quantifying the inter-individual variability of model parameters is an important objective of genetic studies. The comparison of the population parameters computed from individual GB estimates to the $\theta$ values obtained by PB approach (Table 4) showed an overall agreement on the estimation of mean parameters ($\beta$), at least in their magnitude, providing a good overview of the timescales of the system. Conversely, the estimation of parameter variations around the mean value turns out to be more sensitive to the calibration method. Although estimations with PB and GBM approaches shared similar orders of magnitude, quantitative differences were observed in the estimated variance for several parameters. The GBM strategy in particular usually results into a lower parameter variability with respect to the PB approach. The worst results, however, were obtained with the GBS approach. Thus, GBS approach calibration largely over-estimated the expected variance for
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Estimation variability over repetitions of the calibration algorithms. Relative standard deviation of individual estimates of the parameter KHk for the 106 genotypes using GBS/GA, GBM/NSGA-II (10 repetitions) and PB/ MCMC-SAEM (10 repetitions).</figcaption></figure>

<caption>**Table 4**</caption>

Comparison of parameter values estimated by the PB and GB approaches over the real peach progeny. Mean, standard deviation (SD) and confidence intervals (CI) are given for population parameters $\theta_{PB}$ estimated by PB approach; mean parameter values and variance obtained with GB approaches. Bold values highlight the large variances predicted by GBS approach with respect to the other calibration procedures.

<table><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Mean estimate</th><th colspan="2">PB</th><th rowspan="2">GBS<br>Mean estimate</th><th rowspan="2">GBM<br>Mean estimate</th></tr><tr><th>SD</th><th>95 % CI</th></tr></thead><tbody><tr><td>β<sub>LHx</sub></td><td>0.4914</td><td>0.0093</td><td>[0.4747;<br>0.5108]</td><td>0.45</td><td>0.42</td></tr><tr><td>β<sub>RAI</sub></td><td>0.1075</td><td>0.0005</td><td>[0.1065;<br>0.1086]</td><td>0.028</td><td>0.028</td></tr><tr><td>β<sub>KFK</sub></td><td>34.8449</td><td>0.8078</td><td>[33.0419;<br>36.2398]</td><td>34.7</td><td>45.1</td></tr><tr><td>β<sub>KHR</sub></td><td>206.3626</td><td>2.9419</td><td>[200.8797;<br>212.3122]</td><td>236.0</td><td>251.9</td></tr><tr><td>β<sub>OCp</sub></td><td>774.9177</td><td>4.8343</td><td>[765.3165;<br>784.1586]</td><td>943.9</td><td>1001</td></tr><tr><td>β<sub>RSS</sub></td><td>406.2915</td><td>4.1089</td><td>[398.8225;<br>414.7968]</td><td>449.4</td><td>326.0</td></tr><tr><td>β<sub>TAS</sub></td><td>147.5079</td><td>1.4867</td><td>[144.3786;<br>150.1947]</td><td>184.0</td><td>136.7</td></tr><tr><td>β<sub>TPF</sub></td><td>51.0577</td><td>0.7856</td><td>[49.2564;<br>52.4168]</td><td>25.7</td><td>70.38</td></tr><tr><td>β<sub>RSO</sub></td><td>2.6561</td><td>0.0975</td><td>[2.4496;<br>2.8277]</td><td>2.21</td><td>2.66</td></tr><tr><td>Γ<sub>LHx</sub></td><td>1.9430</td><td>0.0355</td><td>[1.8704;<br>2.0099]</td><td>161.7</td><td>2.07</td></tr><tr><td>Γ<sub>RAI</sub></td><td>2.7355</td><td>0.0208</td><td>[2.6967;<br>2.7764]</td><td>72.2</td><td>2.62</td></tr><tr><td>Γ<sub>KFK</sub></td><td>2.0338</td><td>0.0260</td><td>[1.9884;<br>2.0921]</td><td>2.82</td><td>6.51</td></tr><tr><td>Γ<sub>KHR</sub></td><td>1.9708</td><td>0.0268</td><td>[1.9174;<br>2.0210]</td><td>1.05</td><td>1.04</td></tr><tr><td>Γ<sub>OCp</sub></td><td>1.2225</td><td>0.0019</td><td>[1.2189;<br>1.2266]</td><td>1.04</td><td>1.01</td></tr><tr><td>Γ<sub>RSS</sub></td><td>2.0764</td><td>0.0239</td><td>[2.0278;<br>2.1224]</td><td>33.1</td><td>1.76</td></tr><tr><td>Γ<sub>TAS</sub></td><td>2.2615</td><td>0.0258</td><td>[2.2150;<br>2.3154]</td><td>3.82</td><td>1.76</td></tr><tr><td>Γ<sub>TPF</sub></td><td>2.4144</td><td>0.0360</td><td>[2.3521;<br>2.4936]</td><td>9.93</td><td>1.46</td></tr><tr><td>Γ<sub>RSO</sub></td><td>2.7073</td><td>0.0323</td><td>[2.6427;<br>2.7705]</td><td>247.6</td><td>1.36</td></tr><tr><td>σ<sup>2</sup></td><td>0.9824</td><td>0.0341</td><td>[0.9199;<br>1.0484]</td><td>-</td><td>-</td></tr></tbody></table>

parameters LHx, RAI, RSS, TPF, RSO (Cf. bolded values in Table 4). This was probably due to a handful of genotypes for which the estimation process was unstable.

Despite the aforementioned discrepancies, both GB and PB methods were able to capture specific features of the peach population under study. For all three methods, the distribution of *KFk* was clearly bimodal, suggesting the presence of two subpopulations, with contrasting affinity values. This result is in agreement with the prediction of the model (Desnoues et al., 2018) according to which the fructokinase affinity may be responsible for the appearance of two distinct fructose phenotypes, as observed in the peach population. In order to further investigate the relationships between the parameter values and the resulting fructose concentration, we systematically compared the estimations obtained over the two phenotypic populations, with the three calibration methods (Fig. 7). Results supported a significant and robust difference in the estimated values of *KFk* for the two phenotypic groups, corroborating the involvement of the fructokinase affinity in the specification of the fructose content of the fruit. Other parameters emerged as different between the two phenotypic groups when considering the results from the GB calibration methods (see Figure S6 and Figure S7 of the supplementary information), but these results were not supported when using the PB approach. This was the case for example for *KHk* and *TPF*, which were also pointed out as putative mechanisms in the original modeling work by Desnoues et al. (2018), using a GBS approach. These results underline the importance of the choice of the calibration method which may impact the biological interpretation of the results.

## 4. Discussion

In the context of agronomy, models are increasingly used to describe plant growth and metabolism and test their performances in different conditions. Being able to account for intra-specific variability is a key challenge in the field as it can provide useful information for varietal selection and plant breeding. In this paper, a mathematical model of fruit sugar metabolism (Kanso et al., 2020) was calibrated over a population of 106 peach varieties using two contrast strategies that strongly differ in the way data are interpreted and analysed: an individual, genotype-based (GB) approach and a population-based (PB) approach. The former was further decomposed into a single-objective (GBS) and a multi-objective (GBM) formulation, in order to account for possible trade-off among phenotypic traits.

When applied to the model proposed by Kanso et al. (2020), all approaches led to a satisfactory agreement between predictions and data. In term of phenotypic predictions, the GBS formulation generally provided better results, especially concerning the dynamics of sucrose and glucose, the prevailing sugars in peach fruit, closely followed by PB approach. However, when looking at the estimation of the underlying parameter values, the GB strategy suffered from a lack of reproducibility, with several parameters sets giving an equivalent agreement with data. This is a well-known issue of collective fit in system biology: in spite of a good predictive power, parameter estimation may show large uncertainty, in correspondence to sloppy directions of the parameter space (Gutenkunst et al., 2007). Tests over a simulated dataset showed that, in the case of GBS, sloppiness mainly concerned three (out of nine) model parameters, TAS, TPF, RSO. The estimates of these parameters were both inaccurate (see Table 3, Fig. 5 over simulated data) and unstable across repetitions (see Figure S3 of the supplementary information over the peach progeny). The use of a multi-objective (GBM) formulation did not really improve results, leading to an intermediate, but still large variability of estimates. Indeed, although the contribution of individual sugars to the objective function was better represented, the inclusion of 4 additional parameters ($\sigma_i$) reduced the number of data-points available for the parameter estimation. These results agree with previous works on other models which suggested that classical GB formulations are often too simplistic to accurately estimate genotypic parameters, regardless of the algorithm used (Barrasso et al., 2019; Kanso et al., 2020).
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Difference in the estimated model parameter (KFk) between standard ('FRU') and low fructose phenotypes ('SSFRU') and its distribution over 106 genotypes obtained with the three considered methods: GBS, GBM and PB approach.</figcaption></figure>

Constantinescu et al., 2016). Better results were undoubtedly obtained with the population-based approach. The advantage of this approach is that it takes into account the observation of the whole progeny and cast them in the same general framework (common probability distribution). The presence of multiple genotypes with a similar sugar concentration indeed helped to constraint the estimation process, resulting in a high stability and accuracy of individual estimates, on both the simulated (Fig. 5) and the experimental datasets (Fig. 6 and Table 4).

In addition to individual estimates, PB approaches naturally yield information on parameter distribution across genotypes. In particular, it allows for a direct analysis of genetically relevant correlations among parameters, via the estimation of the covariance matrix. To this aim, different covariance structures can be tested, reflecting available information or alternative biological hypotheses. For the peach progeny used in this study, a diagonal covariance matrix, corresponding to the parameter independence, proved to be the best choice according to the AIC criterion (Burnham and Anderson, 2002), compared to a full-rank matrix (full covariance matrix AIC = 2396.1 vs diagonal covariance AIC = 2299.7). This suggests that the genetic control of the biochemical reactions involved in peach sugar metabolism is split over independent regulatory modules (Desnoues et al., 2016, Carreno-Quintero et al., 2013). A diagonal covariance matrix with 10% variance was also assumed in our simulation study. Although this corresponds to a realistic situation, varying these variances would allow us to assess the ability of the algorithms to capture inter-individual differences. In particular, large gamma values could affect the accuracy of individual parameter estimates. When comparing PB results to parameter distribution inferred from individual GB-based estimates, results showed that both strategies overall agreed on the estimation of mean parameter values (Table 4) but the estimation of parameter variance was strongly affected by the presence of outliers in the case of GB approaches. In spite of this limitation, both GB and PB approaches were able to capture specific features of our peach population. The estimation of parameter KFk, corresponding to the affinity of fructokinase, revealed a bimodal distribution, i.e. the observed variability corresponded to the presence of two sub-populations, with different mean values. This result agrees remarkably well with the existence of two distinct phenotypes in the population, with contrasted fructose contents, and with previous analysis (Desnoues et al., 2018) suggesting KFk as the main determinant of such a phenotypic trait.

In light of our results, the PB approach appears a suitable choice for genetic studies, as long as a few conditions are satisfied. First, the advantages of PB strategy with respect to GB approaches may depend on the population size. In Baey et al. (2018), parameter estimation over a population of 34 genotypes gave rise to instabilities, due to convergence issues, and proved sensitive to the parameters setting of the algorithm.

In our case, the large population size ensured a good and stable estimation of parameter values, independently of the starting conditions and the algorithm setting. For the tuning parameters there is no general recommendation in the literature, and in practice one may have to run the algorithm several times to identify reasonable values for these parameters. A future study could include systematic simulation experiments with different population sizes and noise levels. This would provide quantitative insights into the minimum conditions required to ensure the reliability of the estimates and to evaluate the robustness of the method under variable experimental conditions. For populations of small size, GB approaches offer a good alternative in terms of phenotypic prediction but the algorithm has to be repeated several times in order to check the stability of the estimates. In case of problems, efforts should be directed in improving the number of datapoints, in order to further constraint the exploration of parameter space.

Last but not least, the fundamental hypothesis made in this paper for the PB approach, namely that individual genetic parameters are realization of a Gaussian distribution, might not be valid when multi-specific or natural (inbred) populations are considered. In this case, extension of the approach can be considered, with non-Gaussian distribution for the individual parameters, accounting for the presence of sub-populations. However, parameter estimation is more complex in these cases, and other numerical methods may be necessary, such as Bayesian approaches for example.

At term, a reliable estimation of model parameters will offer the possibility to study the genetic control of these parameters, looking for genomic regions (QTL) of interest (Cooper et al., 2009, Yin et al., 2016). This will open new opportunities to decipher the genetic and molecular mechanisms underlying the sugar metabolism in fruit, by combining modelling approaches and lab experimentations. More broadly, such combinations of process-based or crop models with whole-genome approaches are currently used to improve predictions of traits under the influence of GxE interactions (Messina et al., 2018; Reynolds et al., 2021; Jighly et al., 2023).

## Funding

Hussein KANSO was founded by a scholarship of the Lebanese government.

## CRediT authorship contribution statement

Kanso Hussein: Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Baldazzi Valentina: Writing – original draft, Visualization, Validation, Supervision, Software, Methodology, Investigation, Formal
analysis, Data curation, Conceptualization. **Quilot Bénédicte:** Writing – original draft, Visualization, Validation, Supervision, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Baey Charlotte:** Writing – original draft, Validation, Supervision, Software, Methodology, Investigation, Conceptualization. **MEMAH Mohamed-Mahmoud:** Writing – original draft, Supervision, Software, Methodology, Formal analysis, Conceptualization.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper

## Acknowledgements

We would like to thank V. Signoret for her help in maintaining the inter-specific peach progeny. We are grateful to the IE-EMMAH UMR1114 and IE-GAFL UR1052 teams for taking care of the experimental orchard.

## Appendix A. Supporting information

Supplementary data associated with this article can be found in the online version at doi:10.1016/j.eja.2025.127641.

## Data availability

Some data and codes are available on line.

## References

Akdemir, S., 2019. Fritsche-Neto, and Isidro-Sánchez. Multi-objective optimized genomic breed-ing strategies for sustainable food improvement. Heredity 122, 672-683.

Baey, C., Mathieu, A., Jullien, A., Trevezas, S., Cournède, P.-H., 2018. Mixed-effects estimation in dynamic models of plant growth for the assessment of inter-individual variability. J. Agric., Biol. Environ. Stat. 23 (2), 208-232.

Baey, S.Trevezas, Cournède, P.-H., 2016. A non linear mixed effects model of plant growth-and estimation via stochastic variants of the em algorithm. Commun. Stat. - Theory Methods 45 (6), 1643-1669. https://doi.org/10.1080/03610926.2014.930909.

Barrasso, C., Memah, M.-M., Génard, M., Quilot-Turion, B., 2019. Model-based QTL detection is sensitive to slight modifications in model formulation. PLoS ONE 14 (10), e0222764. https://doi.org/10.1371/journal.pone.0222764.

Becker, V., Schilling, M., Bachmann, J., Baumann, U., Raue, A., Maiwald, T., Timmer, J., Klingmüller, U., 2010. Covering a broad dynamic range: information processing at the erythropoietin receptor. Science 328 (5984), 1404-1408.

Bertin, N., Martre, P., Genard, M., Quilot-Turion, B., Salon, C., 2010. Under what circumstances can process-based simulation models link genotype to phenotype for complex traits? case study of fruit and grain quality traits. 02 J. Exp. Bot. 61, 955-967, 02.

Burnham, K.P., Anderson, D.R., 2002. Model selection and multimodel inference: a practical information-theoretic approach. Springer, New York. ISBN 9780387224565.

Carreno-Quintero, N., Bouwmeester, H.J., Keurentjes, J.J., 2013. Genetic analysis of metabolome-phenotype interactions: from model to crop species. Trends Genet. 29 (1), 41-50.

Chou, I.-C., Voit, E.O., 2009. Recent developments in parameter estimation and structure identification of biochemical and genomic systems. Math. Biosci. 219 (2), 57-83.

Coello, C.A.C., Lamont, G.B., Van Veldhuizen, D.A., et al., 2007. Evolutionary algorithms for solving multi-objective problems. 5. Springer.

Conn, A., Gould, N., Toint, P., 1997. A globally convergent lagrangian barrier algorithm for optimization with general inequality constraints and simple bounds. Math. Comput. 66 (217), 261-288.

Conn, A.R., Gould, N.I., Toint, P., 1991. A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds. SIAM J. Numer. Anal. 28 (2), 545-572.

Constantinescu, D., Memmah, M.-M., Vercambre, G., Génard, M., Baldazzi, V., Causse, M., Albert, E., Brunel, B., Valsesia, P., Bertin, N., 2016. Model-Assisted Estimation of the Genetic Variability in Physiological Parameters Related to Tomato Fruit Growth under Contrasted Water Conditions. ISSN 1664462X Front. Plant Sci. 7 (December), 1-17. https://doi.org/10.3389/fpls.2016.01841.

Cooper, M., van Eeuwijk, F. a, Hammer, G.L., Podlich, D.W., Messina, C., 2009. Modeling QTL for complex traits: detection and context for plant breeding. ISSN 1879-0356 Curr. Opin. Plant Biol. 12 (2), 231-240. https://doi.org/10.1016/j.pbi.2009.01.006.

Correndo, A.A., Hefley, T.J., Holzworth, D.P., Ciampitti, I.A., 2021. Revisiting linear regression to test agreement in continuous predicted-observed datasets. Agric. Syst. 192, 103194. https://doi.org/10.1016/j.agsy.2021.103194.

Davidian, M., Giltinan, D.M., 2003. Nonlinear models for repeated measurement data: An overview and update. ISSN 10857117 J. Agric., Biol., Environ. Stat. 8 (4), 387-419. https://doi.org/10.1198/1085711032697.

Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multiobjective genetic algorithm: Nsga-ii. IEEE Trans. Evolut. Comput. 6 (2), 182-197.

Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from incomplete data via the em algorithm. J. R. Stat. Soc.: Ser. B (Methodol. 39 (1), 1-22.

E. Desnoues, Y. Gibon, V. Baldazzi, V. Signoret, M. Génard, and B. Quilot-Turion. Profiling sugar metabolism during fruit development in a peach progeny with different fructose-toglucose ratios. BMC Plant Biology, 14(1):336, Nov 2014. ISSN 1471-2229.

Desnoues, E., Baldazzi, V., Génard, M., Mauroux, J.-B., Lambert, P., Confolent, C., Quilot-Turion, B., 2016. Dynamic qtls for sugars and enzyme activities provide an overview of genetic control of sugar metabolism during peach fruit development. J. Exp. Bot. 67 (11), 3419-3431.

Desnoues, E., Génard, M., Quilot-Turion, B., Baldazzi, V., 2018. A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-toglucose ratio phenotype. Plant J. 94 (4), 685-698.

Donnet, S., Samson, A., 2013. A review on estimation of stochastic differential equations for pharmacokinetic/pharmacodynamic models. Adv. Drug Deliv. Rev. 65 (7), 929-939.

Efron, B., Tibshirani, R., 1993. An introduction to the bootstrap. Chapman & Hall, London, p. 1993.

Goldberg, D.E., 1989. Genet. Algorithms Search. Optim., Mach.

R.N. Gutenkunst, J.J. Waterfall, F.P. Casey, K.S. Brown, C.R. Myers, and J.P. Sethna. Universally sloppy parameter sensitivities in systems biology models. PLoS computational biology, 3(10):1871-1878, oct 2007. ISSN 1553-7358. doi: 10.1371/journal.pcbi.00030189.

Hass, H., Masson, K., Wohlgemuth, S., Paragas, V., Allen, J.E., Sevecka, M., Pace, E., Timmer, J., Stelling, J., MacBeath, G., et al., 2017. Predicting ligand-dependent tumors from multidimensional signaling features. NPJ Syst. Biol. Appl. 3 (1), 1-15.

W.K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. Oxford University Press, 1970.

Hightly, A., Weeks, A., Christy, B., J O'Leary, G., Kant, S., Aggarwal, R., Hessel, D., Forrest, K.L., Technow, F., Tibbits, J.F., Totir, R., Spangenberg, G.C., Hayden, M.J., Munkvold, J., Daetwyler, H.D., 2023. Integrating biophysical crop growth models and whole genome prediction for their mutual benefit: a case study in wheat phenology. J. Exp. Bot. 74 (15), 4415-4426.

Kanso, H., Quilot-Turion, B., Memah, M.-M., Bernard, O., Gouzé, J.-L., Baldazzi, V., 2020. Reducing a model of sugar metabolism in peach to catch different patterns among genotypes. Math. Biosci. 321, 108321.

Kuhn, E., Lavielle, M., 2004. Coupling a stochastic approximation version of em with an mcmc procedure. ESAIM: Probab. Stat. 8, 115-131.

E. Kuhn and M. Lavielle. Maximum likelihood estimation in nonlinear mixed effects models. Computational statistics & data analysis, 49(4):1020-1038, 2005.

Lamsal, S.M., Welch, S.M., Jones, J.W., Boote, K.J., Asebedo, A., Crain, J., Wang, X., Boyer, W., Giri, A., Frink, E., Xu, X., Gundy, G., Ou, J., Arachchige, P.G., 2017. Efficient crop model parameter estimation and site characterization using large breeding trial data sets. Agric. Syst. 157, 170-184.

Lamsal, S.M., Welch, J.W., White, K.R.Thorp, Bello, N.M., 2018. Estimating parametric phenotypes that determine anthesis date in Zea mays: Challenges in combining ecophysiological models with genetics. PLoS ONE 13 (4).

Laperche, A., Devienne-Barret, F., Maury, O., Le Gouis, J., Ney, B., 2006. A simplified conceptual model of carbon/nitrogen functioning for qtl analysis of winter wheat adaptation to nitrogen deficiency. Theor. Appl. Genet. 113 (6), 1131-1146.

Maiwald, T., Hass, H., Steiert, B., Vanlier, J., Engesser, R., Raue, A., Kipkeew, F., Bock, H. H., Kaschek, D., Kreutz, C., et al., 2016. Driving the model to its limit: profile likelihood based model reduction. PloS One 11 (9), e0162366.

Martre, P., Quilot-Turion, B., Luquet, D., Memmah, M.-M.O.-S., Chenu, K., Debaeke, P., 2015. Chapter 14 - model-assisted phenotyping and ideotype design (Academic Press, San Diego, second edition edition). V. O. Sadras D. F. Calderini, Ed., Crop Physiol. (Second Ed.) 349-373 (Academic Press, San Diego, second edition edition).

Messina, C.D., Jones, J.W., Boote, K.J., Vallejos, C.E., 2006. A Gene-Based Model to Simulate Soybean Development and Yield Responses to Environment. ISSN 1435-0653 Crop Sci. 46 (1), 456. https://doi.org/10.2135/cropsci2005.04-0372.

Messina, C.D., Technow, F., Tang, T., Totir, R., Gho, C., Cooper, M., 2018. Leveraging biological insight and environmental variation to improve phenotypic prediction: Integrating crop growth models (CGM) with whole genome prediction (WGP). European J. Agron. 100, 151-162.

Nakagawa, H., Yamagishi, J., Miyamoto, N., Motoyama, M., Yano, M., Nemoto, K., 2005. Flowering response of rice to photoperiod and temperature: a qtl analysis using a phenological model. Theor. Appl. Genet. 110 (4), 778-786.

Owen, J., Fiedler-Kelly, J., 2014. Introduction to Population Pharmacokinetic / Pharmacodynamic Analysis with Nonlinear Mixed Effects Models. Wiley.

Perez, R.P., Pallas, B., Le Moguédec, G., Rey, H., Griffon, S., Caliman, J.-P., Costes, E., Dauzat, J., 2016. Integrating mixed-effect models into an architectural plant model to simulate inter- and intra-progeny variability: a case study on oil palm (Elaeis guineensis Jacq. ISSN 0022-0957 J. Exp. Bot. 67 (15), 4507-4521. https://doi.org/10.1093/jxb/erw203.

Píneiro, G., Perelman, S., Guerschman, J.P., Paruelo, J.M., 2008. How to evaluate models: Observed vs. predicted or predicted vs. observed? Ecol. Model. 216 (3-4), 316-322. https://doi.org/10.1016/j.ecolmodel.2008.05.006.
Prudent, M., Lecomte, A., Bouchet, J.-P., Bertin, N., Causse, M., Génard, M., 2011. Combining Ecophysiological Modelling and Quantitative Trait Locus Analysis to Identify key Elementary Processes Underlying Tomato Fruit Sugar Concentration. J. Exp. Bot. 62, 907–919.

Quilot, B., Génard, M., Kervella, J., Lescourret, F., 2004. Analysis of genotypic variation in fruit flesh total sugar content via an ecophysiological model applied to peach. Theor. Appl. Genet. 109 (2), 440–449.

Quilot, B., Kervella, J., Génard, M., Lescourret, F., 2005. Analysing the genetic control of peach fruit quality through an ecophysiological model combined with a qtl approach. J. Exp. Bot. 56 (422), 3083–3092.

Ravi Kumar, S., Hammer, G.L., Broad, I., Harland, P., McLean, G., 2009. Modelling environmental effects on phenology and canopy development of diverse sorghum genotypes. ISSN 0378-4290 Field Crops Res. 111 (1), 157–165. https://doi.org/10.1016/j.fcr.2008.11.010.

Reymond, M., Muller, B., Leonardi, A., Charcosset, A., Tardieu, F., 2003. Combining quantitative trait loci analysis and an ecophysiological model to analyze the genetic variability of the responses of maize leaf growth to temperature and water deficit. Plant Physiol. 131 (2), 664–675.

Reynolds, M.P., Atkin, O.K., Bennett, M., Cooper, M., Dodd, I.C., Foulkes, M.J., Frohberg, C., Hammer, G., Henderson, I.R., Huang, B., Korzun, V., McCouch, S.R., Messina, C.D., Pogson, B.J., Slafer, G.A., Taylor, N.L., Wittich, P.E., 2021. Addressing Research Bottlenecks to Crop Productivity. Trends Plant Sci. 26 (6), 607–630.

Robbins, H., Monro, S., 1951. A stochastic approximation method. Ann. Math. Stat. 400–407.

Snowden, T.J., van der Graaf, P.H., Tindall, M.J., 2017. Methods of model reduction for largescale biological systems: a survey of current methods and trends. Bull. Math. Biol. 79 (7), 1449–1486.

Yin, X., Chasalow, S.D., Dourleijn, C.J., Stam, P., Kropff, M.J., 2000. Coupling estimated effects of qtls for physiological traits to a crop growth model: predicting yield variation among recombinant inbred lines in barley. Heredity 85 (6), 539–549.

Yin, X., Struik, P.C., van Eeuwijk, F.A., Stam, P., Tang, J., 2005. QTL analysis and QTL-based prediction of flowering phenology in recombinant inbred lines of barley. ISSN 0022-0957 J. Exp. Bot. 56 (413), 967–976. https://doi.org/10.1093/jxb/eri090.

Yin, X., Struik, P.C., Gu, J., Wang, H., 2016. Modelling QTL-Trait-Crop Relationships: Past Experiences and Future Prospects. In Crop Systems Biology. Springer International Publishing, Cham, pp. 193–218.