Journal of Theoretical Biology 341 (2014) 131–142
Contents lists available at ScienceDirect
Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi
A class of compartmental models for long-distance tracer transport in plants Jonas Bühler a, Eric von Lieres b,n, Gregor Huber a,nn a b
Forschungszentrum Jülich GmbH, IBG-2: Plant Sciences, 52425 Jülich, Germany Forschungszentrum Jülich GmbH, IBG-1: Biotechnology, 52425 Jülich, Germany
H I G H L I G H T S
Systematic comparison of models regarding quality of fit to tracer transport data. Model class allows construction of numerous models with adjustable complexity. A designated model filter removes irrelevant and redundant models. Different optimal models were found for exemplary data sets.
art ic l e i nf o
a b s t r a c t
Article history: Received 28 March 2013 Received in revised form 11 July 2013 Accepted 16 September 2013 Available online 25 September 2013
Studies of long-distance tracer transport in plants result in spatio-temporal data sets. Compartmental tracer transport models can be used to quantitatively characterize or compare such data sets derived from different experiments. Depending on the specific experimental situation it might be necessary to apply different models. Here, we present a general class of compartmental tracer transport models which allows a systematic comparison of different models regarding the quality of fitting to the experimental data. This model class is defined by a system of partial differential equations (PDEs) for an arbitrary number of parallel compartments with individual transport velocities and numerous lateral exchange connections. A large number of model instances with adjustable complexity can be derived from this model class by permitting only certain model parameters such as flux velocities or exchange rates between compartments to be non-zero. Since some of these models are either inconsistent or redundant we designed a model filter using combinatory rules in order to keep only valid and unique models. A numerical solver for the PDEs was implemented using finite volumes and a weighted essentially nonoscillatory (WENO) scheme. Several candidate models were fitted to experimental data using a Monte Carlo multi-start strategy to approximate the global optimum within a certain parameter space. Analysis of exemplary tracer transport experiments on sugar beet, radish and maize root resulted in different best models depending on the respective data and the required fit quality. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Mechanistic modeling Model filter Model selection Phloem transport PET data analysis
1. Introduction In a recent publication we presented a mechanistic model for long distance tracer transport in plants including an analytic solution of the respective partial differential equations (PDEs) (Bühler et al., 2011). Fitting such a model to data sets from tracer studies (measured e.g. with positron emission tomography, PET, or magnetic resonance imaging, MRI) results in model parameter estimates which characterize the dynamical behavior of the tracer and can be used for quantitative comparison of different data sets, n
Corresponding author. Tel.: þ 49 2461 61 2168. Corresponding author. Tel.: þ 49 2461 61 1503. E-mail addresses:
[email protected] (E. von Lieres),
[email protected] (G. Huber). nn
0022-5193/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jtbi.2013.09.023
e.g. from plant phenotyping studies (Fiorani and Schurr, 2013). For a recent application of this method, see De Schepper et al. (2013). The model used in Bühler et al. (2011) was chosen on the basis of physiological considerations and consisted of three neighboring compartments roughly representing the main functions of vascular transport pathways (e.g. phloem or xylem) and the adjacent tissues. These functions include: (1) axial transport of tracer, simulated by a plug flow in the first compartment, (2) reversible lateral exchange of tracer between the first and second compartment, (3) storage of tracer in the third compartment, and (4) diffusion in axial direction within the first two compartments, see Fig. 1 in Bühler et al. (2011). The capability of the model to represent experimental data acquired by PET was demonstrated. However, it is not clear whether this model is the ‘best’ model according to the principle of parsimony: the most
132
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
resulting models according to their fitting quality, number of model parameters and parameter certainty.
2. A general mechanistic tracer transport model class 2.1. Model class definition
Fig. 1. Sketch of general mechanistic tracer transport model class. The model class consists of N one-dimensional compartments arranged in parallel along a spatial coordinate x. Within the ith compartment, tracer movement is controlled by flux velocities vi. Lateral exchange of tracer between each of the compartments is controlled by individual exchange rates eij. From this model class, single model instances can be defined by allowing a specific set of model parameters to be nonzero. Tracer enters the system by arbitrary initial and boundary conditions.
exact and at the same time most simple data representation possible (Burnham and Anderson, 2002). For example, considering that the model appeared to be over-parameterized when fitted to some of the studied experimental data sets, it is a legitimate question whether or not a simpler model (e.g. with only two compartments) could fit the same data equally well. On the other hand, only slightly more complex models with several independent transport conduits (representing e.g. a bundle of parallel phloem sieve tubes) instead of just a single transporting compartment may fit the data better in other cases. The ‘best’ model choice will generally depend on the investigated plant species or the analyzed plant organ, e.g. stem, root system, or storage organ. In order to systematically compare different models we set up a generalized class of mechanistic transport models, with an arbitrary number of compartments and flux along one or more compartments as well as lateral material exchange between all compartments. While generalizing the model structure, we still restricted ourselves to pure tracer models with parameters constant in time and space. In addition, the same model restrictions hold as described in Bühler et al. (2011), i.e., geometrical and anatomical details are disregarded, transport is allowed only in positive axial direction, and axial transport and lateral exchange in our models average the specific and complex mechanisms of processes in the plant (Patrick, 2013) such as passive or active axial transport and lateral exchange of labeled molecules. In contrast to Bühler et al. (2011), we excluded diffusion in axial direction for simplicity. Each model instance (hereafter just called model) is defined by a specific system of PDEs characterized by a certain set of non-zero parameters. Consequentially, the model of Bühler et al. (2011) now becomes just one member of the much more general set of possible models defined by the model class. We applied finite volumes and a robust and fast WENO scheme (Shu, 1998) for the numerical solution (forward simulation) of the resulting PDEs, since no analytical solution is available for most models. The framework of model fitting and error analysis (bootstrapping) is similar to Bühler et al. (2011). A Monte Carlo multi-start approach was applied to avoid local optima. The introduced set of possible models is very large even for a small number of compartments N. Therefore, a model filter was designed to eliminate redundant and obviously invalid models. We applied a ‘bottom-up’ approach for evaluating models with an increasing number of model parameters k until the quality of the achieved fits ceased to improve significantly. We ranked the
The basic composition of the model class is sketched in Fig. 1. Tracer is present in N parallel compartments. Each compartment can have individual transport velocities in axial direction (spatial coordinate x). Tracer is allowed to transfer between all compartments, determined by respective exchange rates. This model class is defined by a system of partial differential equations ! ∂ ρ ðx; tÞ ∂ ! ¼ AT V ð1Þ ρ ðx; tÞ ∂t ∂x ! The vector ρ ¼ ðρ1 ⋯ρN ÞT denotes the tracer density distribution (arbitrary units) within each compartment at all spatial positions x and time points t. Here and in the following, we omit ! ! space and time dependence in our notation, i.e. ρ ¼ ρ ðx; tÞ and ρi ¼ ρi ðx; tÞ. The coupling matrix A describes tracer entering compartment j coming from compartment i by an exchange rate eij (first term of Eq. (2)) in units of s 1. The second term of Eq. (2) ensures mass conservation by removing exchanged tracer from the respective compartment. The third term of Eq. (2) describes decay of a radioactive tracer by a material specific rate λ, e.g. λ E5.67 10 4 s 1 for the isotope 11C frequently used in PET measurements on plants. 1 0 1 0 N 0 e12 ⋯ e1N 0 C ∑ e1k B Be C Bk¼1 C ⋮ B 21 ⋱ C B C CB ⋱ A¼B C λI B ⋮ C ⋱ eðN 1ÞN A B C N @ @ A e 0 ∑ Nk eN1 ⋯ eNðN 1Þ 0 k¼1
ð2Þ All diagonal elements eii are zero (eii ¼0 for i¼ 1…N) since there is no tracer exchange of one compartment with itself. The matrix V 0 1 v1 0 B C ⋱ ð3Þ V¼@ A 0 vN contains the flux velocities for each compartment. The exchange rates eij and the elements of V are the parameters influencing the spatial and temporal tracer distribution and are therefore referred to as the model parameters. The complexity of the system of PDEs in Eq. (1) depends on how many of the model parameters are nonzero. Each set of non-zero parameters characterizes a specific model instance (model) of the model class. Current tracer imaging methods typically do not allow a separation of tracer in different tissues. Thus, the compartments cannot be individually observed. For that reason, we compare the measured signal to the sum ρtr of the simulated tracer distribution over all compartments N
ρtr ¼ ∑ ρi i¼1
ð4Þ
Tracer entering the system is described either by a spatial initial condition at time t ¼ t ini !
ρ 0 ðxÞ ¼ ! ρ ðx; t ¼ t ini Þ
ð5Þ
or by a time dependent boundary condition at position x ¼ xb !
ρ b ðtÞ ¼ ! ρ ðx ¼ xb ; tÞ
ð6Þ
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
In general, both initial and boundary conditions have to be specified in order to solve the PDEs in Eq. (1). However, typically there is some freedom in the choice of t ini and xb , such that t ini or xb can be set far enough from the signal to allow either the initial or boundary condition to be set to zero. In this paper we define tracer entering the system in the same way as in Bühler et al. (2011) by employing a Gaussian distribution
ρ0 ðxÞ expð ðx x0 Þ2 =2s2 Þ
ð7Þ
as initial condition for compartment 1 at t ¼ t ini , letting the boundary condition be zero far from the measured signal at x ¼ xb . As in Bühler et al. (2011), the Gaussian initial distribution is just an auxiliary function providing an input of tracer suitable for pulse-labeling experiments in cases where the initial or boundary conditions cannot be deduced from the data. The two parameters of the distribution, width s and shift x0 into negative spatial direction, have to be estimated together with the model parameters but are of no further interest. The same initial condition was used for all compartments with non-zero flux velocity (transporting compartments), multiplied by an individual initial condition scaling (ICS) factor for each compartment, setting the amplitudes of the respective initial conditions into relation. We treat these ICS factors as model parameters, since their presence or absence has an impact on the model dynamics independently of the specific functional form of the initial condition. 2.2. Model filtering For a given number of compartments N, the total number of 2 2 possible models is combinatorially given by 2N ð3=2ÞN : 2N N N (possible configurations of the coupling matrix A) times 3 (each compartment can have zero flux velocity, non-zero flux velocity with non-zero ICS factor or flux velocity with zero ICS factor). This high number (cf. Table 1) makes data analysis with all models impractical for N Z3. However, many of the models are either redundant or obviously do not make sense, e.g. if there is transfer from a compartment not containing any tracer. In order to eliminate such invalid or redundant models a model filter was designed as described in Appendix A. This model filter was implemented in MATLAB (R2012a, MathWorks, Inc.) and performed on models with up to N ¼ 4 compartments; the results are depicted in Table 1. In spite of the large reduction, a calculation and comparison of all remaining models is not feasible for N Z3. Comparing models for a given number of compartments N has the disadvantage of a certain imbalance between the possible numbers of flux velocities nv and exchange rates ne since max(nv) ¼N whereas max(ne)¼ N² N. In order to avoid this imbalance, we used an alternative approach, comparing models up to a given number k of model parameters, each of which is an exchange rate, a flux velocity, or an ICS factor: k ¼ ne þ nv þ nICS where nICS is the number of ICS factors. Regarding model filtering, the same procedure can be applied as before (cf. Appendix A). Table 2 shows the number of models for fixed k after model filtering. A total number of possible models for fixed k is not included in Table 2 since this number is not well defined without limiting the allowed number of compartments N. For k up to 5 we Table 1 Total number of possible models for N compartments vs. number of models after filtering. The reduction factor is the ratio of both numbers. N Total number of models Reduced number of models Reduction factor 1 3 2 36 3 1728 4 331776
1 8 139 7984
3 4.5 12.4 41.6
133
Table 2 Total number of models for a given number of model parameters k after filtering. k is the sum of the number of transporting compartments nv, the number of nonzero exchange rates ne and the number of ICS factors nICS . The 48 models for k up to 5 are shown in Fig. 2. k
Filtered number of models
Cumulative number of models
1 2 3 4 5 6 7
1 1 3 9 34 142 719
1 2 5 14 48 190 909
checked manually that all invalid models were sorted out. As the rules for filtering were built on heuristic considerations there is no guarantee that all invalid models were captured for kZ 6. As will be shown this does not matter for the experimental data at hand. The 48 models remaining after filtering for k up to 5 are depicted in Fig. 2 using a simplified design. Model M13 is the model of Bühler et al. (2011), except for the exclusion of diffusion in axial direction. In the following, the numbering of Fig. 2 is referred to whenever a specific model will be mentioned. Each model originates from a parent model by allowing one or more additional parameters to be non-zero. In this way, the 48 models from Fig. 2 can be arranged as a graph where each layer represents models with a specific number of parameters k ¼1…5 as shown in Fig. 3. This structure easily allows recognizing the origin of each model. Knowing the parents of a model can be useful for model fitting to experimental data since parental minima are appropriate start parameter sets.
3. Numerical methods for model fitting The numerical framework to analyze data sets comprises the following steps: (1) model selection, (2) solution of the PDEs (Eqs. (1) and (4)) for a given set of parameters (forward simulation), (3) parameter estimation by global non-linear least square fitting of the model to the data (inverse problem), and (4) calculation of confidence intervals of the estimated parameters. The numerical methods were implemented in MATLAB R2012a using the Optimization Toolbox. The source code can be requested from the authors. 3.1. Solution of PDEs—Forward problem The forward problem refers to solving the model equations for a known set of parameters. Since a general analytic solution of the system of partial differential equations in Eq. (1) is not available numeric solving methods are required. Depending on the model parameters and the smoothness of the initial or boundary conditions the forward simulation can become a stiff problem. The right hand side of the PDEs was first discretized using finite volumes (FV), because FV methods ensure mass conservation and the implementation is straight forward especially at the system boundaries (Hirsch, 2007). The FV method transforms the discretized PDEs into a system of non-linear ordinary differential equations (ODEs). This ODE system was discretized in time using the standard ODE solver ode15s with variable step size and order for stiff problems included in the standard MATLAB library. The FV approach requires the reconstruction of function values at the cell boundaries. A simple first order upwind scheme, though often used for this reconstruction, is known to cause strong numerical diffusion unless extremely fine grids are used, which results in high computing times. In order to achieve a fast but
134
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
Fig. 2. Consecutively numbered models remaining after filtering for a fixed number of parameters k¼ 1…5. Filtering was done according to the rules presented in Appendix A. Numbering of compartments from left to right as in Fig. 1 is omitted for simplicity. Vertical arrows indicate axial fluxes; horizontal arrows denote lateral exchange between compartments. Where different from zero, the ICS factors for the second and any further compartments with transport are fitting parameters, too, denoted by the symbol ■. The symbol * for compartment 1 highlights the fact that the ICS factor was always fixed to unity for this compartment, thus not constituting a fit parameter.
numerically stable solution to Eq. (1) we implemented a third to fifth order weighted essentially non-oscillatory (WENO) scheme (Shu, 1998; Liu et al., 1994). WENO schemes have already been applied to similar systems of PDEs, albeit in different contexts, see e.g. von Lieres and Andersson (2010). Computation times and numerical accuracy of the numerical solver for model M13 were observed to be similar to the analytic solution of Bühler et al. (2011). 3.2. Parameter estimation—Inverse problem The inverse problem refers to estimating unknown model parameters by fitting a model function to given data. Fitting means to
minimize a cost function (Press et al., 2007). In our case, the cost function χ2 is the sum of the squared differences between the simulation ρtr (Eq. (4)) and experimental data at specific points in time and space (the experimental field of view). The optimization problem is non-linear since the cost function depends non-linearly on the model parameters. For the solution of the inverse problem, we used the non-linear least squares solver lsqnonlin included in the MATLAB Optimization Toolbox. lsqnonlin is searching for a local minimum starting from a user supplied start parameter set. Depending on the model as well as the data, the fitted parameter set strongly depends on the choice of the start parameter set. In order to determine an approximation to the global minimum
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
135
Fig. 3. Graphical illustration of model dependencies up to k ¼5. Each model in this family tree extends its respective parents by adding one or more parameters to the set of non-zero parameters.
we implemented a Monte Carlo multi-start fitting approach (Binder and Heermann, 2010; Martí, 2003; Rinnooy Kan and Timmer, 1989): repeated local minimization using a set of random start parameters drawn out of a subspace of the parameter space which can be assumed to contain the global minimum. Applying this procedure the global minimum was not found with absolute certainty but approximated with a probability increasing with the number of repetitions. Certain start parameter combinations were observed not to be in the attracting zone of any minimum, e.g. very low velocities and initial conditions far away from the field of view such that no tracer will enter it in the given time frame. In these cases, all parameters will appear to be insensitive and the parameter estimation will fail to find a minimum other than numeric noise. In order to avoid such cases we pre-filtered the randomly generated start parameter sets. The procedure of pre-filtered multi-start fitting with a repetition number nms consisted of the following six steps: (1) calculating npf nms random start parameter sets where npf is a multiple ensuring a dense sampling of the parameter space, (2) adding already fitted parameter sets from parent models (Fig. 3) to the set of start parameters from step 1 because minima of parent models are known to be good starting points. The additional parameters not included in the parent models were initially set to the lower boundaries of the parameter space (see Table B1), (3) forward simulation of tracer distribution for these start parameters, (4) determining the cost function χ2 for each simulation, (5) sorting the list of start parameter sets according to their respective value of χ2, (6) choosing the nms best start parameter sets for a full non-linear optimization run.
As typical data are smooth and the solution of the forward simulation (Eqs. (1) and (4)) is continuous in the model parameters we expect the cost function to be continuous such that the global minimum has a zone of attraction large enough to be located with a limited number of repetitions nms even if the dimension of parameter space is high. For our simulations we chose npf ¼200 and nms ¼ 200.
3.3. Confidence intervals of fitted parameters After approximating the global minimum, confidence intervals of the fitted parameters were estimated by using a sampling bootstrapping algorithm. A new pseudo data set was created by drawing data points randomly from the original (experimental) data. This altered data set was fitted again by the model using the global minimum from above as a start parameter set resulting in a new estimate of the model parameters. Assuming a Gaussian distribution of the bootstrapped parameters, the width of this distribution can be used as an approximate for the standard deviation (SD) of the global minimum. The 90% confidence interval (CI90) is then given by CI90 ¼1.644 SD (Press et al., 2007; Efron, 1987). Typically the assumption of a Gaussian distribution is valid except for cases where the parameter values approach a boundary of the parameter space (Table B1). In such cases we refrained from calculating the confidence intervals. Bootstrap iterations were repeated until the standard deviations converged. In our case 400 repetitions were performed. Multi-start fitting and data bootstrapping iterations are independent from each other and were parallelized for faster calculation. The basic framework is similar to Bühler et al. (2011).
4. Results For a practical example of model comparison regarding the quality of the fit, we used the same exemplary PET data sets from Jahnke et al. (2009) as analyzed previously in Bühler et al. (2011). Out of the three data sets for sugar beet, radish and maize root, the boundary condition, i.e. the integral amount of tracer entering the field of view, could be approximated by summation of the data along the spatial dimension only in case of sugar beet. For radish and maize, a significant amount of tracer leaves the system at the bottom of the field of view undetected, thus preventing a straight forward determination of the boundary condition from the data without spatial extrapolation of the data. As described in Section 2.1, instead of a boundary condition we chose a Gaussian distribution Eq. (7) as an initial condition. The search for a global minimum (using multi-start fitting as explained in Section 3.2) was successful only if the explored parameter space was restricted in a feasible way. Table B1 in Appendix B shows the boundaries of
136
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
the parameter space applied to our procedure. Using the fitting results of Bühler et al. (2011) as a guideline for the order of magnitude of the parameters, we left enough room for substantial deviations from the previous results. Model fitting was performed for all models up to k¼5. Table 3 shows values of cost function χ2 indicating fit quality for the three data sets of sugar beet, radish and maize. Best fits for each number of model parameters k are highlighted. Values and confidence intervals of the estimated model parameters are listed in Table 4 for models with lowest χ2 in each category k¼2…5. Table 4 also contains values of system loss calculated from the estimated model parameters (for details of loss calculation see Appendix C). For some models loss depends on distance from the point where tracer enters the system, converging to a limit value for large distances, see Fig. C1. For these cases we included the limit values into Table 4. Complete results of the parameter estimation for all models up to k¼ 5 are available as supplementary material. Regarding the sequence of models fitting best for each number of parameters k (Table 3), there is a strong similarity between sugar beet and radish up to k¼4, not surprisingly since both are storage organs. Even so, the overall best fit was obtained with a different model for each data set: models M12 or M13 for sugar beet (preferred to the many k¼5 models with same fit quality because of lower number of parameters), M25 for maize, and M26 for radish. Up to k¼ 4, the best fit improved with increasing k but the improvement was non-uniform: there was a significant reduction of χ2 for k¼2 compared to k¼1, and for k¼4 compared to k¼ 3, but only a comparably smaller improvement of k¼ 3 over k¼ 2, even at the expense of at least one parameter with increased uncertainty (Table 4). k¼5 showed no significant improvement over k¼ 4, some confidence intervals of the parameter estimates became large and an increasing number of parameters was fitted to zero. For maize, the k¼5 model fitting best (M25, see Table 4) did not converge completely as indicated by parameter v2 reaching the upper boundary of 200 mm/min. All this indicates over-parameterization and justifies the termination of the investigation after k¼ 5 since a better fit for higher k becomes very unlikely. Fig. 4C–E shows the deviations between the data and the models of Table 4, again demonstrating that there was no significant improvement of k¼ 5 over k¼4 and k¼3 over k¼ 2. Thus candidates for overall best model are either k¼4 models (best fit) or k¼2 (small number of parameters). For all three data sets model M02 already provided a reasonable fit, as shown exemplarily in Fig. 4A for sugar beet, where the leading
curves were reproduced quite well. Best fitting k¼ 4 models were M12 and M13 for sugar beet as well as radish, M10 for maize. For all three data sets, M12 and M13 delivered almost identical results regarding quality of fit, value of loss and range of confidence intervals. We took M13 as representative (Fig. 4B) since it was the model used in Bühler et al. (2011). For maize, cost function χ2 of M12 and M13 was not quite as low as for M10. Additionally, M12 and M13 for maize suffered from high confidence intervals of over 94% for parameter e12. Comparing the improvement of best fitting k ¼4 models over k¼ 2 (M02), the value of χ2 was much better (by a factor of 6.9) only for sugar beet. For radish and maize the reduction factor was smaller than 2. Confidence intervals were comparable to those of M02 for sugar beet and maize, but much higher for radish. In conclusion, the most probable choice for the best data representation will be model M12 or M13 for sugar beet, M02 for radish because of low confidence intervals and few parameters, M02 or M10 for maize depending on purpose.
5. Discussion 5.1. Identification of the ‘best’ model The criterion for selecting the overall ‘best’ model is an optimal balance between minimal cost function χ2, minimal number of model parameters and maximal certainty of the parameter estimates. The manual model ranking described in the previous section was based on this idea. For the data sets analyzed here, an automated model ranking using the information criteria AIC and AICc (Burnham and Anderson, 2002; Akaike, 1981) did not add more insight as compared to the manual model ranking (data not shown). Nevertheless, for future analysis of large numbers of data sets automated model selection and model inference can potentially be very useful. From the fitting results for the data sets of sugar beet, radish and maize root (Table 4) it is obvious that there is not one optimal model for all experimental situations. Different plant species, different plant parts being observed, different measuring protocols and methods lead to different characteristics of the data sets and are ‘best’ described by different models. However, there are some models which are most probably never optimal, though all
Table 3 Results of model fitting. Models with up to k ¼5 model parameters were fitted to the three PET data sets from Jahnke et al. (2009). The indicator of fit quality, cost function χ2, is the sum of the squared differences between model and data. Candidates for ‘best’ model for each k are highlighted by grey background.
Sugar beet Maize root Radish bulb
k ¼1 M01 0.5595 0.3189 1.2206
k¼ 2 M02 0.0406 0.1219 0.0860
k¼ 3 M03 0.1943 0.1363 0.4766
M04 0.0276 0.1107 0.0844
M05 0.0251 0.1090 0.0829
k ¼4 M06 0.0276 0.1107 0.0844
M07 0.0227 0.0709 0.0547
M08 0.0251 0.1052 0.0829
M09 0.0074 0.0817 0.0566
M10 0.0276 0.0649 0.0543
M11 0.0276 0.1050 0.0573
M12 0.0059 0.0750 0.0477
M13 0.0059 0.0750 0.0477
M14 0.0251 0.1091 0.0829
χ2
Sugar beet Maize root Radish bulb
k ¼5 M15 0.1234 0.1017 0.2921
M16 0.0188 0.0709 0.0514
M17 0.0065 0.0704 0.0566
M18 0.0276 0.0649 0.0506
M19 0.0251 0.1052 0.0829
M20 0.0074 0.0709 0.0480
M21 0.0182 0.0709 0.0509
M22 0.0227 0.0649 0.0543
M23 0.0227 0.0709 0.0466
M24 0.0102 0.0709 0.0487
M25 0.0059 0.0608 0.0477
M26 0.0059 0.0649 0.0456
M27 0.0225 0.0971 0.0573
M28 0.0059 0.0646 0.0477
χ2
Sugar beet Maize root Radish bulb
M29 0.0276 0.0649 0.0465
M30 0.0064 0.0817 0.0566
M31 0.0065 0.0817 0.0566
M32 0.0074 0.0649 0.0469
M33 0.0217 0.0906 0.0573
M34 0.0059 0.0700 0.0477
M35 0.0276 0.0634 0.0461
M36 0.0065 0.0800 0.0552
M37 0.0276 0.1019 0.0575
M38 0.0226 0.0634 0.0462
M39 0.0059 0.0750 0.0477
M40 0.0059 0.0750 0.0477
M41 0.0059 0.0750 0.0477
M42 0.0059 0.0750 0.0477
χ2
Sugar beet Maize root Radish bulb
M43 0.0059 0.075 0.0477
M44 0.0251 0.1090 0.0829
M45 0.0059 0.0744 0.0477
M46 0.0059 0.0744 0.0477
M47 0.0059 0.0744 0.0477
M48 0.0251 0.1094 0.0829
χ2
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
137
Table 4 Best model fits. Fitting parameters and their respective 90% confidence intervals for the best models in each category k ¼ 2…5 for the three PET data sets from Jahnke et al. (2009). Model instance M13 is the model used in Bühler et al. (2011). n.a. indicates that the parameter estimation reached a boundary or the confidence interval exceeded 100%, violating the assumption of Gaussian shape of the parameter distribution. Sugar beet
k¼ 2 M02 Fit
CI90 [%]
k¼ 3 M05 Fit
CI90 [%]
k ¼4 M13 Fit
1.24 7
8.9
1.16
7
7.1
1.99
7
6.1
0.077 7
12.7
0.071
7
10.8
0.338
7
11.8
0.002
7
37.3
0.025
0.135 7 0.061 7 40.9 7 0.0059
CI90 [%]
k¼ 3 M05 Fit
k ¼4 M10 Fit
1.40 7
3.9
2.34
7
8.2
0.015 7
16.9
0.259
7
27.4
0.232
7
14.5
Parameter
Unit
v1 v2 e12 e13 e21 e23 Loss χ2
mm/min mm/min min 1 min 1 min 1 min 1 %/cm a.u.
46.4 7 0.041
Parameter
Unit
k¼ 2 M02 Fit
v1 v2 e12 e13 e21 e31 e32 Loss χ2
mm/min mm/min min 1 min 1 min 1 min 1 min 1 %/cm a.u.
10.2 7 0.122
Parameter
Unit
k¼ 2 M02 Fit
v1 v2 e12 e13 e21 e23 Loss χ2
mm/min mm/min min 1 min 1 min 1 min 1 %/cm a.u.
Maize root
Radish bulb
3.6
CI90 [%]
8.3 9.7 3.1
CI90 [%] 7 7 7 7
7
4.9 4.6 17.2 19.8
0.109
22.5 0.065
CI90 [%]
k¼ 3 M05 Fit
CI90 [%]
k ¼4 M13 Fit
1.28 7
6.6
1.24
7
6.2
2.16
7
23.4
0.067 7
11.2
0.065
7
10.8
0.423
7
54.7
0.001
7
88.5
0.237 0.071 36.4 0.048
7 7 7
33.1 22.1 4.0
41.0 7 0.086
13.4
1.72 1.31 0.008 0.036
CI90 [%]
4.2 0.083
obviously invalid models have already been filtered out. In order to fit specific tracer data at least one of the following features are required: (1) storage possibility, (2) back-transfer into a transporting compartment, and (3) more than one transporting compartment with different transport velocities. For example, models with a higher number of parameters but pure circling transfer processes, like M14 and M48, failed to describe any of the data better than M05. As another example, no model without a storage compartment was able to reproduce the sugar beet data sufficiently well (e.g. M37, M44). Some of the well-fitting models mimicked a storage compartment by setting one exchange parameter to zero (e.g. models M39–M43). Sugar beet and radish are storage organs while maize roots are known to mainly transport the tracer (Jahnke et al., 2009). Consistent with these properties, the loss values of the best models for sugar beet and radish are high (ca. 41%/cm) whereas for maize the loss is only approximately half of that (Table 4). For maize roots the data at hand also integrated tracer from small side roots which act as little sinks distributed along the main transport pathway and could not be disentangled in the original 3D PET data sets. For this reason, we expect the loss of transporting maize roots to be even smaller. Regarding storage it is apparently rather unimportant where the storage compartment is fed from, as for example models M12 and M13 delivered almost identical results. For all models of Table 4, velocity and loss have the smallest confidence intervals indicating that these parameters are identified best. In light of the high certainty of the estimates, the difference of velocities and loss
k¼5 M26 Fit
CI90 [%]
0.0 1.95 0.129 0.061 0.322
7 7 7 7 7
n.a. 5.5 17.1 8.2 11.5
40.8 0.0059
7
3.1
k¼5 M25 Fit
CI90 [%]
1.90 200
7 7
15.5 n.a.
0.117
7
66.2
0.171 0.021
7 7
34.2 26.7
10.2 0.061
CI90 [%]
k¼5 M26 Fit
CI90 [%]
1.58 1.12 0.003 0.113 0.033
7 7 7 7 7
10.3 23.6 n.a. 19.3 n.a.
40.5 0.046
7
8.5
values between different models for k ¼2 and k ¼4 is significant, emphasizing the need for careful selecting a best model out of the set of valid models, particularly if one wants to quantitatively compare different data sets. In our analysis of maize data, models M02 and M10 remained as candidates for ‘best’ model without a clear decision. This shows that it can be difficult to identify a unique ‘best’ model based on the criteria mentioned above (minimal cost function, minimal number of parameters, and maximal certainty of estimate). In such cases additional factors will be important to choose a model for characterizing the data set: the purpose of the experiment, the required fit quality, the physiological relevance of the candidate models, and the capability of the models to fit to data from repeated or similar experiments, as well as to data before and after a treatment. Yet, the possible deductions are restricted by the informational content of the data. For example, regarding the maize data, a physiological situation can be possible where tracer from one leaf is distributed into phloem pathways with different propagation velocities through vascular connections. This situation is represented by, e.g., model M10 and M18. Our results for maize data show that adding further transport pathways does not necessarily lead to a better fit (M18 compared to M10) and that models with several transporting compartments can neither be clearly identified as ‘best’ model (compared to M02) nor can they be ruled out. Hence, future studies will include experimental design strategies for maximizing the information contained in the measured data
138
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
Fig. 4. Residuals of best model fits. Fitting results of best models for k ¼ 2…5 for PET data from Jahnke et al. (2009). (A) and (B) PET data (symbols) and fitted data (lines) for sugar beet using models M02 (most simple) and M13 (best χ2 fit). (C)–(E) Differences of sugar beet, maize and radish data to the respective fits for the best models of Table 4.
sets, which will improve not only the parameter confidence but also the stability of model ranking. 5.2. Possible model extensions and performance improvements Remaining deviations between the best fitting model and the data result from either missing the global minimum or from features in the data which are not covered by the general
assumptions of our model class. Given the fact that the transporting tissues of plants are not spatially homogenous it is clear that our assumption of spatially constant transport properties is a strong simplification and might result in respective model errors. Moreover, our assumption of temporal invariance would be violated in experimental situations where manipulations of the plant are conducted within the measurement. Thus, it would be worthwhile to extend the model by allowing model parameters to
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
be a function of the spatial and/or temporal coordinate. Of course this would generate additional parameters that need to be fitted if they cannot be measured, possibly leading to over-parameterization. The robustness of the parameter estimation with variable velocities and exchange rates is subject of current investigations. As already mentioned by Bühler et al. (2011) for the data at hand, fitting additional diffusion parameters did not result in a better fitting quality but generally lead to less certain parameter estimations. For this reason, in the model class presented here diffusion was neglected but could easily be included in the governing equations: ! ∂ ρ ðx; tÞ ∂ ∂2 ! ¼ AT V þ D 2 ρ ðx; tÞ ∂t ∂x ∂x Using individual diffusion coefficients di for each compartment would also require an adaptation of the model filter and still would result in a tremendously enlarged set of valid and unique models. As mentioned in Thorpe et al. (1998) for common 11C leaf labeling studies the back transport of tracer within the xylem can be neglected due to the fast tracer decay. For this reason, we used only positive velocities in our models. On the other hand, depending on the plant architecture the local phloem is known to allow simultaneous bidirectional transport in parts of the plant (Thorpe et al., 1998, 2011) as can be deduced from experiments in which leaves of different nodes were labeled and the respective tracer distribution in the plant was measured (Jahnke et al., 1998). The analysis of such experiments would require models also allowing negative transport velocities including a more intricate definition of boundary conditions. Searching the global minimum causes a high computational effort, and the applied Monte Carlo strategy cannot guarantee global optima. Especially for some k ¼5 models, a reasonable minimum could only be found by using already known minima from parent models as starting values for the optimization. With increasing number of parameters k the investigated parameter space becomes larger and an optimal density of the Monte Carlo multi-start parameter sets is difficult to determine. On the other hand, our approach of pre-scanning the parameter space as described in Section 3.2 helped to speed up the search, and the consistency of our results (Table 3) indicates that the remaining deviations between model and data were caused by model errors, not by missing the global minimum.
139
the data sets of sugar beet and radish if starting with only N¼2 compartments: in this case M19 is the complete model containing all possible parameters (k¼5 including ICS factors). Out of this set of parameters all additional parameters of M19 compared to M02 were estimated as (almost) zero, thus effectively reducing the model. We tried the same approach for N¼ 3 compartments (k¼ 11 parameters) but failed in determining an acceptable minimum. Apparently the combination of over-parameterization and a larger parameter space lead to a very small attraction zone of the (global) minimum so that the optimizer required starting values from the next vicinity of the minimum. We conclude that the ‘top-down’ approach is not a viable option. Consequently, the ‘bottom-up’ approach presented in this paper is the most comprehensive and systematic way of determining the optimal model for a given data set.
Acknowledgements The authors gratefully acknowledge Siegfried Jahnke for providing the data and critical reading of the manuscript, Hanno Scharr and Ralf Metzner for helpful discussions, and Nico Becker for implementing the WENO scheme.
Appendix A. Rules for model filtering In order to be able to eliminate invalid and redundant models we first sorted the compartments according to whether they allow tracer transport (i.e. have non-zero flux velocities) or not: nv rN denoting the number of transporting compartments, we assumed the first nv diagonal elements of the flux velocity matrix V (Eq. (3)) to be nonzero. Without loss of generality the ICS factor of compartment 1 can be set at unity since ICS is a relative scaling between the compartments. The remaining nICS r nv 1 non-zero ICS factors are assumed to be arranged in continuous order, starting with compartment 2. Each model can then be characterized by nv, nICS , and by the structure of A (Eq. (2)). The filtering of valid models from the model class was done iteratively for n¼1 … N by applying rules to the configuration of the n n coupling matrix An . The rules are based on considerations regarding the potential presence or absence of tracer
5.3. Comparison to other approaches Our overall approach is somewhat similar to the two step procedure of input-output analysis (Minchin et al., 1996) with a selection of the ‘best’ model followed by detailed analysis of the data using this model. The main advantages of the method presented here are the possibility to fit data from pulse labeling experiments without having to provide boundary conditions and the systematic and straightforward error analysis. Also the presented mathematical model, albeit highly abstracted, is closer to the underlying physical system than a purely data driven approach and therefore more beneficial for the physiological interpretation of the model parameters. In the previous paper Bühler et al. (2011) adding or removing parameters to or from the initial model was suggested for adapting it to a specific experimental situation. The systematic approach to model comparison presented here, i.e. the combination of a generalized model class with careful model filtering, implements this extension in a very general way, ensuring that every possible meaningful model variation is covered. Alternatively, one could think of a ‘top-down’ approach to model selection: starting the fitting procedure with a model containing the full range of parameters for a sufficiently high number of compartments, followed by successive deletion of those parameters which have been estimated as zero or stay undetermined (i.e. by showing high confidence intervals). This works quite well for
Fig. A1. Model filter. Illustration of restrictions for structure of coupling matrix ATn . n is the total number of compartments, nv is the number of transporting compartments. In order to filter out obviously invalid models, we kept only models where rules (1)–(5) apply to the structure of An. Diagonal elements do not need to be considered as evident from Eq. (2). Rule (5) is not shown in the illustration.
140
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
in the compartments. The first four of the rules are illustrated in Fig. A1. We kept only models where (1) there is tracer present in all compartments because any model with empty compartments is covered by a smaller n. This rule implies that each of the columns of ATn for non-transporting compartments and for transporting compartments with zero ICS factor contains at least one nonzero element; (2) there is only one storage compartment. Two or more storage compartments could be combined into one. This rule implies that all but maximally the last one of the rows of ATn for nontransporting compartments contain at least one nonzero element; (3) there is at least one transfer from a transporting into a nontransporting compartment—otherwise non-transporting compartments are dispensable since no tracer is entering. Without
Table B1 Boundaries of model parameter space for optimization and multi-start intervals used for the analysis of PET data of sugar beet, radish, and maize from Jahnke et al. (2009). Parameter Unit
vi eij ICSi x0 s
mm/min 1 min 1 Dimensionless mm mm
Lower Upper Lower boundary boundary boundary for multistart fitting
Upper boundary for multistart fitting
1e 6 1e 8 1e 8 20 3
5 1 1e3 200 30
200 1e3 1e3 300 100
0.5 1e 8 0.5 20 5
loss of generality the first transfer is assumed to be from the first compartment or from the first transporting compartment without ICS factor into the first or second non-transporting compartment, as this always can be achieved by interchanging some compartments; (4) if there is more than one non-transporting compartment, there is at least one transfer back into a transporting compartment—otherwise the non-transporting compartments could be combined into one. Without loss of generality the first back-transfer can be assumed to be from the first nontransporting compartment; (5) as a generalization of rules (3) and (4), there is no subset of non-transporting compartments with exchange amongst each other but (i) no tracer entering or (ii) no back-transfer into a transporting compartment. In mathematical terms, even after any permutation of rows and columns (nv þ 2) … n, An is not allowed to have the form of a block diagonal matrix consisting of two blocks M1 and M2 where M1 ¼ m1 m1 has size m1 Znv and M2 ¼m2 m2 has size m2 ¼(n m1) Z2 ! M1 0 ATn a 0 M2 In addition, there is no subset of compartments without ICS factor (transporting or non-transporting) with no tracer entering, implying that even after permutation of rows and columns ðnICS þ 2Þ… n, ATn is not allowed to have the form of a block lower triangular matrix. ! M1 0 ATn a M3 M2 where M1 ¼ m1 m1 has size m1 ZnICS þ1, M2 ¼m2 m2 has size m2 ¼(n m1)Z 2 and M3 ¼m2 m1.
Fig. C1. Spatial dependencies of loss. Loss depending on position dx1 for models up to k ¼ 5 with two transporting compartments both contributing (by direct or indirect lateral) to permanent storage of tracer. For maize some exchange parameters of models M20, M24, M26 and M32 were fitted to zero, reducing these models to simpler models. In these cases, the respective loss-equations are no longer valid.
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
In a final step, redundant models were eliminated using a graph theoretic approach (Bondy and Murty, 2008). Replacing all non-zero entries of A and V with 1 and setting the diagonal elements of A to zero, the matrix (Aþ2V) can be interpreted as adjacency matrix of a finite, directed graph with no parallels, N vertices, nv self-loops and ne arcs where ne is the number of nonzero elements in A. Thus each model can be characterized by a specific graph. In this context, redundancy reduction means searching for graph isomorphism: keeping only one representative of every isomorphism class eliminates all redundant graphs. The graph isomorphism problem is known to be computationally extremely expensive (Babai, 1995). We facilitated the task by first classifying a model according to its vertex degrees. Then the search for isomorphism was done by comparing all possible permutations of vertices only amongst graphs with equal degree sequence.
Appendix B. Restriction of parameter space For the three data sets of sugar beet, radish, and maize root from Jahnke et al. (2009) feasible restriction of the model parameter space is shown in Table B1.
Appendix C. Loss calculation System loss is the ratio of tracer permanently fixed along a region of interest (i.e. a certain spatial distance x2 x1 ) to the tracer entering this region. Loss can be calculated in the same way as outlined in Bühler et al. (2011): setting the decay rate λ ¼ 0 and assuming a constant flow of tracer j0 into compartment 1 at position x0 (boundary condition ρ1 ðx0 Þ ¼ j0 ) leads to a steady state at t-1 in all compartments not accumulating tracer. Thus Eq. (1) reduces to a system of first order linear differential equations which can be solved analytically. Resulting loss is of the form Loss ¼ 1 expð ϑðx2 x1 ÞÞ
ðC:1Þ
for many models where there is loss at all, i.e. where there is a storage compartment. Here ϑ includes the relevant transport and exchange properties of the respective model: e12 e13 ; ϑM07 ¼ ϑM12 ¼ ϑM38 ¼ ; v1 v1 e þe e12 e23 ; ϑM10 ¼ ϑM22 ¼ 12 13 ; ϑM13 ¼ v1 v1 ðe21 þ e23 Þ e e þ ðe12 e23 =ðe21 þ e23 ÞÞ ϑM35 ¼ ϑM45 ¼ 14 ; ϑM43 ¼ 13 ; v1 v1 e13 e24 e13 e34 ; ϑM47 ¼ ϑM46 ¼ v1 ðe21 þ e24 Þ v1 ðe32 þ e34 Þ
ϑM02 ¼
For models M09, M20, M24, M26, M32, M34 and M36 loss depends explicitly on position, not just distance x2 x1 as in
Eq. (C.1). These are the models up to k ¼5 with two transporting compartments both contributing (by direct or indirect lateral transfer) to permanent storage of tracer. Using the abbreviations Δx1 ¼ x1 x0 and Δx2 ¼ x2 x1 the expressions for loss of these models are with the abbreviation R ¼ ððe12 v2 þ e21 v1 Þ2 þ e13 v2 ðe13 v2 þ 2e12 v2 2e21 v1 ÞÞ1=2 . LossM32 ¼ 1
LossM34 ¼ 1
e23 vv12 e13 exp e12 vþ1 e13 ðΔx1 þ Δx2 Þ e12 exp ev232 ðΔx1 þ Δx2 Þ e23 vv12 e13 exp e12 vþ1 e13 Δx1 e12 exp ev232 Δx1 e
e12 v2 þ e13 v2 þ e21 v1 2v1 v2
n
Δx2 e12 v2 þ e13 v2 þ e21 v1 R e12 v2 þ e13 v2 þ e21 v1 þ Re
e12 v2 þ e13 v2 þ e21 v1 R e12 v2 þ e13 v2 þ e21 v1 þ Re
R v1 v2
R ð2 2v1 v2
Δx 1
Δx1 þ Δx2 Þ
R
e2v1 v2
Δx2
o
1
2
with the abbreviation R ¼ ððe12 v2 þ e21 v1 Þ þ e23 v1 ðe23 v1 2e12 v2 þ 2e21 v1 ÞÞ1=2 . v2 exp ev242 Δx1 ev131 ðΔx1 þ Δx2 Þ ee13 exp ev242 Δx2 24 v1 LossM36 ¼ 1 v2 exp ev242 ev131 Δx1 ee13 24 v1 Fig. C1 shows the loss of these models using the parameter estimates for the three data sets of sugar beet, maize root and radish (see tables in Supplementary material) as function of Δx1 for fixed Δx2 ¼ 1 cm. As can be seen from the above formulae the spatial dependence of loss results from (at least) two exponential functions with different drop in spatial direction. Contributions with fast drop will become less important with longer distance to the origin (here denoting the point where tracer enters the system). In the end, i.e. far away from the origin, only the contribution with the slowest drop will dominate, leading to a state where loss becomes independent of the space coordinate. Eq. (C.2) exemplarily shows how the value of loss in the limit of infinite distance to the origin depends on the relation of the exponents governing the spatial tracer distribution, reducing to the format of Eq. (C.1): 8 e > for ev232 4 ev121 < 1 exp v121 Δx2 lim ðLossM09 Þ ¼ ðC:2Þ > Δx1 -1 for ev121 4 ev232 : 1 exp ev232 Δx2 For our data sets, in some cases the limit was reached only within long distances of Δx1 4 1 m, see Fig. C1. Even so, for simplicity we took the limit as representative for the system behavior in these cases.
Appendix D. Supplementary material Supplementary Tables S1–S3 include fitting parameters and their respective deviations for all models up to k = 5 for the sugar beet, maize and radish PET data, respectively, from Jahnke et al. (2009). The calculation of the standard deviations was done
ev121 ðΔx1 þ Δx2 Þ ev121 exp ev232 ðΔx1 þ Δx2 Þ LossM09 ¼ 1 e23 e12 e12 e23 v2 exp v1 Δx1 v1 exp v2 Δx1 n o e23 e12 e12 e23 v2 e23 v2 exp v1 ðΔx1 þ Δx2 Þ þ v1 e12 v1 BC 2 1 exp v2 ðΔx1 þ Δx2 Þ n o LossM20 ¼ 1 e23 e12 e12 e23 v2 e23 v2 exp v1 Δx1 þ v1 e12 v1 BC 2 1 exp v2 Δx1 exp ev131 Δx1 þ Δx2 þ vv21 BC 2 exp ev232 Δx1 þ Δx2 LossM24 ¼ 1 exp ev131 Δx1 þ vv21 BC 2 exp ev232 Δx1 n o e v þe v þe v R 12 2 2v13 v2 21 1 Δx2 e12 v2 þ e13 v2 e21 v1 þ R 2v Rv ð2Δx1 þ Δx2 Þ v2 þ e13 v2 e21 v1 R 2v1 v2 Δx2 1 2 1 2 e ee12 e e12 v2 þ e13 v2 þ e21 v1 þ Re 12 v2 þ e13 v2 þ e21 v1 R LossM26 ¼ 1 R e12 v2 þ e13 v2 e21 v1 þ R v1 v2 Δx1 v2 þ e13 v2 e21 v1 R ee12 e12 v2 þ e13 v2 þ e21 v1 þ Re 12 v2 þ e13 v2 þ e21 v1 R e23 v2 exp
141
142
J. Bühler et al. / Journal of Theoretical Biology 341 (2014) 131–142
regardless of the assumption of a Gaussian distribution of the bootstrapped parameters. Supplementary data associated with this article can be found in the online version at http://dx.doi. org/10.1016/j.jtbi.2013.09.023. References Akaike, H., 1981. Likelihood of a model and information criteria. Journal of Econometrics 16, 3–14. Babai, L., 1995. Automorphism groups, isomorphism, reconstruction. Handbook of Combinatorics, vol. 2. Elsevier, Amsterdam. Binder, K., Heermann, D.W., 2010. Monte Carlo Simulation in Statistical Physics— An Introduction. Springer, Heidelberg. Bondy, J.A., Murty, U.S.R., 2008. Graph Theory. Graduate Texts in Mathematics, 244. Springer, New York. Bühler, J., Huber, G., Schmid, F., Blümler, P., 2011. Analytical model for long-distance tracer-transport in plants. Journal of Theoretical Biology 270, 70–79. Burnham, K.P., Anderson, D.R., 2002. Model Selection and Multimodel Inference. Springer, New York. De Schepper, V., Bühler, J., Thorpe, M., Roeb, G., Huber, G., van Dusschoten, D., Jahnke, S., Steppe, K., 2013. 11C-PET imaging reveals transport dynamics and sectorial plasticity of oak phloem after girdling. Frontiers in Plant Science 4, 200. Efron, B., 1987. Better bootstrap confidence intervals. Journal of the American Statistical Association 82, 171–185. Fiorani, F., Schurr, U., 2013. Future scenarios for plant phenotyping. Annual Review of Plant Biology, 64. Hirsch, C., 2007. Numerical Computation of Internal & External Flows. ButterworthHeinemann, Oxford. Jahnke, S., Schlesinger, U., Feige, B.G., Knust, E.-J., 1998. Transport of photoassimilates in young trees of Fraxinus and Sorbus: measurement of translocation in vivo. Botanica acta: Berichte der Deutschen Botanischen Gesellschaft 111, 307–315.
Jahnke, S., Menzel, M.I., van Dusschoten, D., Roeb, G.W., Bühler, J., Minwuyelet, S., Blümler, P., Temperton, V.M., Hombach, T., Streun, M., Beer, S., Khodaverdi, M., Ziemons, K., Coenen, H.H., Schurr, U., 2009. Combined MRI-PET dissects dynamic changes in plant structures and functions. The Plant Journal 59, 634–644. Liu, X.D., Osher, S., Chan, T., 1994. Weighted essentially non-oscillatory schemes. Journal of Computational Physics 115, 200–212. Martí, R., 2003. Multi-start methods. In: Glover, F., Kochenberger, G.A. (Eds.), Handbook of Metaheuristics, vol. 57. Springer, New York, pp. 355–368. Minchin, P.E.H., Lees, M.J., Thorpe, M.R., Young, P.C., 1996. What can tracer profiles tell us about the mechanisms giving rise to them? Journal of Experimental Botany 47, 275–284. Patrick, J.W., 2013. Fundamentals of phloem transport physiology. In: Thompson, G.A., van Bel, A.J.E. (Eds.), Phloem: Molecular Cell Biology, Systemic Communication, Biotic Interactions. John Wiley. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 2007. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, New York. Rinnooy Kan, A.H.G., Timmer, G.T., 1989. Global optimization. In: Nemhauser, G.L., et al. (Eds.), Handbooks in Operations Research and Management Science, vol. 1. Elsevier, pp. 631–662. Shu, C.-W., 1998. Essentially non-oscillatory and weighted essentially nonoscillatory schemes for hyperbolic conservation laws. In: Quarteroni, A. (Ed.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, vol. 1697. Springer, Heidelberg, pp. 325–432. Thorpe, M.R., Walsh, K.B., Minchin, P.E.H., 1998. Photoassimilate partitioning in nodulated soybean I. 11C methodology. Journal of Experimental Botany 49, 1805–1815. Thorpe, M.R., Lacointe, A., Minchin, P.E.H., 2011. Modelling phloem transport within a pruned dwarf bean: a 2-source-3-sink system. Functional Plant Biology 38, 127–138. von Lieres, E., Andersson, J., 2010. A fast and accurate solver for the general rate model of column liquid chromatography. Computers and Chemical Engineering 34, 1180–1191.
