Model-Based Design of Long-Distance Tracer Transport Experiments in Plants

Studies of long-distance transport of tracer isotopes in plants offer a high potential for functional phenotyping, but so far measurement time is a bottleneck because continuous time series of at least 1 h are required to obtain reliable estimates of transport properties. Hence, usual throughput values are between 0.5 and 1 samples h⁻¹. Here, we propose to increase sample throughput by introducing temporal gaps in the data acquisition of each plant sample and measuring multiple plants one after each other in a rotating scheme. In contrast to common time series analysis methods, mechanistic tracer transport models allow the analysis of interrupted time series. The uncertainties of the model parameter estimates are used as a measure of how much information was lost compared to complete time series. A case study was set up to systematically investigate different experimental schedules for different throughput scenarios ranging from 1 to 12 samples h⁻¹. Selected designs with only a small amount of data points were found to be sufficient for an adequate parameter estimation, implying that the presented approach enables a substantial increase of sample throughput. The presented general framework for automated generation and evaluation of experimental schedules allows the determination of a maximal sample throughput and the respective optimal measurement schedule depending on the required statistical reliability of data acquired by future experiments.

Introduction

Studying long-distance transport in plants is of high interest for the investigation of functional traits under the influence of diverse environmental factors (Van Bel, 2003; Jahnke et al., 2009). Non-invasive methods using short-lived radioisotopes have been established to detect the transport of radioactive tracer in vivo (Jahnke et al., 1998, 2009; Minchin and Thorpe, 2003; Alexoff et al., 2011; Garbout et al., 2012; Weisenberger et al., 2013; Hubeau and Steppe, 2015; Nakanishi, 2017). For example, after labeling a plant with ¹¹CO₂, transported ¹¹C can be detected and localized within the supplied plant organs from outside with positron emission tomography (PET) or scintillation detectors. These methods yield spatially and temporally resolved data of the tracer distribution within the plant, which can be analyzed by mathematical methods in order to estimate transport properties such as tracer transport velocities and leakage along the transport pathway (Tyree, 1975; Minchin and Thorpe, 2003; Suwa et al., 2008; Bühler et al., 2011, 2014, 2017; Converse et al., 2015). Further development of PET scanners dedicated for plant research (Streun et al., 2014, 2016) offers a high potential for plant phenotyping (Fiorani and Schurr, 2013; Hubeau and Steppe, 2015), but for this purpose the sample throughput needs to be increased. In view of the limited availability and high operating costs of radioisotope production and PET devices dedicated to plant research, a suitable way to increase throughput is by reducing the time of data acquisition per plant sample. Up to date, common tracer experiments in plants typically take at least 1 h of data acquisition plus sample handling times, limiting the throughput to not more than 1 sample per hour (Troughton et al., 1977; Jahnke et al., 2009; De Schepper et al., 2013; Converse et al., 2015). Increasing the throughput by simply shortening the data acquisition time will not work because information about transport properties as well as storage is typically spread over a large part of the time series (Bühler et al., 2011). Instead, the idea of the present study is to reach a higher sample throughput by introducing temporal gaps in the data acquisition of each sample, thus enabling to label multiple plants one after another with a temporal offset and to subsequently measure these plants in a rotating scheme. The resulting interrupted time series will always contain less information about the transport properties than a complete measurement. Thus the task was to identify experimental schedules regarding the distribution of gaps in data acquisition which maximize the throughput and at the same time minimize the loss of information. For this kind of evaluation, a model-based data analysis (Bühler et al., 2011, 2014) needs to be applied, because purely data driven methods (Minchin and Thorpe, 2003; Suwa et al., 2008; Converse et al., 2015) require uninterrupted time series. For the compartmental models established by Bühler et al. (2014), increasing gaps in the time series lead to higher uncertainties in the fitted model parameters. Thus the uncertainties can be used to rank experimental schedules according to the amount of lost information. To do this in a systematic way, a case study was set up, using a compartmental model and added noise to create data sets mimicking typical experimental data. These data sets were then used to evaluate a huge number of different experimental schedules. The working hypothesis was that the proposed scheme of interrupted time series created by rotational measurement should allow a substantial increase in sample throughput while maintaining a sufficient quality of data analysis.

Numerical Methods

Model Definition

The general class of compartmental tracer transport models was introduced in Bühler et al. (2014). This model class consists of up to N parallel compartments where all compartments are interconnected and allow exchange of tracer with individual exchange rates e_ij from compartment i to compartment j. Within the direction of spatial transport, tracer moves due to respective convection v_k and diffusion d_k in all compartments k, (k = 1… N). Accordingly, the underlying system of partial differential equations (PDEs) is defined as

$\begin{array}{rc}\frac{\partial {\rho }_k}{\partial t}=-v_k\frac{\partial {\rho }_k}{\partial x}+d_k\frac{{\partial }^2{\rho }_k}{\partial x^2}+{\displaystyle \sum _{j=1}^N{a}_{kj}}{\rho }_j& (1)\end{array}$

with ρ_k = ρ_k(x, t) as tracer density distribution in space and time. The exchange between compartments i and j with exchange rate e_ij as well as the decay of tracer with decay rate λ is combined in a_ij as

${\left(a_{ij}\right)}_{N\times N}=\left(\begin{array}{cccc}0& e_{12}& \cdots & e_{1N}\\ e_{21}& \ddots & & ⋮\\ ⋮& & \ddots & e_{(N-1)N}\\ e_{N1}& \cdots & e_{N(N-1)}& 0\end{array}\right)-\left(\begin{array}{ccc}{\displaystyle \sum _{n=1}^N{e}_{1n}}& & 0\\ & \ddots & \\ 0& & {\displaystyle \sum _{n=1}^N{e}_{Nn}}\end{array}\right)-\lambda I.$

Tracer entering the system is described either by a spatial initial condition or by a time dependent boundary condition (Bühler et al., 2014). From this general model class single model instances can be derived by allowing only a certain set of model parameters to be non-zero. In Bühler et al. (2014), a model filter was presented which sorts out all redundant models as well as useless models that would have non-functional compartments. For up to 5 model parameters (without diffusion) a family of 48 individual models was established, sorted by complexity and numbered consecutively from M01 to M48. In the present study, mainly model M13 (Figure 1) was used for numerical calculations since it has been shown to be universally applicable to typical experimental data (De Schepper et al., 2013; Bühler et al., 2014).

FIGURE 1

Figure 1. Sketch of model M13 from Bühler et al. (2014). There is transport in the first compartment with velocity v₁, exchange and back-exchange from the first to the second compartment with exchange rates e₁₂ and e₂₁, respectively, and storage in the third compartment with exchange rate e₂₃.

This model consists of three compartments with four model parameters that characterize the transport properties. Convection is controlled by v₁ and takes place in the first compartment only. Additionally, tracer can exchange (e₁₂) and back-exchange (e₂₁) from the first to the second compartment. Moreover, tracer can be stored in the third compartment (e₂₃).

Forward Simulation and Inverse Problem

The PDE in Equation (1) was spatially discretized using finite volume schemes as described in Bühler et al. (2017). Here, a linear fifth order upwind scheme was used (Shu, 2009). This discretization method produces quite accurate results when the initial condition is sufficiently smooth, and, at the same time, is very quick to solve (Bühler et al., 2017). The spatial discretization results in a system of ordinary differential equations (ODEs) which then was solved with a standard fifth order Runge-Kutta time solver (Dormand and Prince, 1980). Finally, the result of a forward simulation is the sum of tracer over all compartments,

$\rho \left(x,t\right)={\displaystyle \sum _{k=1}^N{\rho }_k\left(x,t\right).}$

The procedure of fitting the model to a given data set was described in Bühler et al. (2014). A nonlinear least squares optimization method was used for estimation of model parameters (lsqnonlin, Matlab R2016a, The MathWorks, Inc.).

Parameter Uncertainties

The model fitting yields estimates of the model parameters and the Jacobian matrix of the parameters at the optimum. Asymptotic standard errors of the parameter estimates were estimated from the square roots of the diagonal elements of the covariance matrix Cov of the model parameters. Cov was approximated by Cov = s²(J^T J)−1, using the variance s² of the fit and the Jacobian J of the estimated model parameters at the minimum (Johnson and Faunt, 1992). Though this method to estimate standard errors neglects the covariances between the estimated parameters, it is preferable here to the Monte Carlo bootstrap method used in Bühler et al. (2014), because it is computationally much less expensive and also applicable to the very small data sets considered in this study. Example comparisons showed that the results of both methods are very similar for large data sets used in this study.

Implementation

All numerical routines have been implemented in MATLAB (R2016a, MathWorks, Inc.). The code is publicly available at https://github.com/ForschungszentrumJuelich/mdate.

Experimental Design

When measuring multiple plants, each individual sample can be labeled with a temporal offset to its predecessor, and the plants are subsequently measured in a rotating scheme (Figure 2). The temporal distance (gap size d) between successive measurements of the same plant depends on the total number of plants measured, the length of the measurement window w, as well as the handling time T_h needed for transporting samples in and out of the measurement setup.

FIGURE 2

Figure 2. Consecutive measurement of multiple plant samples. For the pulse labeling experiment of plant 1 the data acquisition begins at time point T_start and is interrupted regularly to allow interlaced measuring of multiple samples.

An experimental design is then defined by the following experimental parameters: (1) starting time point T_start, (2) window width w, (3) number of windows N_w, (4) handling time T_h, and (4) gap length between windows d. The gap d between two respective measurement windows is defined as a multiple of w + T_h, ensuring that all further samples can be measured subsequently using the same experimental design. This procedure allows investigating the experimental design for one sample only which then represents all consecutive sample measurements. The construction of the experimental designs is depicted in Figure 3.

FIGURE 3

Figure 3. Schematic construction of measuring scenarios. Data are sampled with windows varying in starting time, frequency, length and distance between windows. The length of the gaps d between the windows are multiples of window width w plus handling time T_h. In this way, all designs allow an interlaced measurement of multiple samples and an optimized use of measurement equipment.

One of the main properties of each experimental design is the mean sample rate SR which represents the number of samples that can be measured per time unit. This sample rate is given by SR = (N_w · (w + T_h))⁻¹. The cumulated time of sample measurement can be calculated as N_w · w for each design.

Case Study

Reference Data

For this study a general reference data set was created with model M13 and artificial noise. A Gaussian distribution ${\rho }_0(x)=exp(-{(x-x_0)}^2/2{\sigma }^2)$ served as a continuous initial tracer distribution function at time t₀ = 0 with x₀ as distance from the beginning of the experimental field of view at x = 0 and σ as width of the input curve. One of the major advantages of the spatial initial condition is the ability to deal with temporal gaps in the analyzed data and with tracer leaving the field of view, both cases in which a reconstruction of a temporal boundary condition is not possible. Also, implementing a spatial initial condition allows high convergence orders of the PDE solver at the boundaries of the field of view (Bühler et al., 2017). The parameters x₀ and σ as well as the amplitude of the simulated tracer distribution within the field of view need to be fitted along with the model parameters. For fitting the model to experimental data, also the temporal position of the initial condition, t₀, might need to be estimated, because a temporal shift of experimental data could affect the fit quality. Nevertheless, estimating t₀ requires a continuous interpolation of the data to the shifted time grid. This interpolation can cause small numerical oscillations of the non-linear optimization near the minimum which can affect the comparison of different experimental designs. For this reason, t₀ was set to zero and not fitted in this study.

For the calculation of the reference data set, the parameter values of model M13 were set to v = 2.0 mm min⁻¹, e₁₂ = 0.3 min⁻¹, e₂₁ = 0.1 min⁻¹, e₂₃ = 0.05 min⁻¹, and σ = 5 mm, x₀ = 60 mm for the initial condition. The choice of model parameters was based on the fit parameters of model M13 to PET data of sugar beet in Bühler et al. (2014). The forward simulation was evaluated at 11 spatial positions from x = 0 to 100 mm in steps of 10 mm and with a temporal resolution of 1 min for a total time of 3 h. In order to simulate a certain model error, this data was repeatedly manipulated by adding normally distributed noise with a standard deviation of 7e–3. This value was chosen because the resulting standard errors of the estimated parameters are in the same range as the values for sugar beet in Bühler et al. (2014). The reference data set with one of the noise patterns used for all further calculations is shown in Figure 4.

FIGURE 4

Figure 4. Example reference data set based on forward simulation of M13 with parameters ν = 2.0 mm min⁻¹, e₁₂ = 0.3 min⁻¹, e₂₁ = 0.1 min⁻¹, e₂₃ = 0.05 min⁻¹, σ = 5 mm, x₀ = 60 mm and added normally distributed noise with a standard deviation of 7e-3.

Specification of Experimental Designs

The experimental designs described in section Experimental Design were applied on the reference data, resulting in reduced data sets, which are referred to as “designs” in the following. In order to make the number of possible designs feasible for calculations, the following side conditions were applied. Window widths w were considered in a range of values from 1 to 120 min in steps of 1 min. The starting time point T_start was varied in steps of 1 min from 20 to 40 min, i.e., starting data acquisition at the time of tracer arrival at the field of view and later. These conditions limited the maximal number of windows N_w to 26 which were iterated in steps of 1. The gap length d was then increased in multiples of the sum of window width w and handling time T_h, until a maximum of around 120 min after starting the experiment was covered. For all calculations the handling time T_h was set to a constant value of 1 min. This procedure resulted in 12,547 different designs. Figure 5 shows the resulting sample rates SR for a subset of N_w = 1…10 and w = 1…20. For all calculations in the following, only sample rates ≥1.0 h⁻¹ will be discussed. Areas with lower sample rates were not considered and were therefore labeled in white in Figure 5.

FIGURE 5

Figure 5. Exemplary view of the sample rate SR which is the number of samples measurable per hour for different experimental designs with varying number of windows and window widths.

Results and Discussion

All 12,547 designs were fitted by model M13, taking as starting points the model parameters which were used for generating the reference data. The fitting results allowed the calculation of variances of the model parameter estimates for each design. Subsequently, the designs were filtered with regard to specific quality criteria, resulting in a set of possibly best designs for a given throughput. The filtered designs were further analyzed for sensitivity with regard to experimental variations in arrival time of tracer in the field of view. Finally, a general procedure for the generation of optimal experimental schedules was deduced.

Error Estimation of Designs

The noise in the data caused the parameter estimates always to differ from the original underlying parameter set. The fitting was repeated five times for each design with different random noise and the fitting results were averaged. The definition of all designs and the respective fitting results can be found in Supplementary Table 2. The sum of the relative standard errors, SE_sum, of the four most relevant model parameters (velocity v and exchange rates e₁₂, e₂₁, and e₂₃) was used as a measure of the overall uncertainty of each design. This measure has the advantage to be independent of the different absolute parameter values. The standard errors are expected to roughly correlate with the cumulated time of measurement i.e., the longer the time of data acquisition, the lower the standard errors and, hence, the lower SE_sum. Standard errors for the parameters x₀ and σ of the initial function were not taken into account because the initial function is of purely auxiliary nature (Bühler et al., 2011). Figure 6 shows this relationship in detail for all designs up to a maximal SE_sum of 80%. Designs with a low number of windows N_w ≤ 5 are highlighted respectively by colored rings.

FIGURE 6

Figure 6. Uncertainty of model parameters SE_sum vs. cumulated time of measurement for all design with SE_sum ≤ 80%. Designs with a small number of windows are highlighted by colored rings, respectively.

In Figure 6 some groups of continuous relations become apparent especially for designs with 1 and 2 windows. Nevertheless, of special interest are only designs that are located on the Pareto curve of a minimal SE_sum for each unique value of cumulated time of measurement. Additionally, all respective Pareto elements with a limited number of windows N_w ≤ 5 were considered, too. These minimal uncertainty designs were filtered out of the set of 12,547 designs and plotted vs. the respective sample rate SR in Figure 7. Here, designs with SE_sum higher than 70% are not regarded and, again, designs with a low number of windows N_w ≤ 5 are highlighted by colored rings. The maximal value of 12 samples per hour relates to a total measurement time of 5 min per sample.

FIGURE 7

Figure 7. Plot of all potentially best designs from Figure 6 after filtering by the minimal sum of parameter uncertainties SE_sum for each unique value of cumulated measurement time. Additionally designs were included with a low number of windows between 1 and 5. These are highlighted by colored rings, respectively.

For each number of windows, approximately linear relations of model uncertainties SE_sum to sample throughput are recognizable in Figure 7 for higher values of SR. Deviations from a clear linear trend could result from random patterns in data. The mean slope of each relation decreases with increasing number of windows up to N_w = 5. Designs with N_W > 5 do not significantly improve the relation between model uncertainties SE_sum and throughput SR any further. Designs with 2 windows cannot compete with any other designs. For low throughput with values < 1.2 samples h⁻¹, designs with 1 window always prove to be an optimal choice. Above this value the selection of an optimal design becomes complex and designs with N_w ≥ 3 will be preferable.

Selected Designs

The following quality criteria for designs can be applied in order to filter optimal designs from the set of all possible designs. Any potentially best design

(1) must provide a minimal value of SE_sum for its respective values of throughput SR and number of windows N_w;

(2) is preferred over another design if it has a lower number of windows N_w while throughput SR and uncertainty SE_sum remain similar;

(3) allows a unique fitting to data for a preferably large range of starting parameter sets;

(4) enables the identification of the underlying model;

(5) shows a certain stability in SE_sum regarding shifts in T_start relative to the arrival time of tracer in the field of view.

Per definition, all pre-filtered designs shown in Figure 7 already fulfill quality criterion (1). By applying also quality criteria (2) to (4) to the pre-filtered designs, a set of 12 designs was selected covering sample rates from 1 to 12 samples h⁻¹. The design definitions and corresponding values of SR and SE_sum are listed in Table 1. SE_sum correlates almost linearly with the noise amplitude (see Supplementary Figure 1) and, thus, the noise amplitude has no significant influence on the comparison of designs regarding parameter uncertainty. Validity of criterion (3) was confirmed by repeated successful fitting, starting with different parameter values far away from the minimum. Criterion (4) was ensured by comparing the results for fitting the designs with different selected model instances from the general model class, see Supplementary Material S5. SR ≥ 1 was assumed in order to ensure a significant improvement of throughput compared to the full reference data set. Designs with even higher sample rates than SR = 12 h⁻¹ were not considered, because their corresponding high model uncertainties SE_sum of more than 50% severely limit their practical applicability.

TABLE 1

Table 1. Resulting selected designs from application of quality criteria (1) to (4).

For the designs in Table 1, the relative standard errors of the single model parameters adding up to SE_sum increase uniformly with SE_sum (see Supplementary Table 2), but there are deviations from the proposed linear increase of SE_sum with increasing sample rate SR for designs 7 vs. 6 and 9 vs. 8. This results from applying criterion (2), by which designs for SR = 3.75 and SR = 5 with slightly lower SE_sum, but much higher N_w were dismissed in favor of designs with N_w ≤ 5. Design 11 is the only exception from the approximate rule of N_w ≤ 5 for the selection in Table 1, included here because there was no other design matching the criteria in the range of 8 ≤ SR < 12 (see Figure 7 and Supplementary Table 2). The rule of N_w ≤ 5 was motivated by the fact that the specific handling time T_h = 1 min used in the calculations is a very rough estimate only. Handling time could be in the range of seconds for a, yet to build, efficient automated sample handling system and up to several minutes for manual sample handling. Also it is generally desirable to minimize disturbances of the samples by movements as well as to minimize manual handling of radioactively labeled samples by the experimentalist. Together, these arguments lead to the restriction of N_w applied in this study.

The application of the selected designs to one of the reference data sets is visualized in Figure 8, illustrating the temporal distribution of the respective measurement windows. All these possibly best designs show certain common characteristics. For example, all designs start data acquisition between 23 and 25 min. This suggests that the rise of the leading curves in the reference data, corresponding to the arrival time of data in the field of view, contains important information which should always be covered. Also, almost all designs cover at least one of the peaks as well as the steepest areas of descent of one the dominant curves. The tail of the tracer profiles does not seem to contain essential information, since none of the possibly best designs includes data points after 93 min. Given that the time point T_start, which marks the beginning of data acquisition, seems to be of crucial importance for the performance of all designs, criterion (5) will be investigated more thoroughly in the following.

FIGURE 8

Figure 8. Visualization of exemplary reference data of selected experimental designs from Table 1. SR, sample rate; CTM, cumulated time of measurement.

Optimal Starting Time Points

For experimental setups with series of plants there will always be variations in the time point of tracer arrival in the field of view which are caused by structural biological differences between investigated samples. Therefore, experimental designs are preferable if they are less sensitive to variations in the starting time of data acquisition which relates to the time of tracer arrival in the field of view [criterion (5)]. In order to analyze these dependencies more closely, the model uncertainties SE_sum of all designs from Table 1 were plotted in Figure 9 with starting time points T_start varying from 20 to 40 min in steps of 1 min.

FIGURE 9

Figure 9. Plot of parameter uncertainty SE_sum for varying starting time points T_start starting from the 12 possibly best designs of Table 1. The horizontal dotted line shows mean value of SE_sum, the vertical dotted line indicates time point of minimal SE_sum. SR, sample rate; CTM, cumulated time of measurement.

These plots reveal that designs with increasing throughput not only show higher minimal values of SE_sum, but also consistently higher averages of SE_sum over the considered period of starting points and increasing sensitivity of SE_sum to variations of T_start. For designs 9–12, the high sensitivity together with the lack of a plateau around the minimum is putting the practical applicability in doubt. In general, the time point of minimal SE_sum is not necessarily the best starting time point with regard to the sensitivity. For example, for designs 6–11, starting the data acquisition just 2–3 min earlier would cause significantly larger model uncertainties than starting some minutes after the minimum. This is caused by the fact that starting time points smaller than 20 min lead to data acquisition before tracer enters the field of view. Thus, these data windows contain less or no information at all about the tracer distribution, which causes an immediate increase of parameter uncertainties. In these cases, starting data acquisition some minutes later than at the minimum value of SE_sum would be a good strategy to avoid larger model uncertainties caused by biological variations in arrival time. In order to identify optimal starting points, pre-experiments with at least one non-interrupted measurement have to be performed in a specific experimental situation.

General Procedure for the Identification of Optimal Designs

Pre-experiments will also be needed to determine the level of accuracy needed to enable identification of significant differences in an experimental study. A general procedure to arrive at optimal schedules will then consist of the following steps: (i) determination of a maximal uncertainty SE_sum of the parameter estimates depending on the required accuracy level in the experimental study; (ii) analysis of possible designs based on non-interrupted test-measurements, in the same way as presented in the preceding sections; (iii) selection of an optimal design with a maximal number of plants measureable per hour, depending on the maximal uncertainty SE_sum. There will be no need to test all of the possible designs investigated in the present fundamental study. An adaption of the presented candidate designs to the actual shape of the tracer data is expected to be sufficient. As a general guideline, only such designs qualify as candidates for an optimal experimental design which cover the time point of tracer arrival in the field of view as well as the maxima and declining sections of the dominant first curves of the data set. Also, designs with 1 and 2 windows do not need to be considered for higher throughput scenarios. These constraints substantially reduce the number of candidates and enable a faster determination of optimal design for any given experimental setup.

For the case study, modeled instead of real data were used, which were known to represent real PET data of tracer transport in sugar beets very well (Bühler et al., 2014). The use of modeled data avoided obscuring the results by an (albeit tiny) model error, thus facilitating the systematic analysis. Since the presented case study represents only one exemplary experimental situation, two additional case studies with different transport properties were investigated as well. The results are available in Supplementary Materials S3, S4. In the case study based on tracer transport in maize root (Figure S3.1, parameter values from Bühler et al., 2014), the analysis was based on the simpler model M02 of the model class. This model consists of a transport—and a storage compartment only, without back-exchange of tracer. In the case study representing tracer transport in oak stems (Figure S4.1, parameter values based on De Schepper et al., 2013), model M13 was used, but with a lower spatial resolution, i.e., only 3 spatial positions instead of 11. As expected for different experimental setups, the level of SE_sum and the specific best designs for each sample rate differed for the different case studies. But qualitatively the results were very similar: calculation of SE_sum for all experimental designs (Figures S3.2, S4.2, in analogy to Figure 6) allowed the identification of possible best designs (Figures S3.3, S4.3, Tables S3.4, S4.4, in analogy to Figure 7 and Table 1), which could be further analyzed regarding variations of the starting time point T_start (Figures S3.5, S4.5, in analogy to Figure 9). Thus the general procedure works in all three cases, indicating that it is largely independent of the experimental setup as represented by the tracer transport properties and the spatial resolution of the measurements.

Outlook

The quantitative results presented here depend on the specific properties of the chosen case study. Yet the workflow of design generation and evaluation can be applied analogously for other cases, e.g., with data from different plant species, or other spatial and temporal fields of view (possibly from other modalities of data acquisition than PET). The results could also be transferred to experimental setups using other isotopes than ¹¹C, e.g., ¹⁸F with a half-life of 109.8 min (Partelová et al., 2014; Converse et al., 2015). Preliminary investigations, simulating the transport of ¹⁸F, showed a similar selection of possibly best designs, with a tendency of later optimal starting time points due to the longer half-life of this isotope. The construction of designs described in section Experimental Design could be extended by considering also non-regular patterns with different temporal distances between measurement windows. As for regular patterns, there should be no gaps between interlaced measurements of multiple samples in order to make full use of the measurement equipment. Figure 10 shows examples of symmetrical as well as non-symmetrical non-regular patterns meeting this requirement. From these basic patterns a set of 336 different designs was derived (Supplementary Table 2) and analyzed analogously to the regular pattern above. Though none of these designs showed any improvement over the corresponding regular designs with same throughput, this still might be the case for other non-regular patterns or experimental scenarios not considered so far.

FIGURE 10

Figure 10. Examples of non-regular consecutive measuring patterns of multiple plant samples, denoted by numbers, in analogy to Figure 2. Scenarios (A,B) are symmetric, (C,D) are non-symmetric. The process of window splitting used to create (B) out of (A) can also be applied for (C,D).

Conclusions

In this study the experimental design of interlaced labeling experiments of multiple samples for phenotyping internal transport mechanisms was investigated. A representative case study with a multitude of different designs was set up and evaluated with regard to accuracy of model-based analysis. For a range of different sample throughput scenarios from 1 to 12 samples h⁻¹, there was always an optimal design. It is remarkable that in high throughput scenarios the model parameters could still be reconstructed with quite a small number of data points.

The presented procedure of design construction and evaluation serves as a framework for developing optimal experimental schedules by increasing sample throughput and at the same time keeping a required statistical reliability or, the other way around, minimizing the statistical uncertainty for a required throughput. For practical applications with a specific experimental setup, the maximal number of plants measureable per hour will not only depend on the required accuracy, but also on sample handling time. Before starting a series of interlaced plant experiments it will be necessary to perform thorough preliminary investigations on the plant species or genotypes of interest and test the presented procedure of design evaluation on complete data sets in order to find respective best case designs.

Author Contributions

JB, EvL, and GH designed the study; JB and GH performed the research; all authors contributed to writing the manuscript and approved the submitted version.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Siegfried Jahnke for useful discussions of PET issues and critical reading of the manuscript. JB wants to thank Jakob Vanhoefer for helpful discussions regarding data analysis.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00773/full#supplementary-material

Supplementary Figure 1. Dependence of parameter uncertainty SE_sum on noise. Noise level is defined as the standard deviation of added normally distributed noise.

Supplementary Table 2. Data table including all design definitions and fit results.

Supplementary Material S3. Results of additional case study based on maize root transport properties.

Supplementary Material S4. Results of additional case study based on oak stem transport properties.

Supplementary Material S5. Model identifiability.

References

Alexoff, D. L., Dewey, S. L., Vaska, P., Krishnamoorthy, S., Ferrieri, R., Schueller, M., et al. (2011). PET imaging of thin objects: measuring the effects of positron range and partial-volume averaging in the leaf of Nicotiana tabacum. Nucl. Med. Biol. 38, 191–200. doi: 10.1016/j.nucmedbio.2010.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bühler, J., Huber, G., Schmid, F., and Blümler, P. (2011). Analytical model for long-distance tracer-transport in plants. J. Theor. Biol. 270, 70–79. doi: 10.1016/j.jtbi.2010.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bühler, J., Huber, G., and von Lieres, E. (2017). Finite volume schemes for the numerical simulation of tracer transport in plants. Math. Biosci. 288, 14–20. doi: 10.1016/j.mbs.2017.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bühler, J., von Lieres, E., and Huber, G. (2014). A class of compartmental models for long-distance tracer transport in plants. J. Theor. Biol. 341, 131–142. doi: 10.1016/j.jtbi.2013.09.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Converse, A. K., Ahlers, E. O., Bryan, T. W., Hetue, J. D., Lake, K. A., Ellison, P. A., et al. (2015). Mathematical modeling of positron emission tomography (PET) data to assess radiofluoride transport in living plants following petiolar administration. Plant Methods 11:18. doi: 10.1186/s13007-015-0061-y

PubMed Abstract | CrossRef Full Text | Google Scholar

De Schepper, V., Bühler, J., Thorpe, M., Roeb, G., Huber, G., van Dusschoten, D., et al. (2013). 11C-PET imaging reveals transport dynamics and sectorial plasticity of oak phloem after girdling. Front. Plant Sci. 4:200. doi: 10.3389/fpls.2013.00200

CrossRef Full Text | Google Scholar

Dormand, J. R., and Prince, P. J. (1980). A family of embedded Runge-Kutta formulae. J. Comput. Appl. Math. 6, 19–26. doi: 10.1016/0771-050X(80)90013-3

CrossRef Full Text | Google Scholar

Fiorani, F., and Schurr, U. (2013). Future Scenarios for Plant Phenotyping. Ann. Rev. Plant Biol. 64, 267–291. doi: 10.1146/annurev-arplant-050312-120137

PubMed Abstract | CrossRef Full Text | Google Scholar

Garbout, A., Munkholm, L. J., Hansen, S. B., Petersen, B. M., Munk, O. L., and Pajor, R. (2012). The use of PET/CT scanning technique for 3D visualization and quantification of real-time soil/Plant interactions. Plant Soil 352, 113–127. doi: 10.1007/s11104-011-0983-8

CrossRef Full Text | Google Scholar

Hubeau, M., and Steppe, K. (2015). Plant-PET scans: in vivo mapping of xylem and phloem functioning. Trends Plant Sci. 20, 676–685. doi: 10.1016/j.tplants.2015.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Jahnke, S., Menzel, M. I., Van Dusschoten, D., Roeb, G. W., Bühler, J., Minwuyelet, S., et al. (2009). Combined MRI-PET dissects dynamic changes in plant structures and functions. Plant J. 59, 634–644. doi: 10.1111/j.1365-313X.2009.03888.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jahnke, S., Schlesinger, U., Feige, G. B., and Knust, E. J. (1998). Transport of photoassimilates in young trees of Fraxinus and Sorbus: measurement of translocation in vivo. Bot. Acta 111, 307–315. doi: 10.1111/j.1438-8677.1998.tb00714.x

CrossRef Full Text | Google Scholar

Johnson, M. L., and Faunt, L. M. (1992). Parameter estimation by least-squares methods. Meth. Enzymol. 210, 1–37. doi: 10.1016/0076-6879(92)10003-V

PubMed Abstract | CrossRef Full Text | Google Scholar

Minchin, P. E. H., and Thorpe, M. R. (2003). Using the short-lived isotope 11C in mechanistic studies of photosynthate transport. Funct. Plant Biol. 30, 831–841. doi: 10.1071/FP03008

CrossRef Full Text | Google Scholar

Nakanishi, T. M. (2017). Research with radiation and radioisotopes to better understand plant physiology and agricultural consequences of radioactive contamination from the Fukushima Daiichi nuclear accident. J. Radioanal. Nucl. Chem. 311, 947–971. doi: 10.1007/s10967-016-5148-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Partelová, D., Uhrovčik, J., Lesn,ý, J., Hornik, M., Rajec, P., Kováč, P., et al. (2014). Application of positron emission tomography and 2-[18F]fluoro-2-deoxy-D-glucose for visualization and quantification of solute transport in plant tissues. Chem. Pap. 68, 1463–1473. doi: 10.2478/s11696-014-0609-8

CrossRef Full Text | Google Scholar

Shu, C. W. (2009). High order weighted essentially nonoscillatory schemes for convection dominated problems. SIAM Rev. 51, 82–126. doi: 10.1137/070679065

CrossRef Full Text | Google Scholar

Streun, M., Beer, S., Daemen, J., Dorscheid, R., Erven, A., Jokhovets, L., et al. (2014). “phenoPET: A dedicated PET scanner for plant research based on digital SiPMs (DPCs),” in IEEE NSS/MIC 2014, Conference Recreation (Seattle, WA), M11–M18.

Google Scholar

Streun, M., Borggrewe, K., Chlubek, A., Dautzenberg, M., Degenhardt, C., Dorscheid, R., et al. (2016). “phenoPET – Results from the plant scanner,” in IEEE NSS/MIC/RTSD 2016, Conference Recreation (Strasbourg, SK), M12–M18.

Google Scholar

Suwa, R., Fujimaki, S., Suzui, N., Kawachi, N., Ishii, S., Sakamoto, K., et al. (2008). Use of positron-emitting tracer imaging system for measuring the effect of salinity on temporal and spatial distribution of 11C tracer and coupling between source and sink organs. Plant Sci. 175, 210–216. doi: 10.1016/j.plantsci.2008.03.022

CrossRef Full Text | Google Scholar

Troughton, J. H., Currie, B. G., and Chang, F. H. (1977). Relations between light level, sucrose concentration, and translocation of carbon 11 in Zea mays leaves. Plan Phys. 59, 808–820. doi: 10.1104/pp.59.5.808

PubMed Abstract | CrossRef Full Text | Google Scholar

Tyree, M. T. (1975). “Horwitz-type models of tracer distribution during unidirectional translocation,” in Phloem Transport, NATO Advanced Study Institutes Series (Series A: Life Sciences), Vol. 4, eds S. Aronoff, J. Dainty, P. R. Gorham, L. M. Srivastava, and C. A. Swanson (Boston, MA: Springer), 475–494.

Google Scholar

Van Bel, A. J. E. (2003). The phloem, a miracle of ingenuity. Plant Cell Environ. 26, 125–149. doi: 10.1046/j.1365-3040.2003.00963.x

CrossRef Full Text | Google Scholar

Weisenberger, A. G., Kross, B., Lee, S. J., McKisson, J., McKisson, J. E., Xi, W., et al. (2013). Nuclear physics detector technology applied to plant biology research. Nucl. Instr. Meth. Phys. Res. A 718, 157–159. doi: 10.1016/j.nima.2012.08.097

CrossRef Full Text | Google Scholar