Background and Assumptions

We review the prior work to the extent necessary to provide the context for our work. To overcome the obstacle of opaque soil, we grow rice and maize plants in gel with nutrients inside transparent cylindrical containers [13], monitoring the development of the root system for several weeks. To construct the 3D shape of a root system, we place the container on a rotating table connected to a camera and a computer, and we take 40 images at 9° increments. Currently, 20 of these images are used to reconstruct the shape of the root system represented as a collection of voxels in 3D space [24]. Letting t be the time of the acquisition, we write 𝕍_t for the reconstructed shape. Depending on the purpose, we will think of 𝕍_t as a subset of ℝ³, a collection of voxels, or a graph whose nodes are the voxels and whose arcs connect neighboring voxels. Drawn as straight line segments ending at voxel centers, the arcs have length 1, , or times the length of a voxel edge.

In order to analyze the growth of a root system, the same plant is reconstructed at different moments in time. We thus get a time-series of shapes 𝕍_t, for t = 0,1, …, T, which provides a much richer source of information about the root system than a single such shape. However, harvesting this information is more difficult since we need to compare, quantify similarities and differences, and resolve occasional contradictions. To do this automatically within our software platform, we make use of the following three assumptions:

• A₁: A root system grows in only two ways: by generating a new branch at a fork, and by elongating an existing branch at the tip.
• A₂: A root system is connected and contractible. In other words, the branches do not form loops, and they do not surround voids in space.
• A₃: The time-series is dense enough to observe the correct root hierarchy: when a side branch first appears at a fork, it is shorter than the parent branch.

While we have biological reasons to believe that A₁ and A₂ are true, both assumptions are sometimes violated by the reconstructed shapes. Our software can tolerate mild violations of all of the three assumptions but will fail when the reconstructed shapes contain gross contradictions to our expectations.

Depth, Seed, and Alignment

In this and the next two subsections, we explain how to construct the growth record of a root system, which provides a decomposition into a hierarchy of branches together with the time of growth. As a first step, we equip each reconstructed shape with a function φ_t;𝕍_t → ℝ that maps each voxel to its geodesic distance from a seed area, a concept we will explain shortly. We call this distance the depth of the voxel, and we call φ_t the depth function of 𝕍_t. After fixing a set of voxels S_t ⊆ 𝕍_t as the seed area, we compute φ_t using Dijkstra’s shortest path algorithm [25]. Running this algorithm on the graph of voxels is straightforward, but the result crucially depends on the choice of S_t, which is inspired by the biological idea that the root grows from a particular first location. It will be more important that the location of this area is consistent among all shapes in the time-series than to correctly identify the biological beginnings of the root system. We therefore compute S_T, and propagate its location to the other shapes after computing an initial alignment of all shapes in the time-series. We will explain the alignment first and return to computing the seed area thereafter.

Note that Assumption A₁ implies 𝕍_s ⊆ 𝕍_t for all 0 ≤ s ≤ t ≤ T. To align the two shapes, we can therefore use an algorithm that requires that the first shape be contained in the second. Our algorithm of choice is the 4-point congruent sets (4PCS) algorithm described in [22]. To apply it, we think of 𝕍_s as a collection of points in ℝ³, namely the centers of its voxels. The algorithm selects four approximately coplanar points from 𝕍_s and searches for a congruent configuration of four points in 𝕍_t. Let M be the matrix of the rigid motion that moves the 4-point configuration in 𝕍_t to that in 𝕍_s. While the 4PCS algorithm is fast, the alignment it produces can often be improved, and we use the iterative closest point (ICP) algorithm originally described in [26] for this purpose. It computes a second matrix M′, and we set M_st = M′ ⋅ M so that M_st(𝕍_t) is aligned with 𝕍_s. To align all T + 1 shapes in the time-series, it suffices to compute T rigid motions, and we choose to compute M_0t for t = 1,2, …, T. To align V_t with V_s, we then use . Suppose now that we have computed the seed area in the last shape, S_T ⊆ 𝕍_T. Using the rigid motions, we get S_t ⊆ 𝕍_t as the set of voxels closest to M_tT(S_T), for 0 ≤ t ≤ T.

To detect the seed area in 𝕍_T, we perform principal component analysis (PCA) on a small neighborhood of a radius R around every voxel. Covariance matrices with three similar eigenvalues indicate a similar spread of points along the principal axes, which in 3D space corresponds to a spherical neighborhood. We find the biggest connected component of voxels with roughly spherical neighborhoods and form the seed area by adding voxels with offset at most R from this component. For the results reported in this paper, we used R = 3 and we suggest that this parameter is set to a value from the interval [r,2r], in which r is the average radius of the branches. If the seed was not detected—for example when the set of voxels with roughly spherical neighborhoods is empty or larger than an allowed fraction of the total volume—we choose a voxel in the center of the topmost horizontal slice of the root system as the seed area.

Branch Hierarchy

As a second step in the construction of the growth record, we decompose each shape in the time-series into branches, which we organize hierarchically. This decomposition is inspired by the path decompositions of trees used in the design of efficient data structures; see e.g. [27]. Each branch starts at its tip at the bottom and extends upward, ending right before its fork, which lies somewhere between the two ends of the parent branch; we therefore represent each branch by the triplet of these entities: = (Tipⁱ, Forkⁱ, Parentⁱ). An exception is the topmost branch, which does not have a parent. Among all possible branch decompositions, we are interested in the one defined by the depth function.

To construct the decomposition, we traverse the voxels in the order of decreasing depth, adding an arc right after its second endpoint. With each visited voxel, we create a new component, and with each visited arc, we either form an additional connection within the component, or we merge two components into one. After visiting all voxels and arcs, they are part of a single component representing the entire root system, and we recover the branches and the relations between them by analyzing the history of events. To make this concrete, we equip each voxel, v, with a pointer to another voxel, ρ(v). In the beginning, each voxel points to itself. After completion of the algorithm, each voxel points to the tip of its branch, except for the tip itself, which points to the tip of the parent branch. A tip is therefore identified by having other voxels point to it, a branch is represented by its tip and consists in addition of all voxels pointing to this tip, and the parent branch is given by the pointer stored at the tip.

The only interesting event during the course of the algorithm is when it encounters an arc that connects two components. Let its endpoints be the voxels u and v, and write p = ρ(u) and q = ρ(v) for the corresponding tips, noting that p ≠ q by assumption. Assume without loss of generality that q is deeper than p: φ(p) < φ(q). If the difference in depth between u and p is smaller than some fixed threshold, then we declare the entire branch to be a part of the branch of q by updating the pointer: ρ(p) = q. Otherwise, we consolidate the branch of p by declaring v as its fork and the branch of q as its parent.

As mentioned earlier, each branch is represented by its tip, p, and consists of all voxels u with ρ(u) = p, together with p. To avoid possible confusions arising from the two different uses of the ρ-pointer, we save the pointer to the parent branch by setting Parent(p) = ρ(p), and thereafter set ρ(p) = p to indicate that it belongs to the branch it represents. In addition, we store a pointer Fork(p) to the fork in the parent branch. The thus linked branches form what we call a branch hierarchy. At this moment, the hierarchy is determined by depth: a branch has a tip that is necessarily less deep than the tip of the parent branch. Ultimately, we would like the hierarchy be determined by age. In the next section, we will explain how the hierarchy can be modified to achieve this goal.

Time Function

We are now ready to compute the time function of the last shape in the time-series, τ;𝕍 → ℝ with 𝕍 = 𝕍_T. Compressing the information of the time-series into a single shape, it maps every voxel of 𝕍 to the moment in time the voxel was created. The time function is piecewise constant, with possibly large steps in particular but not exclusively at forks. We note that the time function is the central piece of the growth record that enables the dynamic analysis of the root system.

To prepare the computation, we initialize the time function to τ(v) = T for every voxel v ∈ 𝕍. Traversing the other shapes in the time-series from back to front, we use the Euclidean distance between voxels to match the branches of 𝕍_t with those of 𝕍. It is therefore important that the shapes in the time-series are aligned and registered into a common coordinate system, which we henceforth assume. Specifically, for each tip in 𝕍_t, we find the closest voxel v ∈ 𝕍, and we match the branch of the tip with the branch 𝕁 in 𝕍 that contains v. Finally, we set τ(u) = t for all voxels u of 𝕁 whose depth satisfies φ_T(u) ≤ φ_T(v).

After running the algorithm, the set of voxels u ∈ 𝕍 with τ(u) ≤ t would ideally be the set of voxels in 𝕍_t, for every 0 ≤ t ≤ T. Of course, this will never be the case in practice, but there is a reason other than the inevitable noise in the system that may lead to discrepancies between the computed and the ideal time function, namely possible inconsistencies between the decompositions into branches between 𝕍_t and 𝕍. Instead of detecting inconsistencies between the branch decompositions, we repair the time function. Suppose v ∈ 𝕍 is a fork at which the classification into main branch and side branch is incorrect. By Assumption A₃, we get a conflict between the depth and the time. Specifically, we observe neighbors u of v for which (1) If we detect such a fork, we repair the inconsistency by an operation we call a switch. Letting p and q be the tips of the two branches at v, with v = Fork(p) but v ≠ w = Fork(q), we note that q = Parent(p). To perform the switch, we reassign the part of the branch of q between v and w to the branch of p; see Fig 1. In addition, the time at any voxel u within this part is set to the minimum of the current value and the value at v.

Root Traits

Using the growth record explained above, we compute a suite of traits that describe a root system captured in a time-series. Denoting the growth record by 𝕍 = 𝕍_T, we make essential use of the decomposition into branches, which we denote as . Recall that is the sublevel set of the time function on the branch. We call the voxel with maximum depth the tip of the branch at time t. We compute the following traits for every branch and at every moment of time 0 ≤ t ≤ T:

• volume: the number of voxels in ;
• depth of tip: the geodesic depth of the tip of ;
• location of tip: the three Cartesian coordinates of the tip;
• length: the difference between the geodesic depth values of the tip and the fork;
• switch event: a logical variable that indicates whether or not is longer than its parent branch;
• tortuosity: the ratio between the length of a branch and the Euclidean distance between tip and fork;
• average radius: the square root of the volume divided by π times the length;
• angle to gravity: the angle between the major PCA axis of and the vertical direction;
• angle to parent branch: the angle between the major PCA axes of and of its parent branch;
• number of children: the number of branches that have as parent branch.

There is an ambiguity about which part of a branch should be taken into account when computing the root emergence angle or direction of growth. In order to avoid this uncertainty, we use the major PCA axis of the branch which shows the main direction of the spread of all the voxels forming this branch.

Only two of the measurements relate the root system to the surrounding, namely the location of the tip and the angle to gravity. Both are useful to study the reaction of the root system to nutrients and other local environmental conditions. The measurements we output for each branch can of course be accumulated to obtain overall traits of the root system, such as the total volume, the total length, the total average radius, and so on. Similarly, we can take differences to compute dynamic traits, as for example the amount of growth from time t to time t + 1. Beyond the traits listed above, we compute the following characteristics of the seed area:

• seed volume: the number of voxels in the seed area of the root system;
• seed orientation: the angle between the major PCA axis of the seed area and the vertical direction.

Comparing the seed volume with the total volume, we get the volume fraction of the seed area. The concrete representation of the seed area is useful in turning the geometric hierarchy of the branches into a biologically more meaningful one. The former is defined by assigning the index 0 to the topmost branch and increasing the index by 1 whenever we pass from a branch to a child branch. A biologically more meaningful classification is into

• crown roots, which emerge above the seed area;
• seminal roots with forks located within the seed area;
• the primary root corresponding to the topmost branch;
• lateral roots, which are the children of other branches.

Unfortunately, crown roots often emerge from the stem above the reconstructed space, which results in disconnected branches. The reconstruction software used in this study [24] enforces connectivity and tries to find the shortest path between the biggest component of the root system and a disconnected branch. This shortest path does not always join the branch in a biologically correct manner, making the classification into biological classes more challenging. At the moment, we explore other ways to differentiate between crown, seminal, and lateral roots. In particular, elongation rates might provide a good clue for detecting crown roots as their initiation is developmentally delayed but typically proceeds faster than of the other root types, see Fig 2.

Datasets and Validation

We use the software to analyze four different time-series of 3D shapes, three obtained by imaging growing rice and maize plants, and one artificially created to validate the software. The latter consists of three models forming a caricature of a growing root system; see Fig 3. Created from a collection of discretized algebraic curves, the models are obtained by sweeping balls of varying radii along the curves to get branches with prescribed thickness. Observe that the time-series satisfies Assumptions A₁ to A₃ stated in Section 1. After designing the models, we 3D-printed them from resin, imaged the three shapes, reconstructed the 3D shapes from the 2D images, and finally processed the result with our software. Table 1 compares the hand measurements of branch length, depth of tip, volume, average radius, and tortuosity with the values computed by our software. The agreement in branch length and depth of tip is remarkably high, with a somewhat larger relative error only for very short branches. The tortuosity of a branch depends on the estimation of its length and agrees fairly well with our hand measurements, with only 6% average relative error. The volume and therefore the average radius are slightly underestimated by the software because the reconstruction of the 3D shapes tends to lose some of the thickness due to small misalignments of the 2D images. Nonetheless, the total volume estimated by our software differs on average only by 5.2% from the volume computed as a product of the mass of the printed model and the resin density. In order to validate the values of the angle to parent branch, we represent the each branch as a line segment passing through the tip and the fork and compute the angle between them. We note that the software computes the angle between the branches using PCA and it is impossible to measure this value by hand. Nevertheless, we observe a similar behavior in the values measured by hand and computed with the software: the angles decrease as branches elongate. Unfortunately, there is no a common way to measure branch angles, and DynamicRoots offers a new approach to estimate this trait.

The first of the three natural datasets consists of 478 time-series of 3 shapes each, reconstructing root systems of growing rice plants on Days 12, 14, 16 after planting. Using static traits, this data has been used in [16] to study the genetic control of shape and development of the root architecture. We intend to extend this study with the dynamic traits and the information about the branch hierarchy available from our software. The second natural dataset consists of 4 time-series of 25 shapes each, reconstructing root systems of growing rice plants imaged every two hours from 8am to 4pm for five days starting from Day 7 after planting. The third natural dataset consists of 3 time-series of 23 shapes each, reconstructing root systems of growing maize plants images from Day 1 to Day 9 after planting.

The relatively high temporal resolution in the last two datasets leads to detailed characteristics of the root development not available in the first dataset.

Traits

The traits described in Section 1 are output in a file for each time-series; see the supplementary material (S1 Table), which includes the file for the maize plant shown in the top row of Fig 2. For complete disclosure, our software also outputs files containing the color-coded geometry of the root system, representing the time function and the branch hierarchy, among other information; see Fig 2.

We use the remainder of this subsection to demonstrate how DynamicRoots can be used to analyze growth dynamics. As a first example, we look at the dynamic addition of volume to the shapes for the two displayed time-series. Each vertical bar in Fig 4 represents the root system at a moment in time, with the colored segments representing branches in the decomposition. The displayed behavior reveals a remarkable dynamics in the growth process: the volume of the maize root almost doubles in the course of just a few hours, with the bulk of the additional volume contributed by newly formed branches (see Fig 4 top, between 72 and 96 hours).

As a second example, we look at the branching dynamics, and in particular at the difference between the number of branches formed during day- and night-time. Fig 5 shows the number of branches as graphs over time for four rice plants. Observing the root systems during five consecutive days, we note that there are more branches formed during day-time than during night-time. The rice plants grew in the 12-hour cycle with the temperature of 28°C and 25°C at day and night. In [28] the authors studied diurnal elongation rates in Arabidopsis thaliana and observed higher elongation rates at night time than during the day. In [29] the authors demonstrated that growth conditions have a great impact on the diurnal growth patterns in sorghum and rice plants, they did not find differences in day-night elongation rates when the temperature was constant. In our data, we did not see differences in elongation rates at day and night time for plants in Fig 5. While more experiments are needed to understand the branching dynamics of root systems, it is clear that DynamicRoots can be of valuable assistance in collecting and analyzing the data.

Performance

We distinguish between performance characteristics that pertain to the speed and the correctness of the software. In addition to the traits discussed above, the software outputs performance statistics, including the running time, the average distance between aligned 3D shapes, and the number of switch events. For example, the average running time for one time-series from our first dataset, including file format conversion, graphical file output, shape alignment and processing, is 51 seconds on an Intel Core 3.4GHz desktop computer.

We use the remainder of this subsection to illustrate how we make use of the performance measurements. A delicate step in our algorithm is the alignment of the 3D shapes in a time-series. We underline the fact that the alignment performs quite accurately due to the property of the 4PCS algorithm which can robustly align shapes even with small overlap, noise, and outliers [22]. To shed light on the quality of the alignment, we compute the average distance between the aligned shapes and provide the result as an indicator of the quality of the achieved alignment. More precisely, we sample 10% of the voxels in the first shape, 𝕍₀, and we compute the Euclidean distance to the closest voxels in 𝕍₁ to 𝕍_T. For a high quality alignment, the average distance will be small, namely less than the edge length of a voxel. In the first dataset, there are several time-series for which the accuracy is compromised by technical difficulties with the imaging platform. For these series, the average distance we observe is significantly larger than for series that do not suffer from these imaging artefacts: e.g. 5.03 versus 0.84 times the length of an edge. In our analysis of the first dataset, we have used a threshold of 10.00 edge lengths for the average distance of the alignment beyond which a shape was excluded from the computations.

Another measurement that varies significantly with the quality of the alignment is the number of switch events. On average, we observe 9.04 switches between shapes in the first dataset. However, in the time-series affected by the mentioned image artefacts, this number increases to 58.96. The number of switches by itself can be an interesting biological trait as it measures how often a branch out-grows its parent branch. However, if the number is significantly higher than the average, then this may be reason for concern since it might be caused by technical difficulties in the acquisition pipeline rather than biological reality.