
Methodology

Open

 Published:
Tailormade composite functions as tools in model choice: the case of sigmoidal vs bilinear growth profiles
Plant Methodsvolume 2, Article number: 11 (2006)
Abstract
Background
Roots are the classical model system to study the organization and dynamics of organ growth zones. Profiles of the velocity of root elements relative to the apex have generally been considered to be sigmoidal. However, recent highresolution measurements have yielded bilinear profiles, suggesting that sigmoidal profiles may be artifacts caused by insufficient spatiotemporal resolution. The decision whether an empirical velocity profile follows a sigmoidal or bilinear distribution has consequences for the interpretation of the underlying biological processes. However, distinguishing between sigmoidal and bilinear curves is notoriously problematic. A mathematical function that can describe both types of curve equally well would allow them to be distinguished by automated curvefitting.
Results
On the basis of the mathematical requirements defined, we created a composite function and tested it by fitting it to sigmoidal and bilinear models with different noise levels (MonteCarlo datasets) and to three experimental datasets from roots of Gypsophila elegans, Aurinia saxatilis, and Arabidopsis thaliana. Fits of the function proved robust with respect to noise and yielded statistically sound results if care was taken to identify reasonable initial coefficient values to start the automated fitting procedure. Descriptions of experimental datasets were significantly better than those provided by the Richards function, the most flexible of the classical growth equations, even in cases in which the data followed a smooth sigmoidal distribution.
Conclusion
Fits of the composite function introduced here provide an independent criterion for distinguishing sigmoidal and bilinear growth profiles, but without forcing a dichotomous decision, as intermediate solutions are possible. Our function thus facilitates an unbiased, multipleworking hypothesis approach. While our discussion focusses on kinematic growth analysis, this and similar tailormade functions will be useful tools wherever models of steadily or abruptly changing dependencies between empirical parameters are to be compared.
Background
Kinematic growth analysis aims at the quantitative description of spatial growth patterns to provide a basis for the study of developmental mechanisms [1, 2]. As the term kinematic indicates, this approach focuses on the movement of parts of a growing organ relative to each other. The concepts for kinematic growth analysis were laid out halfacentury ago for the root, which was assumed for simplicity to grow unidirectionally (i.e. pure elongation or axial growth) [3–5]. The assumption of pure elongation keeps the mathematics manageable. Consequently, the vast majority of studies applying kinematic growth analysis to physiological problems have focused on pure elongation in suitable organs such as roots [6–8] and grass leaves [9–11].
The key parameter is the velocity field, the spatial distribution of the velocities at which these displacements occur. For unidirectional expansion, the velocity field reduces to a velocity profile. The derivative of the velocity profile has often been referred to as the relative elemental elongation rate. This rate is relative because it is a measure of growth that is independent of the size of the growing entity, with the dimension of reciprocal time, and it is elemental because it represents a calculusbased description of infinitesimal elements of tissue [2]. However, an elemental rate is by definition a relative rate. To avoid this redundancy, we drop the relative and use simply elemental growth (or elongation) rate to refer to the spatial derivative of a velocity profile. This usage also helps to avoid confusion between an elemental growth rate, which describes motion within a spatial system of reference, and the relative growth rate, a wellestablished concept in classical growth analysis which describes changes in size over time. Here we are concerned primarily with the former type of rate, although the analytic tool we introduce might prove useful in classical growth analysis as well.
Velocity profiles along growing root tips and leaves have been widely reported to be sigmoidal; concomitantly, the corresponding elemental growth rate profiles reported were bellshaped, with a single, smooth peak. Recently, for the root, determination of velocity profiles at greatly increased temporal and spatial resolution has produced distributions that appeared to be bilinear ([12]; for a comparison of sigmoidal and bilinear velocity profiles, see Fig. 1). Accordingly, the corresponding elemental growth rate profiles had "stepstool"shapes, showing two relatively stable plateaus (Fig. 1B). Intriguingly, profiles along growth zones of anatomic parameters such as cell length, which under steady state conditions and in the absence of cell division are geometrically similar to the corresponding velocity profiles [13], have occasionally been interpreted as being bilinear [14–16]. The recent report provides experimental support to the suspicion that sigmoidal velocity profiles might, at least in some cases, be artifacts resulting from measurement error, averaging, or insufficient spatial or temporal resolution [12].
The distinction between sigmoidal and bilinear velocity profiles is important biologically. As a root cell traverses a sigmoidal growth zone, it increases and then decreases its elongation rate steadily and smoothly from one end of the zone to the other. In contrast, as a cell traverses a bilinear growth zone, it elongates at one steady rate for part of the zone, rapidly increases to a new rate for the rest of the growth zone, and then stops (Fig. 1C). In the sigmoidal case, regulation of elongation rate is expected to be continuously variable whereas in the bilinear case, the regulation should establish two distinct rates of elongation as well as the positions where the transitions occur. In reality, no growth zone will be exactly bilinear because elemental growth rates cannot be perfectly constant and change instantaneously in a mathematical sense. It is an open question whether some growth zones are purely sigmoidal but it appears possible that all or most root growth zones have a mixture of sigmoidal and bilinear characteristics.
The analytical power of kinematic growth analysis rests on the fact that knowledge of the velocity field enables one not only to calculate the local rates of deposition of any parameter of interest, whether it be cells, cell wall material, solutes, or water, but also to compute timecourses of these and other parameters as experienced by cells that traverse the growth zone [17–20]. To exploit this analytical power in full, the profiles of velocity and elongation rate need to be rendered as continuous functions. For this reason, kinematic growth analysis is associated with curvefitting, where some function (or group of functions) is fitted to the raw velocity data.
Because at present there is no mechanistic model for the regulation of the velocity field within a growth zone, the choice of a function to fit is arbitrary. As velocity profiles generally resemble sigmoid curves, authors have applied sigmoid functions previously established in classical growth analysis (such as the Gompertz or Richards functions; [21]) as well as versions developed specifically for kinematic growth analysis [22, 23]. However, such functions will smooth out regions of the profile whose behaviour departs from the sigmoid. Alternatively, one may use a piecewise approach where polynomials are fit to small, overlapping subsets of the data [2, 7]. This approach excels at capturing local behaviour but nevertheless smoothes abrupt transitions [12] and can be difficult to apply to noisy data. Alternatively, van der Weele et al. [12] fitted linear regression lines to the velocity data, but this approach yields profiles with discontinuous derivatives and allows no possibility for any sigmoid character. Describing stem elongation over time, Fisher et al. [24] introduced a 3phasic equation which avoided discontinuous derivatives, but which was based on the assumption that a linear growth phase did in fact exist.
To analyze velocity profiles without forcing them to be either sigmoids or to contain a straight line, we formulated a function which describes sigmoidal and bilinear profiles equally well. The result of fitting such a function provides an independent criterion to distinguish between the two types of profile, as well as an estimate of the transition point positions and their degree of abruptness. The function is a composite of terms chosen to satisfy defined requirements. We suggest that this kind of tailormade function created for specified purposes can be a useful tool for solving various analytical problems in growth research.
Theory
Requirements of the function
The desired function must demonstrate the following properties of velocity profiles:

At x = 0, y = 0. That is, velocity is defined to be zero at the point of reference. Note that this definition inverts the intuitive frame of reference. With respect to the plant's environment, the root tip has the maximal velocity and velocity falls to zero at a point that defines the end of the growth zone. In contrast, when the root tip is chosen as the reference point, velocity is zero at the tip and rises to a maximum at the point defining the end of growth zone. This inversion of the reference frame provides a host of mathematical advantages and is ubiquitously employed in kinematic analyses [1].

It must be able to describe a series of three intervals of linear relationship between x and y: the slopes of the first and second intervals will be positive with the second slope greater than the first, and the slope of the third interval will usually be zero. The transitions between linear domains must be continuous.

It also must be able to describe a sigmoidal relation between x and y. In other words, the function should be able to mimic conventional growth equations up to practical identity.
As we will see, a function that satisfies these requirements can be assembled using six coefficients. Two are needed to determine the nonzero slopes of the first two linear domains; we call these linear factors b_{1} and b_{2}, respectively, with the subscript indicating which of the two linear domains the coefficient controls. Two coefficients are required to determine the positions (c_{1} and c_{2}) of the two transitions, and the final pair of coefficients (d_{1} and d_{2}) determine the extent of the transitions. The extent of the transitions defines the linear versus sigmoidal character of the profile: the more extensive the transitions relative to the overall length of the profile, the more sigmoidlike the curve. In the following, we describe the function by breaking it down into its component building blocks.
Assembly of the stepstool function
The basic element of the function is an exponential term which can be used to create a smooth transition between two linear domains. Let us consider this term first in a form that resembles the socalled "expolinear" expression used in classical growth analysis [25, 26]:
where b, c, and d are constants and x is the variable. When c >> exp(b d x), then y reduces to the constant (ln c)/d, indicating that at small values of x, there is a linear domain that has zero slope. On the other hand, when c << exp(b d x), then y approaches (b x). Therefore, the function defines a transition between a line of zero slope at small x and one of slope b at large x. The transition is characterized by its location and width, which are determined by the coefficients c and d.
However, with equation (1), there is no direct correspondence between c or d and the location of the transition. Such a direct correspondence is desirable, though, because it would facilitate the process of initial value estimation, important for practical curvefitting. Looking again at equation (1), we note that the center of the transition zone is located at the x value at which c = exp(b d x). Thus, if we replace the coefficient c in equation (1) by the expression exp(b c_{1}d), the new coefficient c_{1} will have the value of x at the center of the transition. We now rewrite equation (1) using the conventions for coefficient identification defined above, and arrive at
A sample y_{I} with rather narrow transition zone centered on x = 0.3 is shown in Fig. 2.
Obviously, the assumption that the first linear domain has zero slope does not necessarily hold for real velocity profiles. To allow for nonzero slopes, a linear term, b_{1}x, has to be added to y_{I}:
where, for simplification, we define:
K = ln (exp [b_{2} c_{1} d_{1}] + exp [b_{2} d_{1} x]) (4)
In y_{II}, the slope of the first linear domain is b_{1}, while that of the second equals b_{1} + b_{2} (Fig. 2). As we will use it later in forcing the complete function to attain the value 0 at x = 0, we note that the value of y_{II} at x = 0 is
To create a second transition between the second linear domain and a third one with zero slope (which corresponds to the nongrowing parts of the root proximal of the growth zone), we employ the same tools. Taking another look at equation (2), we see that the variable x appears only in the term b_{2}d_{1}x in the second exponential expression, and that the first exponential expression contains the same term with x replaced by the constant, c_{1}. The simple function of x, f(x) = b_{2}x, in the second exponent can be replaced by any function of x, call it g(x), and the first exponent can be exchanged for g(x) with x replaced by a new constant, c_{2}. The result will always be a transition between a zeroslope linear domain and g(x). For example, we could insert y_{II} as g(x), and subtract the resulting term from y_{II} itself:
With c_{2} and d_{2} appropriately adjusted, we thus create a second transition zone right of the first one (Fig. 2; y_{III}). For x values significantly greater than x at the position of this second transition, this results in a third linear domain with zero slope.
At x = 0, equation 6 becomes:
where, for simplification, we define
L = ln (exp [b_{2} c_{1} d_{1}] + exp [b_{2} c_{2} d_{1}]) (8)
By definition, velocity equals zero at the point of reference (conventionally the root apex, i.e. x = 0 in our model), and we would like to see this feature in our function. Therefore we subtract y_{III,0} from y_{III}:
y_{IV} = y_{III}  y_{III,0} (9)
This function passes through the origin and asymptotically approaches a constant maximal value after the second transition zone (Fig. 2). In contrast to most classical growth equations, there is no parameter in our function which directly corresponds to this asymptote; it rather is determined indirectly by the geometry of the curve which results from the cooperation of all of the function's coefficients. The coefficients control the graph's shape in specific ways: b's determine the slopes of the linear domains, which are b_{1} (first linear domain) and b_{1} + b_{2} (second linear domain). The c's are the x values at the positions of the first and second transition zone (c_{1} and c_{2}, respectively), while d's define the width of the transition zones (the greater d, the narrower the transition zone). When equation 9 is used to describe root growth zones, negative values of any b's and c's will be meaningless; similarly, d's will be positive and c_{2} > c_{1} will always hold. When fitting the function to experimental data, it is advisable to implement these restrictions of possible coefficient values in the automated fitting algorithm to avoid unexpected results.
The complete function can be shifted along the yaxis by adding a constant q:
y_{V} = y_{IV} + q (10)
When fitting datasets which do not cover the region close to x = 0 (i.e. the root apex), this may be useful although it formally is a violation of the assumption that velocity is zero at x = 0.
Derivatives of the stepstool function
The derivative of y_{IV}, equalling that of y_{V}, is:
where
and
This derivative gives us the profile of elemental growth rate along the growth zone. The second derivative of y_{IV} also is helpful as it facilitates the identification of stationary points in the growth rate profile such as the position of its maximum:
Results
First, we explored the flexibility of the stepstool function, y_{IV}, and found that it could model any desired stage in the transition from bilinear to sigmoidal growth profiles (Fig. 3A–D). As illustrated in Fig. 1, one of the biologically relevant differences between sigmoidal and bilinear velocity profiles is the length of the period in which growing cells remain at their maximal growth rate. To quantify this parameter, we computed the timecourse of elemental growth rate and expressed the flatness, F, of the curve as the ratio of the periods in which the root elements grew at greater than 90% and at greater than 50% of their maximal rate. As expected, flatness decreased with increasing sigmoidal character of the profile (Fig. 3F; is shown to the right).
Our function also could produce asymmetric profiles (Fig. 3E, F). The ability to describe asymmetric sigmoidal profiles is the hallmark of the Richardsfunction, a classical growth equation with four coefficients [27], which added flexibility to the socalled functional approach of growth analysis [28]. Figure 3E, F illustrates that our function achieved the aim of mimicking flexible growth equations such as Richards'.
Having established a stepstool function with promising features, we tested its behavior in automated curve fitting procedures. To this end, we created MonteCarlo datasets [29]: various levels of Gaussian noise were added to several model velocity profiles, such as those shown in Figure 3, with a number of data points that was comparable to that of our empirical datasets. Then, the function was fitted using the MarquardtLevenberg algorithm included in the nonlinear regression tool of SigmaPlot, a scientific graphics and data analysis package widely distributed among biologists. Briefly, this algorithm searches for coefficient values for a bestfit of a given function as defined by the leastsumofsquares criterion. The search procedure is an iterative improvebyguessandtry process which does not establish the correct result, but rather provides an estimate of an acceptable solution. The significance of the estimated coefficients and their mutual dependences need to be assessed by statistical analyses of the fitting results, which are provided automatically by modern software packages such as the one used here. Particularly in the case of complex, versatile functions, the fitting algorithm is sensitive to the initial coefficient values from which it starts, and a smart guess of these values is an essential step [30]. As expected, the stepstool function was sensitive in this respect (Fig. 4). However, when reasonable coefficient values had been established by an initial round of manual curve fitting (see Additional File 1 for details), the results of the subsequent automated fitting procedure were reproducible, robust against small changes in coefficient initial values, and insensitive to Gaussian noise in the datasets, for all bilinear and sigmoidal velocity profiles examined (Fig. 5 gives two examples).
On this basis, we proceeded to test the function on real root growth data. We selected three datasets produced by the highresolution methodology recently introduced [12]; visual inspection suggested that these datasets possessed sigmoidallike (Gypsophila elegans; Fig. 6), intermediate (Aurinia saxatilis; Fig. 7), and strong bilinear (Arabidopsis thaliana; Fig. 8) characteristics. Fitting of the stepstool function invariably started with a manual adjustment of the coefficients to provide a "smart guess" of initial values for the subsequent automated fitting procedure. As the fitting of the stepstool function was intended to provide an independent criterion to decide whether a given dataset had more or less sigmoidal properties, we routinely included the following check for reliability of the fitting results. As the samples in Fig. 3A to 3D demonstrate, the degree of stepstoollikeness of the derivative of a fitted curve depends mostly on the coefficients d, which determine the abruptness of the transitions; the greater d_{1} and d_{2}, the more abrupt the transitions and the less sigmoidlike the curve. To see whether the degree of sigmoid likeness indicated by the automated fit was robust, we repeated the fitting procedure twice, starting with initial values of the coefficients d either doubled or halved. In all our examples (Figs. 6, 7, 8), the differences between the fitting results were insignificant, indicating that the curves obtained represented the best fit of the stepstool function to the datasets.
Inspection of the dataset obtained from a G. elegans root (n = 3678) revealed a pronounced sigmoidal character; fitting of the stepstool function supported this conclusion (Fig. 6A; the derivative of the function fitted is shown in Fig. 6B). Because the dataset showed no indication of a nonzero slope for the linear domain close to the root apex, a truncated version of the stepstool function (y_{IV}, equation 9) was used in which b_{1}, the coefficient controlling the slope of the first linear domain, was set to 0. Thus, the fitted function had only five coefficients instead of the six of the complete version. To test for the appropriateness of our decision, we also fitted the complete version of the function. The resulting curve was practically the same as was the sum of squared residuals, with the estimated b_{1} having a value close to 0 and an enormous variance, indicating that b_{1} had virtually no effect on the goodnessoffit (all other coefficients were significant and noncorrelated). Thus, the fitting statistics confirmed that there was no justification to include the sixth coefficient, b_{1}, in the fitted function for this data set.
To test the accuracy of the fit, we used DurbanWatson statistics (a test for correlation between residuals [31]) and the Levene median test for the homogeneity of variances of residuals [32], as implemented by the software package employed. The fit of the stepstool function failed in these tests, which could be explained by the presence of significant trends in the dataset that were not accounted for by the fitted function. To discover such trends, we examined a standardized residual plot (Fig. 6C). Residuals are the differences between the predicted and measured values; standardization divides them by the standard deviation of the particular set of residuals [33]. The plot showed that the residuals were distributed around the expected values in a nonstochastic manner; oscillatory deviations of various wavelengths from the predicted value were evident. This qualitative observation can be quantitatively assessed by a socalled runs analysis; here, a run is defined as a series of observations consisting of at least one measurement, in which all observations are either above or below the predicted values. For datasets comprising 3678 observations and fit to the correct model, residuals distributed stochastically would be expected to show 1840 runs (standard deviation = 30) [31]. However, the fit of the stepstool function in Fig. 6A produced only 181 runs, confirming the significance of the trends observed in the residual plot (Fig. 6C).
For comparison, we fitted the Richards function to the dataset (blue curve in Fig. 6A) and also plotted its derivative (Fig. 6B) and distribution of residuals (Fig. 6D). For this fit, the coefficient of determination was only slightly lower than for the stepstool function (r^{2} = 0.995 compared to r^{2} = 0.997). In contrast, the sum of squared residuals, the primary measure of goodnessoffit used by the fitting algorithm, was almost doubled (0.173 compared to 0.094). Intriguingly, the estimate of the location of the maximal elemental elongation rate differed notably between the Richards and the stepstool function, because the profile defined by the Richards fit was leftskewed (Fig. 6B). Comparison of the standardized residual plots (Fig. 6C, D,) shows that a longwavelength deviation from the measured values was more pronounced in the fitted Richards function, especially around the position of maximal growth rate. We conclude that the stepstool function, while not able to describe all minor trends present in this highresolution dataset, provides a more accurate description of the data than does the Richards function.
For this data set, we were surprised by the inferior fits from the Richards function, given that it has a pronounced sigmoidal character and therefore falls within the purview of classical, sigmoidal growth functions such as Richards'. For the stepstool function, the fitting result statistics indicated highly significant contributions to the prediction of the independent variable by all of the five coefficients of the truncated function used. This finding justified the inclusion of a fifth coefficient (as compared to the Richards function, which has four) on statistical grounds; we further verified the conclusion by comparing the stepstool and Richards models by the corrected Akaike's information criterion (AIC_{C}; [30]). This information theorybased criterion provides information on whether an improvement of goodnessoffit due to the inclusion of additional parameters is significant. This stringent test is preferable to the popular Ftest, particularly when comparing nonnested models, as done here [30] ([34] gives an indepth discussion of the informationtheoretical basis; for a more condensed introduction, see [35]). The AIC_{C} score (equation 17) of the stepstool model was smaller than that of the Richards function (38881 compared to 36640; ΔAIC_{C} = 2241), confirming that the fiveparameter stepstoolfit provided a better description of this particular dataset than the fourparameter Richards fit. To appreciate the robustness of this conclusion, note that a ΔAIC_{C} of 4.6 implies that the model with the lower score is 10 times more likely to be true than its competitor [[30]; equation [18]]: with a ΔAIC_{C} of 2241, the factor by which the stepstool model is more likely than the Richards function to be a correct description of this dataset is 10^{487}!
The second dataset, this time from an Au. saxatilis root (n = 3366), did not extend to position x = 0 (Fig. 7A), exemplifying a situation that frequently occurs because the collection of reliable data from the root meristem where velocities approach zero is difficult. Moreover, the data available seemed to include an initial, zeroslope linear domain. Therefore, we fitted a truncated version of the stepstool function lacking the coefficient b_{1} as in the previous example, but including the constant, q, to enable nonzero values of the function at position x = 0 (equation 10). The fitted curve described the data well (Fig. 7A), and the analysis of variances of the estimated coefficient values showed that all six coefficients, including q, made significant contributions to the prediction of the dependent variable.
The shape of the fitted curve resembled that seen in the previous example, but the graph of its derivative appeared wider and flatter (Fig. 7B). As a consequence, the flatness value was increased to almost 0.4, indicating that in this root, cells grew near their maximal growth rate for a relatively longer proportion of their phase of expansion. Again, there were nonstochastic components visible in the residual plot, showing up as longwavelength and shortwavelength oscillatory deviations (Fig. 7C). The fit of a fivecoefficient version of the Richards function (including a constant, q, as in the stepstool function; Fig. 7A) produced a coefficient of determination slightly lower than that of the stepstool function (0.997 compared to 0.999), but the sum of squared residuals was twice as high (0.269 compared to 0.136) and the nonstochastic oscillations in the residual plot were more pronounced (Fig. 7D). Thus, the stepstool function provided a more accurate description of the dataset than the classical sigmoidal model, which was confirmed by its lower AIC_{C} score (34038 as compared to 31745).
In the third example, an A. thaliana root, the bilinear character of the velocity profile was unambiguously visible in the raw data plot (Fig. 8A; n = 2729). As expected, the superiority of the description provided by fitting the stepstool function (the complete sixcoefficient version, equation 9) as compared to the fourparameter Richards function was obvious from the plots of the fitted velocity curves alone (Fig. 8A; coefficients of determination were 0.999 and 0.992, and sums of squared residuals were 0.084 and 0.509 for the stepstool and the Richards function, respectively). As expected, the AIC_{C} scores for the stepstool and the Richards model (28337 and 23424, respectively), as well as the tests of significance and independence, showed that the inclusion of two additional parameters in the stepstool fit was justified. The biphasic character was clearly expressed in the derivative of the stepstool function (Fig. 8B); the flatness value of 0.86 indicated that cells in this root expanded near their maximal growth rate for a substantial part of their elongation phase. Comparison of the residual plots (Fig. 8C, D) confirmed that the sigmoidal Richards function provided an inappropriate description of this dataset. However, as in the cases discussed before, the fit of the stepstool function failed to pass the DurbanWatson and Levene median tests, indicating nonstochastic factors in the distribution of residuals (Fig. 7C) and consequently, the existence of significant trends in the dataset unresolved by the stepstool fit.
Discussion
Recent technical advances have facilitated the demonstration of velocity profiles along root growth zones that consist of two distinct, nearly linear domains rather than being sigmoidal [12]. In contrast to the sigmoidal models, which dominate textbooks and previous research reports, bilinear growth profiles imply that cells switch between two distinct expansion modes when going through their period of elongation (Fig. 1). Because of this implication, the distinction of bilinear from sigmoidal growth zones is essential in kinematic growth studies. The stepstool function was created as a tool to characterize the shape of velocity profiles without imposing either a sigmoidal or a bilinear character. Because we currently do not possess a physiological hypothesis from which to derive a mechanistic, quantitative model of the regulation of velocity across a growth zone, the coefficients in the stepstool model refer to readily observed features of the profile, including slope, transition position and abruptness, rather than to underlying physiological processes. Biological meaning enters our analyses via the geometry of the fitted curve, which ultimately reflects the characteristics of the timecourses of physiological parameters experienced by cells traversing the growth zone (Figs. 1, 3).
Some velocity profiles may be truly sigmoidal while others may be bilinear. Insofar as the possible results of fitting our function include both alternative extremes, our approach is a multipleworkinghypotheses one, which conceptually complements conventional model testing based on statistical quantification of the goodnessoffit (as exemplified by the comparative evaluations of the stepstool and the Richards function; for a general discussion of multipleworkinghypotheses approaches in biological modelling, see [29]). In this context, the increased number of parameters in the stepstool function as compared to conventional growth models such as the Richards function should not be viewed as merely a means of improving the goodnessoffit. Rather, it is the cost of gaining an additional criterion for model choice by fitting an equation capable of describing two competing models equally well. Because there is no exclusive answer as to whether a given dataset is either sigmoidal or bilinear, a pragmatic measure of "bilinearity" has to be defined. As one possibility, we here introduced F, the flatness of the timecourse of elemental growth rate.
Our analyses of noisy MonteCarlodatasets and three exemplary experimental profiles demonstrated the success of our approach: the stepstool function provided better fits in all cases than did the Richards function; noteworthily, the latter is superior to other sigmoid functions due to its flexibility and therefore is considered a standard in growth analysis [27, 28, 36]. The flexible application of the stepstool function, which includes "smart guessing" of initial parameter values as well as the addition or deletion of coefficients depending on the properties of a particular dataset, certainly requires some familiarity with the mathematical basis of nonlinear curve fitting. However, this should not be a major obstacle, given the availability of powerful, userfriendly software and desktop computing power that was unthinkable of at the time the Richards function was introduced. In general, the flexibility of tailormade modular functions created to assist the solution of specific analytical problems is an asset in cases in which no mechanistic hypothesis has been developed, and from which explanatory mathematical models with biologically meaningful interpretations of parameters could be derived.
The application of the stepstool function may not be limited to kinematic growth analysis, as the necessity to distinguish between continously and abruptly changing curves is a notorious source of problems in classical plant growth analysis [24] as well as in the analysis of procaryotic [37–39] and eucaryotic [39–41] cell expansion. Formally similar problems exist in the field of biochemical [42, 43] and transport kinetics [44]. Moreover, we currently are applying extended, trilinear versions of the stepstool function in the description of polyphasic stomatal movements and the spatiotemporal quantification of volume changes in contractile forisomes [45].
With respect to the kinematic analyses discussed in this study, the finding of a velocity profile that is unambiguously a sigmoid implies that under certain conditions cells may regulate elemental elongation rate rather smoothly as they traverse the growth zone. This possibility merits further study; however, the purpose of the present work is to test the behavior of the stepstool function, and it will be up to subsequent investigation to determine whether a sigmoid profile for G. elegans indeed represents root growth in that species.
The application of the stepstool function will also form the basis for further investigations into the causes of the nonstochastic trends evident in the residual plots (subfigures D in Figures 6, 7, 8). It cannot be excluded that such deviations may occur as artifacts in the highresolution computational analysis applied to serial images of growing roots [12]. However, we consider it more likely that the oscillatory deviations reflect the behaviour of cells or cell groups that grow slightly faster or slower than the average growth curve suggests. This interpretation is substantiated by a previous report suggesting that local minima and maxima of growth rate travel along the growth zone similarly as cells do [46]. Future studies will aim to establish the nature of these local deviations of growth rate.
Conclusion
We present here an empirical function, with six coefficients, that is able to fit onedimensional velocity profiles regardless of whether they are bilinear or sigmoidal. Furthermore, the values of the coefficients provide analytical estimates of the key parameters of the profile, including the local slope, the position where the change of slope is centered, and a parameter characterizing the abruptness of the change in slope. The function is appropriate for data sets with thousands of points, as generated by modern, digital methods of velocity estimation, and fitting the function is robust to Gaussian noise. Finally, we demonstrate examples of velocity profiles from real roots that are sigmoidal, bilinear, and intermediate. We suggest that this function will facilitate kinematic analysis of growth and that our strategy for constructing the function may prove to be useful in general for quantitative biology.
Methods
Plant material and growth measurement
The plant material and growth measurement have been described in full previously [12]. Briefly, seedlings of Arabidopsis thaliana (L.) Heynh. Columbia background, Aurinia saxatilis (L.) Desv., and Gypsophila elegans Bieb. were grown on the surface of a nutrientagar medium in vertical Petri dishes under continuous light. A plate was put on the stage of a horizontal microscope and a stack of nine images collected of the primary root with 10 seconds between images. A series of four to eight stacks, depending on the root, was captured, from the root tip to an area with mature root hairs. The velocity field was recovered from these image sequences by RootflowRT, software that recovers dense velocity fields for deformable motion based on principles of optical flow. Velocities are obtained for essentially each pixel in the image and the component parallel to the local tangent of the root's midline is averaged perpendicular to the midline to produce the onedimensional velocity fields used here. Further details available in [12] and [47]. RootflowRT may be downloaded from [48].
Curve fitting
All mathematical operations including the generation of model velocity fields, creation of MonteCarlo datasets by adding Gaussian noise to these models, and numerical integration were performed using SigmaPlot (version 7.101, SPSS Inc., Chicago IL, USA) and TableCurve 2D (version 5.1, Systat Software Inc., Richmond CA, USA). Plotting of graphs and automated curve fitting was carried out with SigmaPlot; standard procedures implemented in this software were used to statistically analyze the fitting results. These tests included estimations of the standard errors, coefficients of variation, dependencies, tstatistics and P values of the parameters fitted, which were used to judge parameter significance and noncorrelatedness. Special attention was paid to the residuals, which were tested routinely for correlation by DurbinWatson statistics, normal distribution around the regression, and constant variance by the Levene median test.
The construction of the stepstoolfunction is described in detail in the Theory section. For comparison, the Richards function (Richards, 1959) was also fitted to experimentally determined velocity profiles (a, b, c, d are constants, x is the variable):
y = a (1  exp [b  cx])^{d} (15)
It should be noted that the graph of this function invariably intersects the ordinate at a(1  exp [b])^{d}; a fifth parameter has to be added if flexible intersections are required. The derivative of the Richards function is
To determine whether an improvement of the goodnessoffit due to the inclusion of additonal parameters in a model was significant, we calculated the corrected Akaike's information criterion (AIC_{C}) according to [30]:
where n is the number of datapoints, k is the number of parameters fitted plus one, and SS is the sum of squared residuals. The difference of the AIC_{C} scores (ΔAIC_{C}) of two competing models is a function of the relative likelihood of each of them to be a true description of the particular dataset to which they have been fitted; the model with the lower AIC_{C} is more likely to be true by a factor known as the evidence ratio:
References
 1.
Silk WK: Quantitative descriptions of development. Annu Rev Plant Physiol. 1984, 35: 479518.
 2.
Erickson RO: Modeling of plant growth. Annu Rev Plant Physiol. 1976, 27: 407434.
 3.
Erickson RO, Sax KB: Elemental growth rates of the primary roots of Zea mays. Proc Am Phil Soc. 1956, 100: 487498.
 4.
Goodwin RH, Avers CJ: Studies on roots. III. An analysis of root growth in Phleum pratense using photomicrographic records. Am J Bot. 1956, 43: 479487.
 5.
Hejnowicz Z: Growth and differentiation in the root of Phleum pratense. II. Distribution of cell divisions in the root. Acta Soc Bot Polon. 1956, 25: 615628. (in Polish)
 6.
Sharp RE, Silk WK, Hsiao TC: Growth of the maize primary root at low water potentials. Plant Physiol. 1988, 87: 5057.
 7.
Beemster GTS, Baskin TI: Analysis of cell division and elongation underlying the developmental acceleration of root growth in Arabidopsis thaliana. Plant Physiol. 1998, 116: 15151526.
 8.
Peters WS, Felle HH: The correlation of profiles of surface pH and elongation growth in maize roots. Plant Physiol. 1999, 121: 905912.
 9.
Schnyder H, Nelson CJ, Coutts JH: Assessment of spatial distribution of growth in the elongation zone of grass leaf blades. Plant Physiol. 1987, 85: 290293.
 10.
BenHajSalah H, Tardieu F: Temperature affects expansion rate of maize leaves without change in spatial distribution of cell length. Plant Physiol. 1995, 109: 861870.
 11.
Fricke W, Peters WS: The biophysics of leaf cell growth in saltstressed barley: a study at the cell level. Plant Physiol. 2002, 129: 374388.
 12.
van der Weele CM, Jiang HS, Palaniappan KK, Ivanov VB, Palaniappan K, Baskin TI: A new algorithm for computational image analysis of deformable motion at high spatial and temporal resolution applied to root growth. Roughly uniform elongation in the meristem and also, after an abrupt acceleration, in the elongation zone. Plant Physiol. 2003, 132: 11381148.
 13.
Silk WK, Lord EM, Eckard KJ: Growth patterns inferred from anatomical records. Plant Physiol. 1989, 90: 708713.
 14.
Brumfield RT: Cell growth and division in living root meristems. Am J Bot. 1942, 29: 533543.
 15.
Paolillo DJ, Sorrells ME, Keyes GJ: Gibberellic acid sensitivity determines the length of the extension zone in wheat leaves. Ann Bot. 1991, 67: 479485.
 16.
Ivanov VB, Maximov VN: The change in the relative rate of cell elongation along the root meristem and the apical region of the elongation zone. Russ J Plant Physiol. 1999, 46: 7382.
 17.
Silk WK, Erickson RO: Kinematics of plant growth. J Theor Biol. 1979, 76: 481501.
 18.
Gandar PW: Growth in root apices. Bot Gaz. 1983, 144: 119.
 19.
Peters WS, Bernstein N: The determination of relative elemental growth rate profiles from segmental growth rates. Plant Physiol. 1997, 113: 13951404.
 20.
Peters WS: Growth and extracellular pH in roots: How to control an explosion. New Phytol. 2004, 162: 571574.
 21.
Hunt R: Plant Growth Curves. 1982, London: Edward Arnold.
 22.
Barlow PW, Brain P, Parker JS: Cellular growth in roots of a gibberellindeficient mutant of tomato (Lycopersicon esculentum Mill.) and its wildtype. J Exp Bot. 1991, 42: 339351.
 23.
Morris AK, Silk WK: Use of a flexible logistic function to describe axial growth of plants. Bull Math Biol. 1992, 54: 10691081.
 24.
Fisher PR, Heins RD, Lieth JH: Quantifying the relationship between phases of stem elongation and flower initiation in Poinsettia. J Amer Soc Hort Sci. 1996, 121: 686693.
 25.
Goudriaan J, Monteith JL: A mathematical function for crop growth based on light interception and leaf area expansion. Ann Bot. 1990, 66: 695701.
 26.
Yin X, Goudriaan J, Lantinga EA, Vos J, Spiertz HJ: A flexible sigmoid function of determinate growth. Ann Bot. 2003, 91: 361371.
 27.
Richards FJ: A flexible growth function for empirical use. J Exp Bot. 1959, 10: 290300.
 28.
Causton DR, Venus JC: The Biometry of Plant Growth. 1981, London: Edward Arnold.
 29.
Hilborn R, Mangel M: The Ecological Detective: Confronting Models with Data. 1997, Princeton: Princeton University Press,
 30.
Motulsky H, Christopoulos A: Fitting Models to biological Data using linear and nonlinear Regression. 2004, Oxford: Oxford University Press.
 31.
Draper NR, Smith H: Applied Regression Analysis. 1998, New York: John Wiley & Sons, 3.
 32.
Conover WJ, Johnson ME, Johnson MM: A comparative study of tests for homogeneity of variances, with applications to the outer continental shelf bidding data. Technometrics. 1981, 23: 351361.
 33.
Glantz SA, Slinker BK: Primer of applied Regression and Analysis of Variance. 2001, NewYork: McGrawHill, 2.
 34.
Burnham KP, Anderson DR: Model Selection and multimodel Inference – A practical informationtheoretic Approach. 2002, New York: Springer, 2.
 35.
Myung JI, Pitt MA: Model comparison methods. Meth Enzymol. 2004, 383: 351366.
 36.
Richards FJ: The quantitative analysis of growth. Plant Physiology: A Treatise. Edited by: Steward FC. 1969, 5A: 376. New York: Academic Press.
 37.
Kubitschek HE: Bilinear cell growth of Escherichia coli. J Bacteriol. 1981, 148: 730733.
 38.
Cooper S: What is the bacterial growth law during the division cycle?. J Bacteriol. 1988, 170: 50015005.
 39.
Mitchison JM: Growth during the cell cycle. Int Rev Cytol. 2003, 226: 165258.
 40.
Cooper S: Length extension in growing yeast: is growth exponential? – Yes. Microbiology. 1998, 144: 263265.
 41.
Mitchison JM, Sveiczer A, Novak B: Length extension in growing yeast: is growth exponential? – No. Microbiology. 1998, 144: 265266.
 42.
CornishBowden A: Abrupt transitions in kinetic plots: an artifact of plotting procedures. Biochem J. 1988, 250: 309310.
 43.
Engel PC, Syed SEH: Abrupt transitions in kinetic plots: an empirical reality. Biochem J. 1988, 250: 310311.
 44.
Nissen P: Multiphasic uptake mechanisms in plants. Intern Rev Cytol. 1991, 126: 89134.
 45.
Knoblauch M, Peters WS: Forisomes, a novel type of Ca^{2+}dependent protein motor. Cell Motil Cytoskel. 2004, 58: 137142.
 46.
Salamon P, List A, Grenetz PS: Mathematical analysis of plant growth. Zea mays primary roots. Plant Physiol. 1973, 51: 635640.
 47.
Palaniappan K, Jiang H, Baskin TI: Nonrigid motion estimation using the robust tensor method. IEEE Computer Vision & Pattern Recognition Workshop on Articulated and Nonrigid Motion. 2004, 2533. Washington DC: IEEE Computer Society.
 48.
http://www.bio.umass.edu/biology/baskin/research_rootflow.htm
Acknowledgements
This work was supported by the U.S. National Science Foundation (award no. IBN 0316876 to T.I.B.). W.S.P. thanks Ch. RothKäppchen for clarifying discussion of conceptual issues.
Author information
Additional information
Competing interests
The author(s) declare that they have no competing interests.
Authors' contributions
T.I.B. contributed the experimental growth data, W.S.P. formulated the stepstool equation and performed the statistical analyses. Both authors conceived of this study while discussing their previous work and produced the manuscript cooperatively.
Authors’ original submitted files for images
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Velocity Profile
 Growth Zone
 Linear Domain
 Richards Model
 AICC Score