In Situ Root System Architecture Extraction from Magnetic Resonance Imaging for Water Uptake Modeling

Root system architecture is a crucial plant component of productivity: plant resistance to drought stress, nutrient acquisition, and plant yield are strongly influenced by the plant’s capacity to develop and adapt its RSA to environmental conditions. Although important progress has been achieved in understanding the molecular and genomic bases of RSA, the interaction of the RSA with the soil matrix is far less understood. This is mainly due to the opaque nature of the soil, which prevents direct observation of root systems and water uptake processes. Novel detailed models that rely on knowledge of the RSA have been developed for simulating nutrient and water uptake or root growth and performing in situ experiments for testing and validating the impact of specific root traits (Dunbabin et al., 2002; Doussan et al., 2006; Javaux et al., 2008; Draye et al., 2010). For these, not only the RSA is needed but also the full network of nodes with their connections. Characterizing the three-dimensional RSA in situ, however, is still a challenging task today. Various methods have been used over the years to obtain root architecture, starting with manual drawing (Kutschera, 1960), light transmission (Garrigues et al., 2006), rhizotron minicameras (Garré et al., 2011), and a three-dimensional imaging platform (Iyer-Pascuzzi et al., 2010; Clark et al., 2011). In the last two decades, however, noninvasive three-dimensional observation techniques with high resolution have been adapted for root–soil interaction research. These are x-ray tomography, which probes the physical density of the medium (Wildenschild et al., 2002), neutron tomography (Esser et al., 2010), x-ray microcomputed tomography (Mairhofer et al., 2012), and MRI (Pohlmeier et al., 2008, 2010; Rascher et al., 2011). Among the root measurement techniques, MRI possesses several great benefits. It can provide three-dimensional images of a strong heterogeneous sample, giving the exact location of the various structures within, with high spatial resolution. Knowledge about the relaxation properties of root tissue and the soil matrix (Pohlmeier et al., 2010, 2013) can be used to differentiate between roots and soil and obtain high-resolution three-dimensional architecture of roots growing in natural soil. Compared with two-dimensional scans or pictures, three-dimensional imaging may also provide more precise information on branching connections, explored soil volume, and branching angles. Ideally, all the connections between root segments should be characterized in a unique way (which is not the case with a two-dimensional image).

Limitations in using high-field MRI for root–soil imaging are mainly due to the presence of air bubbles and para- or ferromagnetic particles in natural soil, which lead to susceptibility artifacts caused by local magnetic field inhomogeneities. This shifts the resonance frequency of the root tissue, which is translated by the subsequent image reconstruction in a shifted location. The roots appear “pitted” or even gaps in the root strand can occur (Menzel et al., 2007). Another difficulty for achieving good contrast is in wet soils near saturation where transverse relaxation times of the soil matrix are relatively long compared with the relaxation in the root issue. In this case, long echo times must be used for fading out the signal from the soil. It should be noted also that information about the course of the root strands is sometimes interrupted and noisy in the other three-dimensional imaging techniques mentioned above, in common with MRI. Therefore, such data cannot be utilized as input for, e.g., inverse modeling or analysis of the RSA. What is needed is an image processing procedure that yields a three-dimensional network of root nodes that represents the RSA with no gaps, because such gaps prevent water flow. In addition, such a reconstruction should allow characterization of the root diameter and exact location distributions.

In this study, a novel automated method was developed and validated to recognize and reconstruct the three-dimensional RSA using a high-resolution (0.39-mm) MRI data set of a 3-wk-old lupin plant grown in sandy soil material (unpublished data, 2012). For validation, this automated method was compared with a semimanual three-dimensional reconstruction procedure in virtual reality and a two-dimensional analysis. The automatic reconstruction was performed using software developed at Bonn University based on a tube similarity measure for blood vessel detection in MRI images (Sato et al., 1998; Frangi et al., 1998). The semimanual reconstruction was performed using the virtual reality system Pi-Casso (Jülich Supercomputing Centre, Forschungszentrum Jülich) developed by S. Wienke and H. Zilken (Winke, 2010) based on the ViSTA software platform developed at RWTH Aachen University (Assenmacher and Kuhlen, 2010). The two-dimensional skeletonization is based on the SmartRoot algorithm developed by Lobet et al. (2011). The comparison between the different methods is based on statistical differences between topological indices. A transient water flow experiment using the Doussan et al. (1998) model was simulated using both three-dimensional root structures (semimanual and automatic) to assess the influence of the observed differences between the root structures on the water potential and water flux within the root xylem.

Our experiment consisted in monitoring the water content distribution dynamics in a pot with a lupin during a drought experiment. A lupin seedling was planted in an 8- by 10-cm Perspex cylinder filled with medium sand as a growing medium and watered from the bottom with Hoagland nutrition solution. Nine days after germination, the cylinder was saturated to a water content of 0.33 cm³ cm⁻³ and sealed at the bottom and top so that water loss was possible only by transpiration. During the following period, MRI images of the root system were taken after 0, 11, and 16 d to monitor root system development during increasing drought stress. From these series, the image obtained at the last date was chosen in this study for further analysis because in this state it was most similar to the photographs taken after finishing the experiments.

The MRI experiments were performed on a 4.7-T vertical-bore superconducting magnet (Magnex Scientific) operated by a Varian console. The nuclear magnetic resonance radio frequency resonator is a birdcage-type resonator with an internal diameter of 100 mm. The signal from the soil was efficiently suppressed by a long echo time of 20 ms and further by applying a threshold of 16% of the maximum signal. Repetition time was 30 s. The field of view was 100 by 100 mm at a matrix size of 256 by 256 pixels, resulting in an in-slice resolution of 0.39 mm and a slice thickness of 1.1 mm; 120 slices were monitored in interleaved mode. The total measuring time was about 4 h for two acquisitions (unpublished data, 2012).

For the automatic reconstruction, we used the algorithm for automatic root reconstruction described by Schulz et al. (2012). The algorithm starts with a three-dimensional MRI image. The unit of a three-dimensional image is a voxel (cf. pixel in two-dimensional images). Each voxel is described by a single gray value representing the signal intensity.

We first enhanced the raw data by emphasizing tubular structures (the roots) using an algorithm developed by Frangi et al. (1998) for blood vessel detection in medical MRI images. This results in another three-dimensional image, where values in voxels are larger if the local structure of voxels has similarity to a tube. The root model we desire should, however, be defined in terms of root segments and branches, not in terms of voxels. Our multistep procedure for the automatic extraction of the root model from the three-dimensional image is summarized here.

Tubular structures (root segments, in our case) can be found in the second-order derivatives of an image (Frangi et al., 1998). Locally, the image structure is similar to a tube if in one direction (where the root grows), the gray value does not change much, while in the other two directions (orthogonal to the first one), the gray value drops off equally fast and rapidly. In our case, the change in gray value is captured by the eigenvalues of a Hessian matrix, while the directions are captured by the eigenvectors. Whether a tube is found depends on the scale. Small roots are found when the Hessian is calculated in a small neighborhood, large roots are found when a larger window is considered. The neighborhood can be varied more efficiently by first convolving the image with a Gaussian distribution with standard deviation σ and then using a fixed neighborhood (e.g., all adjacent voxels) to determine the Hessian. The parameter σ is chosen proportional to the size of the roots that are intended to be found.

The MRI image is first preprocessed by up-sampling to isotropic (equal edge lengths) voxels using cubic spline interpolation. The resulting image is then convolved with a Gaussian distribution at scale σ, and the tube similarity measure by Frangi et al. (1998), a function of the eigenvalues, is computed. We repeat this process for all scales and find the maximum value for each voxel. In the resulting three-dimensional image, each voxel is represented by a number, which is larger if a scale exists at which the neighborhood of the voxel has the appearance of a tube.

The starting point for model extraction is the root collar, i.e., the position in the three-dimensional image where the plant exits the ground. We know the approximate slice (z coordinate) in the three-dimensional image where this happens due to the measuring process. To find the other coordinates robustly, we convolve the slice with a Gaussian with large σ to remove maxima, which are due to measurement noise. The maximum in the smoothed slice then has the x and y coordinates of the root collar.

Every voxel is now connected to the root collar position, but we already know that not all voxels are part of the root structure. The voxel-based tree needs to be pruned by removing branches that represent only part of the soil, leaving the tree representing the root system. For this purpose, we carefully select voxels that are definitely part of the root system. They and the whole path connecting them to the root collar are then retained for the final root system tree. The selection is based on two thresholds. The gray value of a selected voxel must be above some first threshold, υN, where N is the average noise level of the imaging procedure estimated in air above ground. If the signal/noise ratio is low and υ is small, this is not sufficient because all spurious measurements would be part of the root system tree.

We therefore define a second threshold, accepting only voxels for which the mass of their subtree (i.e., all their children, grandchildren, …, etc.) is larger than a threshold μN. The two thresholds can be chosen manually or automatically if the true root system is known (for details, see Schulz et al., 2012). The retained graph can now be represented as a graph data structure instead of an image. The retained voxels are the nodes of this graph, while the predecessor relationship defines its directed edges.

To determine the lengths, surface, and volume of a root, high-precision positioning of nodes is essential. So far, the nodes are positioned at voxel centers where we initialized them. Voxel centers are merely an artifact of the MRI procedure; the true root probably does not pass through voxel centers. We now apply a iterated mean-shift procedure to move the nodes to the center of the root with subvoxel precision. At each node n at position x_n, we estimate a covariance matrix C_n in a radius of 3 mm and determine its eigenvalues λ₁ ≤ λ₂ ≤ λ₃ as well as corresponding eigenvectors v₁, v₂, and v₃. If λ₃ > 1.5λ₂, we can assume that a clear direction is defined and v₃ corresponds to the local root direction. We then move the node to the mean of a neighborhood in the voxel grid, weighted by the tube similarity measure υ. To do so, we choose a four-neighborhood of x_n in the plane spanned by v₁ and v₂, and evaluate υ by linear interpolation.

Nodes where no main principal axis can be determined (λ₃ < 1.5λ₂) are moved to the mean of their immediate neighbors in the root graph. We iterate these steps until positions stop changing significantly.

We determine the radius of a root assuming a (truncated) symmetric Gaussian cross-section profile of roots. First, we fit a radial Gaussian function, M_n(r) = a exp(−b‖r‖²) + υN, to the MRI data in the plane defined by the second and third eigenvector of the covariance matrix C_n at every node n. The parameters a and b are estimated using unconstrained least squares optimization with the Levenberg–Marquardt algorithm starting at α = L(x_n), the interpolated MRI voxel intensity at the position of node n, and 1/5 of the expected maximum root radius r_max for the width, i.e., b = 12.5/r_max². Starting with a thin-root hypothesis biases the local optimization toward thin roots and avoids merging information of two neighboring root segments. The radius is set to b. Because estimates are noisy, especially for thin roots, we smooth the radius using a robust median filter across the closest six nodes in the graph.

Many parts in the described algorithm are inherently parallel. We therefore used a 12- by 2.67-GHz core Intel machine. The runtimes reported are for the 256 by 256 by 120 reference root, which was up-sampled to 256 by 256 by 256. Most processing time was spent for calculation of the tube similarity measure (30 s per scale). The shortest path selection took 3 min. Subvoxel positioning took 15 s, radius estimation 1.5 min. All other steps had runtimes totaling <1 min. The complete data set could therefore be processed in about 10 min. This was significantly shorter than the measuring time for one MRI image of the whole root system, which took about 1 h.

The automatically reconstructed RSA was compared with a semimanual method using a three-dimensional virtual reality (VR) system. The setup of the virtual reality system is shown in Fig. 1. It contained a virtual reality engine and input–output hardware to supply the platform for simulation. The engine reads the input devices, accesses databases, performs updates for real-time calculation of the virtual world and is responsible for presenting the results to the output device. These tasks were performed by a desktop computer in cooperation with the visualization system Pi-Casso, which includes a display, tracking cameras, tracked stereo glasses, a tracked gyro-stick (a three-dimensional mouse), and the space navigator as output–input devices. ViSTA is a software platform that allows integration of VR technology and interactive three-dimensional visualization (Assenmacher and Kuhlen, 2010). After cutting off voxels below a manually determined threshold, the MRI root image is visualized in VR.

The reconstruction of the root was performed manually by creating nodes and root segments using the gyro-stick, which is represented on the three-dimensional working environment by a traced pointer (Fig. 2). The user first defines the parent node, which is the point where the shoot enters the soil domain. This is referred to above as the “root of the tree” in the graph-theoretical sense. Further descendant nodes are then defined and automatically connected to the parent node by segments. The user must also define the thickness of the segments, which later defines their surface area through which water uptake takes place, by the gyro-stick. More nodes and segments are then added until the whole root system is defined. The developed software provides the user with the option to store and later reload the reconstructed root system.

The final RSA is hierarchically structured in a single tree and stored in ASCII format in the corresponding RSWMS input file (Javaux et al., 2008). In addition to the x, y, and z coordinates of each node, the identification number of the node and the branching order, the surface area of each root segment, and the length of the root segment and its mass are also calculated. The surface area of each segment is a morphological characteristic of particular importance because it affects the amount of water absorption. The time needed for the whole procedure was about 3 h. Detailed information about the development of the VR system can be found in Winke (2010).

The two-dimensional skeletonization of the lupin root was based on the SmartRoot image analysis toolbox developed by Lobet et al. (2011). First, the root was extracted from the soil, washed, and spread on a horizontal plate to avoid overlapping of branches as much as possible. Next, the root was photographed. The two-dimensional photo was used to trace the root architecture and topology using SmartRoot. According to Lobet et al. (2011), “SmartRoot is semi-automated image analysis software that combines vectorial representation of root objects with a tracing algorithm that determines the center of the root at a picked position by mouse click and, continues with a stepwise reconstruction of the root segments backward and forward to the tip and the basis of the root. It is implemented as a plug-in of the ImageJ software” (http://imagej.nih.gov/ij/). Detailed explanations on how SmartRoot works can be found in Lobet et al. (2011).

The obtained MRI image of the lupin root is presented in Fig. 3a. Discontinuities along root branches can easily be observed (see also Fig. 2). These large gaps are not realistic; they are the result of high threshold values imposed to differentiate the root structure from the surrounding soil. This effect demonstrates that reconstruction of the RSA is needed to obtain a continuous structure, which is necessary for the simulation and evaluation of hydraulic events. The manually reconstructed root architecture as obtained from the Pi-Casso virtual system is presented in Fig. 3b. The reconstruction allows root system skeletonization, i.e., transformation of an MRI image into a compatible file, which contains spatial coordinates of nodes, branch order, and the information about segment dimensions necessary for simulating root water uptake.

The root skeleton obtained from the Pi-Casso semimanual reconstruction has one ax (main root), 105 branches, 106 tips, and 3308 segments. In comparison, the root architecture obtained by automatic reconstruction is presented in Fig. 3c. The root skeleton obtained by automatic reconstruction has one ax (main root), 226 branches, 227 tips, and 2667 segments. When visually compared with the semimanual RSA, the automatic RSA shows small differences in the connectivity of the branches due to the fact that the automatic reconstruction was based on the Dijkstra algorithm, which searches for the minimum-cost path to connect two neighboring nodes, which is not necessarily optimal according to a biological model.

To quantify the differences between the two reconstructions, various indices were estimated based on the reconstructed RSA and are summarized in Table 1. It can be seen that the obtained RSAs vary in number of branches, segments, and connectivity of the segments but do not vary extensively in the average diameter of the segments (0.137 cm for semimanual and 0.143 cm for automatic), total root surface (127 cm² for semimanual and 124 cm² for automatic), total root length (293 cm for semimanual and 282 cm for automatic), and total root volume (4.7 cm³ for semimanual and 4.7 cm³ for automatic). In general, both RSAs exhibit comparable values for these parameters.

For a better assessment, these parameters were also compared with parameters obtained by tracing the same root from a two-dimensional photo using the SmartRoot image analysis toolbox developed by Lobet et al. (2011). The two-dimensional photo, together with its tracing, is shown in Fig. 3d. It shows overlapped branches, which cannot easily be traced individually, as well as thin roots close to higher order branches that have been considered as one thicker or thinner segment when the root diameter was estimated. The calculated parameters from the two-dimensional skeletonization are also included in Table 1. The total root lengths obtained by all three methods agree quite satisfactorily. This shows that MRI has detected the majority of all roots. The major difference is the smaller average segment diameter obtained from the two-dimensional (0.078 cm) compared with the reconstructed data (0.137 and 0.143 cm for the semimanual and automatic methods, respectively). This results in smaller root surface areas and smaller total root volumes. The reasons for this effect are not yet fully clear. One possible explanation is the partial volume effect of the MRI data, which classifies all voxels as roots irrespective of whether they are filled completely or only partially by a root.

Furthermore, we quantified the observed variations with the depth of the soil domain by describing the root architecture as density functions of both geometric and topological properties (Fig. 4) according to the mathematical approach of Dupuy et al. (2010). The soil domain was discretized and the density profiles of various root system parameters were calculated from the spatial location of the root segments and tips as functions of the grid element volume. The root length density is the total segment length divided by the grid volume and integrated over the axial slice. Analogously, the root surface and root volume densities are defined as integral segment and volume per slice. Root branch and apex densities aid in finding the distribution of the nodes at which branching occurs and the distribution of the apexes throughout the slice volume. The root length density, root surface density, and root volume density have nearly identical values and shapes of the profiles. Also, the shapes agreed very well with the root density profile obtained directly from the evaluation of the MRI images by axial integration (unpublished data, 2012). The distribution of the number of roots per depth of soil differed slightly between the semimanual and automatic structures, with fewer roots for the automatic RSA, especially in the upper part of the soil domain. This could be explained by different threshold values imposed in the semimanual and automatic reconstructions to suppress the signal from the soil. Due to differences in the data processing (a single threshold for semimanual RSA and two thresholds for automated RSA, see above), these threshold values can’t be set equal. Therefore, sparse pixels defining very fine roots with a signal intensity close to the threshold value may fade out together with the soil signal. The shape of the apex densities was similar for both root structures, with the difference that for the automatic RSA the profile was shifted to higher values. This can be explained by the difference in the number of traced root segments. The automated reconstruction generated more root branches, and thus roots apices, by artificially cutting root segments that should be connected. This is not a defining parameter of the roots, however, because in both methods the number of segments can be easily influenced and modified by the user. Branching density profiles show that for the automatic RSA the branching occurred in a more shallow soil layer than for semimanual RSA. This can be a measure of the differences in the connectivity of the branches, which varies for the automatic reconstruction due to the limitations of the Dijkstra algorithm.

To assess how these differences in the root architecture affect the root water uptake, simulations using both root structures were performed using the algorithm developed by Doussan et al. (1998). The results of the simulation are shown in Fig. 5. The water potential and radial flux of the semimanual automatic RSA are displayed, with visible differences both in values as well as in their distribution along the root network. The evolution of the water potential and water flux closely follow the transpiration rate applied at the root collar, but the response decreases in the direction of the root base for both structures as a result of uniformity assumed in the hydraulic conductivity distribution. The differences in the distribution of water potential and radial flux are a clear indication of the differences between the root structures because the most visible differences are localized in areas with differences in the connectivity of the branches. The most obvious difference occurs in the central, top region. In the semimanual RSA, the local fluxes are distributed across more root segments than in the automatic RSA. Hence, in this region, the local flow velocities are with about 0.7 cm d⁻¹, somewhat smaller than in the automatic RSA (about 1.2 cm d⁻¹).

In this study, we tested two methods of reconstructing the three-dimensional RSA of a lupin plant, grown in sandy soil, from high-resolution MRI images, with respect to root system connectivity and water uptake and transport properties. The first method was semimanual reconstruction using virtual reality and the second method was an automatic reconstruction based on local detection of tubular structures and globally enforced connectivity. The methods were validated by a two-dimensional classical skeletonization. The obtained root structures showed variations in the connectivity of segments. Fewer differences were encountered in the root parameters (number of branches and segments, diameters of the segments, total surface, and total volume). We assume this to be the influence of different threshold values imposed in the semimanual and automatic reconstruction of the roots to differentiate the root structure from the soil domain.

From the outcome of the two reconstructions, it was shown that automatic annotation gives objective results of equal quality to the semimanual one, suitable for large-scale experiments and repeatable for the same type of soil and root. Qualitatively, both methods shift the root diameter toward larger values in comparison with the MRI image resolution. To assess the influence of the observed differences between the root structures on root water uptake and xylem pressure head distribution, simulations using the Doussan numerical model were performed in a homogeneous soil and assuming a homogeneous distribution of root hydraulic conductivity.

The outcome of the simulations showed that the distribution of the water flux and xylem water potential is slightly different for each root structure, with lower values for the automatic reconstruction. The main differences in water flux and xylem water potential were localized in areas where the two root structures showed differences in segment connectivity. This was expected because in the automatic annotation algorithm, no experimental information about root segment connectivity is introduced. Incorporating such knowledge is the subject of further research.

We are grateful to Sandra Wienke (RWTH-Aachen) for developing Pi-Casso virtual system software during her master’s thesis; the CROP.SENSe.net research network (Project B4) for financial support, Dagmar van Dusschoten and Horst Hardelauf (FZ-Jülich, IBG2/IBG3) for technical support, Guillaume Lobet (UCL, Belgium) for the SmartRoot code, and Natalie Schroeder and Jan Vanderborght (FZ-Jülich) for helpful discussions.