Analysis of root growth from a phenotyping data set using a density-based model

Major research efforts are targeting the improved performance of root systems for more efficient use of water and nutrients by crops. However, characterizing root system architecture (RSA) is challenging, because roots are difficult objects to observe and analyse. A model-based analysis of RSA traits from phenotyping image data is presented. The model can suc-cessfully back-calculate growth parameters without the need to measure individual roots. The mathematical model uses partial differential equations to describe root system development. Methods based on kernel estimators were used to quantify root density distributions from experimental image data, and different optimization approaches to parameterize the model were tested. The model was tested on root images of a set of 89 Brassica rapa L. individuals of the same genotype grown for 14 d after sowing on blue filter paper. Optimized root growth parameters enabled the final (modelled) length of the main root axes to be matched within 1% of their mean values observed in experiments. Parameterized values for elongation rates were within ±4% of the values measured directly on images. Future work should investigate the time dependency of growth parameters using time-lapse image data. The approach is a potentially powerful quantitative technique for identifying crop genotypes with more efficient root systems, using (even incomplete) data from high-throughput phenotyping systems.


Introduction
The identification of crop genotypes that can efficiently capture soil water and nutrients has become a major focus of research (White et al., 2013b). Recent work has demonstrated, for example, that crops with deeper root systems may have better water uptake efficiency, while shallower root systems often improve phosphorus (P) uptake (Lynch and Wojciechowski, 2015). However, plant roots serve many other functions. They must acquire other essential mineral elements (White et al., 2013a) and provide mechanical support to aerial stems (Dupuy et al., 2005a), while delivering nutrients for soil microorganisms through exudation and biomass turnover (van der Putten et al., 2013). They are able to achieve such functions in a vast heterogeneity of soil profiles and a wide range of soil physical, chemical, biotic and abiotic environments. To grow and thrive in such diverse soil environments, plants adapt their root development to explore patches of nutrients (Wang et al., 2006), to explore cavities and pores to grow through hard soils (Stirzaker et al., 1996), to create specialized structures (Hodge et al., 2009) and to establish mutualistic associations with a range of microorganisms to access nutrients that are not available to them (Richardson et al., 2000). Understanding what makes a root system more efficient is therefore a multifaceted problem, with many aspects remaining poorly understood.
To improve the root systems of crops, it is essential to find ways to observe and characterize their growth patterns. Sophisticated imaging techniques are now commonly used in root studies. These include magnetic resonance tomography, X-ray microtomography, neutron radiography for detailed analysis of root systems in situ, and visible light imaging in rhizotrons and minirhizotrons (Clark et al., 2011). Detailed visualization of root growth processes also reveals complex root geometries (Zhu et al., 2011). These complex root branching structures cannot be described easily as quantitative traits and, as a consequence, this represents a serious limitation to breeding for improved root system architectures.
Since root architectures are complex dynamic systems that have underlying patterns of growth, models may be used to interpret and analyse such root growth data (de Dorlodot et al., 2007). Complex root architecture arises from elementary processes of root elongation and branching. Models allow us to derive complex patterns from simple growth rules, and inverse modelling can be used to estimate what combination of elongation rates and branching interval parameters may have given rise to any given root system structure that obeys the underlying model rules. Models may also help to overcome issues of incomplete image data sets caused by limitations of imaging systems (Wells et al., 2012). However, to estimate model parameters from data and derive quantitative traits that could be used for breeding and improvement of crop root systems, efficient optimization techniques must be developed. Previous studies have shown that the problem is challenging, because models that are too simple may be biased or limited to early developmental stages (Dupuy et al., 2005b) and optimization of models that are too complex can overfit the data and usually lead to multiple sets of parameters that are equally optimal (Pedersen et al., 2010;Garré et al., 2012).
In this study, a framework to extract root growth parameters from sets of images of root systems is presented. Our approach is to develop an optimization procedure to backcalculate the parameter values for elongation, gravitropic and branching rates that would give rise to an observed root density distribution. The framework is based on a set of simple time-delay partial differential equations that describe the evolution of root length density distribution. Density estimation methods that are specific to the hierarchical root system structure have been developed so that models can be efficiently fitted to experimental data. The performance of a range of different optimization techniques was assessed based on the ability to retrieve growth parameters from simulated data. The methods were then applied to analyse the growth of Brassica rapa L. subsp. trilocularis cv. R-o-18 roots.

Materials and methods
Root phenotyping data Data were derived from a previous set of experiments (Adu et al., 2014) in which seedlings of B. rapa L. subsp. trilocularis cv. R-o-18 were grown under standardized conditions on seed germination paper for up to 14 d after sowing and were scanned in five experimental trials using A4 CanoScan 5600F scanners (Canon UK Ltd, Reigate, UK). The young B. rapa root system skeleton consisted of a single primary root and a number of first-order lateral branches (Fig. 1A). The spatial position and angular orientation of roots were extracted from the segmented images of 89 plant root systems obtained via the semi-automated image analysis software SmartRoot (Lobet et al., 2011), which provided length measurements and topological information about the root system. In the source image, after root tracing, root systems are divided by root trajectories that could comprise a finite number of line segments. Each line segment may be of variable length and is represented as two successive points described with their co-ordinates on the Cartesian plane. In order to build density distributions, the length, mid-point, and angle of each root segment were calculated. In total for the 89 plants, primary roots were made up of 22 290 individual line segments and lateral roots of 196 055 segments.
Elongation rates for the primary root and lateral roots were estimated directly from the total length and number of primary and lateral roots for each seedling using an established formula (Hackett and Rose, 1972). The elongation rate for primary root e (0) (cm d -1 ) was determined as the total root length l (0) (cm) of the primary root divided by 14 d, Description of the root system using density distribution functions Estimation of model parameters through optimization is based on the ability to quantify differences between the experimental data and the model output. In order to quantify and then minimize such differences, a suitable mathematical representation of the root system density distribution functions must be derived. Here, the description of the root system, both for the model and for the experimental data, used the notion of generalized densities following the concepts proposed by Dupuy et al. (2010a). Three fundamental quantities were used to represent the root system: the root tip density ρ a , the root length density ρ l , and the root branching density ρ b . The root tip density ρ a (cm -3 ) indicates the number of root tips per unit volume (Fig. 1B). The root length density ρ l (cm -2 ) expresses the overall length of roots per unit volume (Fig. 1C). The number of connections between roots (i.e. branching points) per unit volume is represented by the root branching density ρ b (cm -3 ) (Fig. 1D). The root systems were grown on a flat surface forming a 2-dimensional structure. However, to describe both the spatial and angular distribution of roots, densities were defined in a generalized 3-dimensional space of R 3 , so that for each time t the function ρ a (t,x,y,α) defines the density of root tips at a horizontal position x and a vertical position y from the seed as well as in a direction of growth given by the root angle α, defined as the angle from the vertical axis. The branching order of roots is also defined, such that the primary root is of branching order 0 and laterals are of branching order 1. Root density distribution functions were also labelled according to their branching order with ρ ( ) , a i ρ ( ) , l i and ρ b ( ) i denoting root tip, root length and root branching density of the roots of the ith order, respectively. In the following sections, upper indices will be omitted when an equation is valid for any given branching order.
Time-delay density-based root growth model for describing root density distributions Root density distribution functions vary with time as a result of growth processes where root tips are created through branching and move as a result of elongation in a direction determined by gravitropism. A mathematical model adapted from the previously published model of Dupuy et al. (2010a) must therefore be defined to describe the changes in root densities as a function of these processes. Since the root system architecture is defined on a continuum, such models are derived using the concept of conservation of the root quantity in a volume Ω of R 3 . For any arbitrary representative volume Ω of R 3 which is enclosed by a surface ∂Ω, the flux of root tips across the surface ∂Ω can be expressed as a 3-dimensional vector F. In order to find a suitable expression for the flux F, one must consider that roots at a given location (x,y) and growth direction α move only in one direction u = (sinα,cosα). The velocity of root tips is therefore (esinα,ecosα) where e (cm d -1 ) is the elongation rate. However, at the same time, these same roots also change orientation as a result of gravitropism, which is assumed to be linearly related to the root angle based on semi-quantitative observations. The gravitropic rate is expressed as -gα, where g (d -1 ) is the gravitropic rate constant and the negative sign indicates the vertically downward direction. The flux of root tips is therefore expressed as F=ρ a (esinα,ecosα,-gα). It is helpful to write the flux vector in the form: Since κ is the change in root angle divided by the change in root length (dα/dl), it is also the local curvature of the root (Dupuy, 2011). This relationship expresses the relationship between the elongation rate and the gravitropic rate on the morphology of the root system: both gravitropism and elongation are independent parameters, but it is the ratio between the two that determines the shape of the root system. The expression of the velocity field here differs slightly from previous papers (Dupuy et al., 2010a), because the angle α is defined from the vertical axis instead of the horizontal axis. Solutions to the problems, however, are unchanged. Thus, the conservation of the number of root tips for a given branching order i in an arbitrary representative volume Ω implies: where n is the unit normal vector to the surface element dS, and b (cm -2 d -1 ) is the volumetric branching rate. Herein, the model suggested by Dupuy et al. (2010b) is extended to incorporate a time delay in the emergence of lateral roots, which was experimentally observed in the branching patterns of B. rapa (Adu et al., 2014). Therefore, (2) and the source term b (i) for a given branching order i, using notations with explicit branching orders, is then expressed as a function of the root tip density of the previous branching order i-1:  roots of order i are initiated with a time delay denoted by T (i) . In our study, only one primary root per seedling (branching order i=0) and first-order lateral roots (branching order i=1) emerged during the experiment in the examined B. rapa root systems (i.e. i ∈{ , }). 0 1 Hence, the constant b r ( ) 1 1 − or b r ( ) 0 refers to the branching rate of the primary root (branching order i=0) and b a ( ) 1 signifies the branching angle at which the first-order laterals emerge, with a time delay given by T. The time delay T (instead of T (1) for simplification purposes) was experimentally estimated from Adu et al. (2014). Root length density ρ l and root branching density ρ b are then derived from the root tip density ρ a at a particular position X=(x,y,α) using the following relationships: ρ ρ l e t = ∫ a d and ρ b b t = ∫ d .
Estimation of root tip density and root length density from root image data sets It was necessary to quantify images of roots into experimental density distributions for comparison with the model output.
Kernel density approach: A common approach to determine root densities consists of counting the number or length of roots falling in each cell of a spatial grid (Scott and Sain, 2005). The method is equivalent to making a histogram from a multidimensional data set. However, this method suffers from various limitations. Multidimensional histograms are less accurate because the method does not account for differences in position within a cell, the density estimates obtained are discontinuous and it is very difficult to derive an optimal grid, because such a process must consider not only the resolution of the grid but also the position and orientation of the grid in the multidimensional space. Herein, density distribution functions instead were derived using kernel density estimation. Kernel density estimation can be seen as a generalization of the histogram approach with additional theoretical foundation. Kernel-based estimators differ from histograms in the following ways: (i) A cell is viewed as a function, called kernel function, with nonzero values defined only in the neighbourhood of the centre of the kernel. For example, Gaussian or boxcar functions are classically used to describe the values as a function of the distance from the kernel centre.
(ii) Kernel functions are not arranged in a grid but are adjusted so that they can be centred on each point of the data. (iii) The sum of all the kernel functions provides the estimate for the density distribution function.
A histogram density estimate is therefore a boxcar kernel estimate for which the centres of the kernels are placed at the position of the nearest cell of a grid. The method presented here was adapted not only to include the effect of changes in lengths of root segments but also to account for the specific structure of the spatial distribution of points derived from the digitization of the root system. SmartRoot ( Fig. 2A) was used to produce the list of oriented root segments (Adu et al., 2014). Thus, only the spatial distribution of root segments was used in this study, not the connectivity between them. Even though SmartRoot was used to obtain root segments, other software for image analysis could be used to provide similar data. The list of oriented root segments was then used to estimate a length density distribution function (Fig. 2B). In the experimental set-up, seedlings did not have an identical position in the image. Therefore, co-ordinates of all root segments from one image were translated by a constant vector so that all root systems were aligned to where all the seeds were located, at a unique position The position of the centre of the ith root segment of length l i is denoted by X i (Fig. 2B). Then, the estimator of the root length density ρ l  Χ ( ) is constructed using the following formula: Here, W (cm -2 ) is a kernel function and m indicates the total number of root segments. The kernel function W is usually a positive, even, unimodal real-valued function with a maximum at zero and integrates to 1 (i.e.  For this study, the classical Gaussian kernel was chosen and defined as (Scott and Sain, 2005): where |Σ| defines the determinant of the positive definite matrix Σ, which is the covariance matrix in the d-dimensional generalized space. Since the variables x,y,α are independent and therefore uncorrelated, Σ is a non-zero 3 × 3 diagonal of the form: where k is the scaling factor of the kernel density estimator and s x , s y , s α are normalization factors which account for the different grid spacing along the three dimensions x,y,α. Since k determines the smoothness of the density distribution function, we will also call it the spatial resolution. The density estimation depends on the parameters k, s x , s y , and s α of the kernel function W (Hwang and Lay, 1994).
Optimal scaling factor for the kernel density approach: In order to match the kernel density function to the experimental data, an optimal value for the scaling factor (k opt ) must be found. The optimal scaling factor corresponds to the trade-off between overfitting of the data, which is obtained when values of k are too small, and a coarse representation of the data, when values of k are too large. This optimum was determined using a method called pseudo-loglikelihood cross-validation (Turlach, 1993). The likelihood is a statistical metric that indicates how well an estimate agrees with the data. It is calculated as the product of the probabilities of each of the model predictions on the data sets. The logarithm of the likelihood is more commonly used in numerical applications because its expression is usually simpler. Since the scaling factor k depends on the data, there is no independent observation to determine the likelihood value. One way to overcome this limitation is to use the Leave-One-Out cross-validation method denoted by CV LOO .
In the Leave-One-Out cross-validation method, the density is estimated on all but one of the data points (j in Equation 7) and the log likelihood is determined on that missing point. The calculation is repeated for each point in the data set so that a global estimation of the log likelihood is obtained: where m denotes the total number of mid-points in a root system. Unfortunately, this method does not differentiate between root segments belonging to the same root and root segments of two different roots. In order to address this problem, it is important to adopt a different approach to cross-validation. First, it is useful to use a different numbering for the root segments so that indices of roots and their segments appear explicitly in the summation. We therefore denote X ij and l ij the co-ordinates and length of the jth root segment on the ith root: where n is the number of roots and m i is the number of segments on the ith root, so that Traced roots show a high density of nodes on a single root but low point density between roots. To include this point cloud structure, data points were subgrouped with respect to the roots to which they belong. A V-fold cross-validation (CV V-fold ) method was used to determine k opt . In a V-fold cross-validation (Geisser, 1975), entire roots are omitted ( r i ≠ in Equation 9) for the determination of the pseudo-log-likelihood, and the cross-validation function is expressed as: The optimal value of the scaling factor k opt maximizes the The size of the elements of the finite volume grid was adjusted to correspond to the scale factor k identified by the V-fold crossvalidation. In this setting, an increase in the grid resolution cannot lead to an improved fit since the resolution of the experimental data cannot be improved further.

Estimation of density-based model root growth parameters using optimization and comparison with the experimental data set
The previous steps have established a mechanism to determine a root density distribution from root experimental data. In the next step, the objective was to determine the values of the root growth parameters that best reproduce the root length density distribution evaluated from the experimental data for roots of each branching order i (Equation 4). Experiments took place in a controlled environment, for a short period of time (14 d), on seedlings that had already germinated for 3 d. Therefore, growth parameters were assumed to be constant. Differences between simulated density ( ρ l i ( ) ) and the density distribution function estimated from data ( ) and divided by the volume of the domain L x L y L α . For roots at branching order i, the error E (i) was determined as: The vector of optimal model parameters θ θ opt ( ) i was defined by recurrence: This approach is exact because root length density at branching order i is independent of growth parameters for branching order j, if i<j. Since models are solved numerically, optimization is computationally intensive. However, solving at each branching order separately allows the problem to remain computationally tractable. Root branching order can be derived experimentally using either time-lapse data about root emergence or using classes of root diameter as a proxy for root branching order or elongation rate.
To determine the most efficient non-linear optimization algorithm, four different optimization methods were tested: Powell's method, the Conjugate-Gradient method, the Broyden-Fletcher-Goldfarb-Shanno method (BFGS), and the Nelder-Mead method (Nelder and Mead, 1965).The root density at the start of the experiment was chosen as the initial root distribution for the model. When the true parameters are known, the quality Q of the fit can be defined as the ratio of the estimated parameter θ to the original model parameter θ for the simulation so that optimal parameters were identified when The calculated elongation rates were also compared with those derived by direct estimation based on the increase in root length density for each order of root: Since it is not usually possible to identify root tips from images of plants grown in soil, the number of root tips ρ α a V dxdyd ∫∫∫ was assumed to be constant for primary roots. This value is often known and consistent in seedlings of a given plant genotype.

Sensitivity analysis and confidence intervals for fitted densitybased model root growth parameters of B. rapa
The method described in the previous section provides a single set of optimal root growth parameters that best reproduce the root length density distributions measured experimentally. Further analyses are required to determine how reliable the estimates of these growth parameters are. Therefore, sensitivity analysis was performed on model parameters. Sensitivity was assessed based on a 1% increase in each model base parameter value. The Model Elasticity Value (MEV) of the output error, which is specified as the percentage change in the cost function E (Equation 11) per percentage change in the model parameter p, was determined: A bootstrap method was used to evaluate the confidence intervals of model parameters. Since the spatial distribution of mid-points was highly structured, data sets were re-sampled randomly following the V-fold sampling approach described previously. For each bootstrap sample, the model was parameterized and parameter values were recorded. For each parameter, the bootstrap Standard Error (SE) was defined as: where R is the total number of replications, p r denotes the estimated parameter at replication r and p − is the mean parameter value. Then, 95% confidence intervals (CIs) were evaluated together with the coefficient of variation (CoV) for each of the model parameters.
Analysis of the performance of the optimization process Two types of analyses of the optimization process were performed. First, the performance and testing of optimization algorithms was carried out using root growth simulation data (Fig. 3A) to create a target root length density distribution function resulting from known model parameters. With the simulated data, it was possible to quantify directly the goodness of the fit of model parameters.
Optimal model parameters were determined using various nonlinear optimization algorithms. Once the performance of the algorithms was assessed, the algorithm with the faster convergence rate and improved accuracy was used for analysis of real experimental data. In a second step, the optimization algorithm was applied to real data with the experimental root density at t=0 chosen as the initial root distribution for the simulations. Estimation of model parameters from data, however, includes the additional kernel density estimation step (Fig. 3B). The optimization process then follows the same set of steps. An optimization algorithm (Fig. 3C) determines a set of candidate model parameters that are used to simulate the root system (Fig. 3A), and the results of these simulations are compared with the target root system using the cost function (Fig. 3D). The optimization procedure ( Fig. 3C and Fig. 3D) is iterated until convergence criteria are met. Solutions of the model were obtained numerically using the upwind finite volume method presented in Dupuy et al. (2010b). Numerical simulations were carried out on a 3-dimensional domain The size of the domain was determined to match the experimental set-up presented in Adu et al.
(2014), with L x =15 (cm), L y =22 (cm), and . L α π = ( ) 3 2 rad Spatial resolution of the grid was determined based on a scaling factor k such that dx=k, dy=k, dα=gπk/e, and dt=Ck/e where the Courant-Friedrichs-Lewy number C is set to 0.5 (Courant et al., 1967). The algorithms were implemented in the python SciPy library (http://www.scipy.org/) using a personal computer of 3.1 GHz CPU (IntelCore i5-2400 CPU @ 3.1 GHz) and 4 Gb RAM, comparing the target root length distribution obtained for specific user-defined values of model parameters with the length distribution derived from the numerical model at each grid point (Fig. 3). The relevant codes and experimental data can be found at http://www.archiroot. org.uk.

Kernel density estimates: V-fold cross-validation performed better than Leave-One-Out cross-validation
Density estimations based on V-fold cross-validation and Leave-One-Out (LOO) cross-validation were obtained on the root length density both of primary root and of lateral roots. In each case, comparisons between the two methods were performed visually. The optimal scaling factor k' opt at which the function CV LOO (k) attained its maximum (Fig. 4A) was equal to 0.23 for primary and 0.20 for lateral roots. It produced patchy density distributions for both primary root and lateral roots (Fig. 4B). Estimation of total root length with both methods was very close to those measured experimentally, with values of 15.77 cm and 15.68 cm, respectively, obtained by LOO and V-fold cross-validation for the primary root, and values of 106.45 cm and 105.56 cm, respectively, for lateral roots. The measured total primary root length was 15.54 cm and the total lateral root length was 105.26 cm. The close match between model-estimated and measured total lengths showed only minor edge effects at the boundary of the domain. When the CV V-fold method was employed, the optimal scaling factor k opt was equal to 0.9 for primary (Fig. 4A) and 0.57 for lateral roots. The method showed better compromise Fig. 3. Diagrammatic representation of the suggested method for data-driven testing and estimation of the model parameters. The target root system architecture (RSA) is derived through the available data sets (experimental data) or from models for which model parameters are known (simulated data). Testing consists of using simulation with known parameters (A) to create a target root length density distribution function. When the target density distribution results from simulation, it is possible to compare directly the estimated model parameters to target parameters. Estimation of model parameters from data uses a density estimation step (B). The optimization process then follows the same set of steps. An optimization algorithm (C) determines a set of candidate model parameters that are used to simulate the root system (A), and the results of these simulations are compared with the target root system using a cost function (D). The optimization procedures (C and D) are iterated until a convergence criterion is met. Optimal model parameters were determined using non-linear optimization algorithms. The most efficient method which was proven to minimize the discrepancies between the target and the estimated RSA was then applied to experimental data in order to extract root growth parameters. Fig. 4. Comparative performance analysis of Leave-One-Out and V-fold cross-validation methods on the optimal scaling factor k of the kernel estimator. (A) For primary root, employing the Leave-One-Out grouping method, the cross-validation function CV(k) attains its maximum value at k' opt (dashed), while sampling based on the V-fold method yields the maximum CV(k) at k opt (plain). (B) Visualization of the root length distribution with scaling factor k' opt for the kernel estimator (Leave-One-Out method). (C) Visualization of the root length distribution with scaling factor k opt for the kernel estimator (V-fold method). The color map indicates the root length density and units of the colour bar are cm −1 .
between the bias and variance of the estimator, with density maps being smooth, yet not oversmoothed (Fig. 4C).

Factors that influence density-based model parameterization
Initially, optimization algorithms were applied to simulated density data with known growth parameters. The results showed that the direct estimation of the elongation rate e calculated from the change in the root length over time (Equation 13) was not accurate, and direct estimates of elongation rates decreased over time. The convergence to the true model parameter was increased with improved grid resolution (data not shown), but root fluxes at the boundaries of the finite volume grid prevented accurate estimation of the elongation rate as roots grew out of the grid. Simplex-based optimization techniques, such as Nelder-Mead, performed better than gradientbased methods, such as Conjugate-Gradient or BFGS, which were slower to converge toward optimal parameters. Other widely used methods, such as the Powell's conjugate direction method, proved to be less effective for this type of problem.
The fitting of the model by optimization was influenced by the resolution of the grid and the time of growth simulated in the model. First, there was a minimum grid size, below which optimization algorithms did not converge to the true model parameter (Fig. 5A, B). The results also showed that longer experiments provide better parameter estimation. For the Nelder-Mead algorithm, when the experiment lasted 15 d, the true value of the gravitropic rate constant g of the primary root could only be obtained for grid elements (NY=12) along the depth (Fig. 5A), whereas for a 25 d experiment the same optimal value g was obtained at a lower grid resolution with NY=8 (Fig. 5A). Although good estimates for model parameters could be obtained using a small number of grid points, higher grid resolution is still required for accurate model predictions. For example, irrespective of the duration of the experiment, the optimal values of the branching rate b of the primary root and the elongation rate e of first-order laterals were obtained at NY=5 (Fig. 5E, F). However, the representation of root density distribution at this resolution is too coarse for practical application (e.g. studies of genotypic variation, nutrient uptake, or fertilizer applications).

Calculation time, spatial and temporal resolution
Calculation time increased with the number of elements in the finite element grid on which the model was run. The increase was not linear because the model was run on 3-dimensional grids (Fig. 6A) and the complexity of calculation increased as a power of 3. There was a linear relationship between the minimum time interval at which a change in growth parameter can be detected and the spatial resolution of the grid (Fig. 6B). The scale factor k was determined on different numbers of root systems. This was achieved by removing roots randomly from the database and running the cross-validation algorithm on the reduced data set. Simulation showed that increasing the number of roots reduced the scale factor k and therefore improved the spatial resolution (Fig. 6C, plain line). However, increasing the spatial resolution involves additional resources to acquire more measurements (Fig. 6C, horizontal isolines). Based on the results in Fig. 6B and C, it is possible to determine the optimal combination of time points and replicate them in an experiment (Fig. 6D). The number of time points as a function of the number of roots measured (replication) was calculated on the basis of a 100 d experiment. Increasing the number of roots measured at a given time point (replication) increases the spatial resolution of the model, and this also allows for improved temporal resolution. However, the increase was not linear so that increasing both spatial and temporal resolution usually requires a large increase in the number of replicates. In general terms, this result suggests the number of replicates contributes more to resolution than the number of time points in an experiment (Fig. 6D).

Estimation of root growth parameters using a density-based model
The model was fitted on data from B. rapa roots of 89 plants, 22 290 and 196 055 mid-points for primary and lateral roots, respectively. The primary root elongation rate was 1.19 cm d -1 and the gravitropic rate was 0.2 d -1 . For lateral roots, the branching rate was 3.52 cm d -1 , the root branching angle directly measured on image data was 7π/8, and the delay for lateral root initiation was 3 d. The optimal elongation rate and gravitropic rate of first-order laterals were 0.44 cm d -1 and 0.07 d -1 , respectively. The model-predicted total primary root length was 15.68 cm, while the estimated total root length was 15.54 cm. In addition, according to the model output, the total lateral root length per plant was 105.56 cm, whereas the real data-based estimate was 105.26 cm. Hence, the model slightly overestimated both the primary and lateral root length by 0.9% and 0.3%, respectively.

Accuracy of optimized root growth parameters for the density-based model
Root growth parameters obtained by optimization were compared with growth parameters estimated directly from the image data and the results showed that optimized parameters matched the experimental parameters. The experimentally estimated elongation rate of the primary root was 1.15 ± 0.05 cm d -1 (n=89), which agreed with the corresponding optimized value of 1.19 cm d -1 obtained by fitting the model parameters to the density data (see above). The experimentally estimated rate of increase in lateral root length was 0.45 ± 0.03 cm d -1 (n=89), which agreed with the value of 0.44 cm d -1 obtained by fitting the model to the data (see above). Sensitivity analysis also showed that predictions of the root length density varied differently for the same change in the elongation, branching, and gravitropic parameters. In response to a 1% change in parameter values, such as elongation rate and branching rate for the primary root and the elongation rate of laterals, the cost function rose by 0.947%, 0.996%, and 0.972%, respectively (MEV in Table 1). Conversely, the error was less sensitive to the gravitropic rate for either primary or lateral roots. A change in the model parameters by 1% induced a 0.073% and 0.003% change in the error function (Table 1). Bootstrap analysis of the standard errors showed that all growth parameters were estimated within the 95% CI of the model parameter. Elongation rates and branching rate were very sensitive parameters with CIs and CoV smaller than those obtained for the gravitropic rates (Table 1). Overall, our method leads to a good fit of the mathematical model to the experimental root length density data (Fig. 7A-D) with reasonable parity between predicted and measured root growth parameter values (Fig. 7E).

Using models to extract root growth parameters from experimental image data
Although the growth and development of individual roots is reasonably well understood, root system architectures (RSAs) observed experimentally often show complex patterns that cannot be easily characterized or predicted. RSAs are not merely self-similar fractal structures that simple reiteration of elongation and branching rules might predict. It is often observed, for example, that mortality and turnover are more frequent in thinner and higher order lateral roots (Norby and Jackson, 2000). Elongation and branching are also modulated by environmental and edaphic factors as well as systemic signals, so that more roots can be produced in regions of soils that are favourable to the plant such as localized nutrient patches (Hodge et al., 2009). Root growth can also be extremely variable, with, for example, 5-to 10-fold variations in elongation rates observed for lateral roots (Forde, 2009).
Obtaining suitable information to describe the growth of RSAs is therefore essential, but to date it has been a relatively major (low-throughput) undertaking (Tsegaye and Mullins, 1994;Tsegaye et al., 1995;Pellerin and Pagès, 1996). For these reasons, it has been proposed that optimization methods could be used to invert mathematical models so that growth parameters can be obtained from the models (de Dorlodot et al., 2007). In the past, simple empirical root depth models were used to fit root data extracted from soil cores. These models were rarely related to meaningful growth parameters that could explain the changes taking place in the root length density profiles through time. Reddy and Pachepsky (2001), Heinen et al. (2003), and Bonneu et al. (2012) have proposed the use of advection-diffusion equations to describe temporal changes in root spatial distribution in soil. In each case, the methods proposed could reliably mimic the spread of the root system in soil, but the models used in these studies were very simple and could not provide information on parameters such as the branching rate or the elongation rate of lateral roots.
More sophisticated architectural models have been proposed to explore whether a greater range of parameters can be extracted from mini-rhizotron data (Garré et al., 2012). However, when root architectural approaches are used, it is difficult to obtain the many architectural parameters to compare experimental with modelled root architectures. This is because root architectures, whether measured or modelled architecturally, are unique realizations of a stochastic process and there is no single obvious way to quantify their similarity. There have been various attempts to find ways of comparing RSAs: topological indicators such as magnitude and altitude (Fitter, 1987), fractal dimensions (Fitter and Stickland, 1992) or topological distances have all been proposed for direct comparisons between topologies (Ferraro and Godin, 2000), but they all fail to include a spatial representation of the root system. For these reasons, various forms of root density distribution have been used to compare root systems. For example, Grabarnik et al. (1998) used root length profiles (root length distribution with depth) to describe root distribution. Garré et al. (2012) used root number density (normalized number of roots) to compare mini-rhizotron observations at various depths with model predictions. Relationship between the spatial and temporal resolution is obtained by determining the minimum time interval over which growth parameters can be estimated at a given spatial resolution. Crosses indicate simulated results and the plain line indicates the linear fit to the data. (C) Spatial resolution increases with the number of replicates, but increasing the resolution required an exponential increase in the number of roots measured. Horizontal isolines indicate the number of grid elements required to simulate the growth of root systems at a given spatial resolution. (D) Optimal number of time points as a function of the number of roots measured (replication) calculated on the basis of a 100 d experiment. Isolines indicate the total number of measurements in an experiment. The intersection of an isoline with the curve indicates the optimal sampling strategy for a given number of measurements. In general terms, this result suggests that it is usually more important to increase replication than the number of time points in the experiment. Table 1. Sensitivity analysis of the estimated root growth parameters for primary root (e (0) ,g (0) ) and lateral roots (b (0) ,e (1) ,g (1) ) with the model elasticity value (MEV), the 95% confidence interval (CI) with sample size N = 100 and the coefficient of variation (CoV). Since some form of density estimation is required during an optimization process to compare a model with experiments, why not simulate the root system as a density distribution directly? This was the original idea of Reddy and Pachepsky (2001), Heinen et al. (2003), Bastian et al. 2008, andBonneu et al. (2012). However, because their density model was a generic advection-diffusion equation, it was not possible for them to obtain biologically meaningful growth parameters from their model. In the current study, the mathematical equations have been tailored to represent root systems so that growth parameters appear explicitly in the model. Because of this modification, it was possible to obtain accurate estimates of elongation rate, branching rate and gravitropic rate.

Accuracy of the estimated root growth parameters
This study has shown that by using density-based models, it is possible to extract root growth parameters such as Fig. 7. Measured and predicted root length density for the primary root (A, B). Measured and predicted root length density for lateral roots (C, D). The arrows illustrate mechanisms that may explain differences between experimental and predicted growth patterns. For the primary root, gravitropism induced non-smooth changes in root angles (A) that were not accounted for in the model (B). Lateral roots, on the other hand, seemed to have varying elongation rates through time (C) and these were not represented in the model (D). The color map indicates the root length density, and units of the color bar are cm −1 . (E) Comparison between measured root growth parameters (black) and estimated root growth parameters (grey). elongation rate, branching rate, and gravitropic rate directly from experimental density distribution functions. Although the data used to test the method were simple, they provided very useful information on growth parameters that can easily be obtained from experimental root density distributions. The study showed that elongation rates of primary or lateral roots were estimated with greater accuracy than the branching rate (Fig. 7E). Although direct experimental estimates for the gravitropic rate could not be obtained, the bootstrap analysis revealed that the gravitropic rate is the parameter that is estimated with the least confidence (Table 1), and it is likely that estimates of the gravitropic rate are less accurate than that of the branching rate. Results obtained in this study are difficult to compare with similar studies. Garré et al. (2012) could obtain the elongation rate and initiation angle using the RootTyp model, but there was much more experimental variation in their system compared with that used here because they grew barley in soil and had a limited number of observations given that they used mini-rhizotrons. Other studies by Bonneu et al. (2012) and Heinen et al. (2003) did not obtain estimates of root growth rates from their model inversion.
Most of the optimization algorithms tested converged toward the correct model growth parameters in the simulated data, but there were large differences in the convergence rate and in the time required to determine the optimal parameters. It is very interesting to note that the choice of the optimization algorithm is very specific to the problem studied. For example, Bonneu et al. (2012) used a combination of Nelder-Mead and simulated annealing algorithms to determine the parameters of their model. Reddy and Pachepsky (2001) used a Marquardt-Levenberg algorithm, and Heinen et al. (2003) used a Powell algorithm. Garré et al. (2012) applied the AMALGAM multicriteria optimization method.
The accuracy of the estimates of the growth parameters was dependent on the size of the finite element grid used to simulate the model (Fig. 5). Grids that were too coarse did not provide good estimates for the elongation rate because root tips were not traversing enough compartments of the grid. This made it more difficult for the optimization algorithm to detect the effect of changes in the elongation rate. The growth rate of primary roots was also more difficult to estimate than the growth rate of lateral roots (Fig. 7). This was because lateral roots were more numerous and occupied more soil volume than the primary root and it is, therefore, easier for the optimization algorithm to detect small changes in lateral root length density.
Not all the discrepancies in the estimation of the growth parameters were due to numerical approximations induced by the discretization of the solutions (Fig. 7). Another major limitation to accurate estimation of root growth parameters is the relative simplicity of the model compared with the complex behaviour of a real root system. The results presented here indicate an 18% difference between the branching rate of lateral roots predicted by the model and those measured experimentally (Fig. 7E). This difference could be explained by initial assumptions made in the model; for example, that the branching rate is constant. It is well known that the branching rate of laterals is very sensitive to environmental gradients and can vary during development (Van Norman et al., 2013). One way to overcome this limitation is to increase the complexity of the model with the addition of parameters that vary with time (Javaux et al., 2008). However, more complex models need more data to be parameterized accurately. For any additional model parameter, the size of data required to fit the model appropriately increases disproportionately, a phenomenon that was termed the 'curse of dimensionality' (Bellman, 1961). For example, if a time-varying coefficient is used for branching, then lateral roots must be measured at different times so that a response curve for the branching rate can be fitted. Choosing a suitable model is therefore equivalent to finding a good balance between model complexity and bias based on the data available, also termed the 'bias-variance trade-off' (Scott and Sain, 2005).

Data requirements for successful model inversion
The nature of the data available for optimizing the model will influence the quality of the estimations of model parameters directly. Systems such as aeroponics (Lobet et al., 2011), hydroponics (Mathieu et al., 2015), pouches (Adu et al., 2014) and thin rhizotron boxes (Rellán-Álvarez et al., 2015), in which a large fraction of the root system can be measured, require fewer replicates than experimental systems in which fewer roots are sampled, such as excavations (Bengough et al., 2000) or mini-rhizotrons (Garré et al., 2012), to determine growth parameters at the same level of accuracy (Fig. 6C).
It was also interesting to discover that the temporal resolution at which changes in model parameters can be measured is directly dependent on the spatial resolution of the grid and therefore on the scaling factor and the number of replicates taken at a given time point (Fig. 6B). One possible way to view this phenomenon is to consider the mathematical conditions under which it is possible to detect a change in the number of root tips in a cell of the finite volume grid. The time it takes for a root tip of velocity e to travel through a grid of size DX is exactly DX/e. Therefore, one may suggest that a criterion to be met for detecting a change in the occupation of space by a root tip is eDT/DX>K, with K being a constant to be determined theoretically or empirically. This relationship resembles the inverse of the Courant condition for the stability for finite volume simulation (Courant et al., 1967). However, in this case DT is not the time increment of the simulation, but the time increment between two experimental observations. Further, the constant is not unity and simulations indicate that K=4.8. Simulations also indicate that this condition is very accurate for describing the limits in temporal resolution (Fig. 6B). Based on these results, it is possible to provide a recommendation for sampling strategies to determine root growth parameters using this method of model inversion (Fig. 6D). The total number of root measurements (isolines in Fig. 6D) can be equated with the effort required for a given experiment. The same effort can be invested in more replicates and fewer time points or in more time points and fewer replicates. The optimal combination of replicates and time points is given at the intersection between the time point curve and the isoline. In general terms, it is usually more efficient to increase replication rather than increase the number of time points in an experiment.
Since more replication is required to improve the temporal resolution of estimates of growth parameters, the collection of root data can become a major limitation in root phenotyping analyses. Easing this limitation will require improving the spatial resolution to replicate ratio. One possibility would be to develop density estimation techniques optimized for RSA, for example by using parametric or semi-parametric techniques to estimate root density (Bishop, 2006). Alternatively, significant improvement could be obtained if the model could be fitted globally instead of being fitted by steps. The field of root modelling and phenotyping is evolving at a very fast pace, and there is no doubt that nascent technologies will overcome many of these limitations in the future.