A General Model for Describing the Ovate Leaf Shape

Many plant species produce ovate leaves, but there is no general parametric model for describing this shape. Here, we used two empirical nonlinear equations, the beta and Lobry–Rosso– Flandrois (LRF) equations, and their modified forms (referred to as the Mbeta and MLRF equations for convenience), to generate bilaterally symmetrical curves along the x-axis to form ovate leaf shapes. In order to evaluate which of these four equations best describes the ovate leaf shape, we used 14 leaves from 7 Neocinnamomum species (Lauraceae) and 72 leaves from Chimonanthus praecox (Calycanthaceae). Using the AIC and adjusted root mean square error to compare the fitted results, the modified equations fitted the leaf shapes better than the unmodified equations. However, the MLRF equation provided the best overall fit. As the parameters of the MLRF equation represent leaf length, maximum leaf width, and the distance from leaf apex to the point associated with the maximum leaf width along the leaf length axis, these findings are potentially valuable for studying the influence of environmental factors on leaf shape, differences in leaf shape among closely related plant species with ovate leaf shapes, and the extent to which leaves are bilaterally symmetrical. This is the first work in which temperature-dependent developmental equations to describe the ovate leaf shape have been employed, as previous studies lacked similar leaf shape models. In addition, prior work seldom attempted to describe real ovate leaf shapes. Our work bridges the gap between theoretical leaf shape models and empirical leaf shape indices that cannot predict leaf shape profiles.


Introduction
The ovate leaf shape is a common geometric form among many plant species, especially in the Lauraceae. It therefore exists in almost all floras. Similar linear leaf forms with comparable structural characteristics occur in other widely spread flowering plant families. However, there is no strict distinction between these two leaf shapes, although the difference between these two kinds of shapes comes from the variation in the ratio of leaf width to leaf length [1][2][3]. Previous studies show that a simplified version of a polar coordinate equation (proposed by Gielis [4]) can describe all shapes of bamboo leaves, including narrow linear leaves and even some with approximate ovate shapes [5,6].
However, variations in the general ovate shape occur as a result of subtle changes in the curvature along the leaf margin from the leaf apex to the point of the maximum leaf width. These variations limit the utility of the simplified Gielis equation in generically describing the ovate leaf shape [1,7].
In order to address these biological variations in leaf form mathematically, we employed the beta equation [8] and the Lobry-Rosso-Flandrois (LRF) equation [9,10]. These equations have two significant advantages: (i) they provide a better goodness of fit than other equations, and (ii) their four-parameter space reflects a better close-to-linear behavior for a nonlinear fit [11,12]. These two equations can generate curves that look like one side of an ovate leaf (as defined by the midrib running along the length of the leaf lamina). Importantly, the four parameters in each of these two equations can be shown to correspond to half of the maximum leaf width, the point on the leaf length axis associated with the leaf maximum width, the leaf base that connects the lamina and the petiole, and the leaf apex. These parameters, especially the first two, are especially meaningful for describing the shapes of ovate leaves. Specifically, the ratio of the distance from leaf apex to the point on the leaf length axis associated with leaf maximum width and total leaf length can be used to study the evolution and variation in leaf shape, even among distantly related species [13]. In addition, Shi et al. [14] proposed more general forms for the beta and LRF equations by introducing an additional parameter for each equation, which makes both equations more flexible and morphologically versatile.
Although some parametric leaf shape models exist, these models focus on linear leaves, especially those of crops and bamboos [5,6,15]. Here, we further modify the beta and LRF equations and apply their modified forms to generate ovate leaf shapes and test their utility by examining the leaves of seven species (including a variety) of the genus Neocinnamomum (Lauraceae). These species have typical ovate leaves that are morphologically representative of the general ovate leaf form of many other species (e.g., Prunis). In addition, to test whether an additional parameter is valuable in the modified beta and modified LRF equations, we used 72 individual leaves from the species Chimonanthus praecox (Calycanthaceae) to examine the goodness of fit and its trade-off with model complexity. Previously, it was shown that similar models were useful for fitting temperature-dependent developmental rates of arthropods [12] and the ontogenetic growth of plants [8,14]. If our four models could be demonstrated to be valid in describing the actual ovate leaf shapes of a plant group, it would show that these models have a wider scope in applications, and further suggest that there probably exists a similar mechanism for ontogenetic development and growth.
In this study, we aim to build several new leaf shape models that can describe ovate leaves and to find whether there is a best one relative to other candidate models. This work does not recommend a leaf shape index such as the leaf roundness index and leaf dissection index [16][17][18][19][20], because a leaf shape index cannot be used to fit the edge data of ovate leaves. Here, the mathematical model that we finally recommend can describe the ovate leaf shape with model parameters that can be well understood. In addition, the model parameters can be used to quantify leaf shape (see the Materials and Methods section for details). Although leaf shape can be affected by environmental factors [21,22], we did not attempt to check the influences of environmental factors on leaf shape, which exceeded the scope of the current study, but deserves future study.

Models
The original beta and Lobry-Rosso-Flandrois (LRF) equations, which have been used to describe the influence of temperature on the growth rate of crops and bacteria and the developmental rates of insects, produce similar left-skewed curves in the range between the lower and upper developmental thresholds (Figure 1a). Each has four parameters with the same or very similar geometric meanings [8][9][10]. Specifically, the and the LRF equation is: where y represents the developmental rate at temperature x; y c represents the maximum developmental rate; x c represents the temperature associated with the maximum developmental rate y c ; x 1 and x 2 represent the lower and upper developmental thresholds that are the left and right intersections of the curve and the x-axis (Figure 1a). The beta equation can produce left-or right-skewed curves, whereas the LRF equation produces only left-skewed curves. When using the leaf base as the starting point and the leaf apex as the ceiling point, another form of the LRF equation should be used [14], i.e., Equation (3) is actually symmetrical with Equation (2) around the vertical axis x = (x 1 + x 2 )/2 ( Figure 1b). To increase the flexibility of curve fitting, Shi et al. [14] proposed generalized versions of the beta and LRF equations, i.e., Modified Beta equation: Modified LRF equation: where δ is an additional parameter, the numerical value of which influences the curvature of the curve. To reduce the model's complexity in the current study, x 1 was fixed to be zero. An ovate leaf shape can be generated by defining a symmetrical equation along the x-axis for Equations (4) and (5): f (x) = y(x) and g(x) = −y(x), where f (x) and g(x) represent the upper and lower parts of the leaf shape curve, respectively. Figure 1c shows simulated leaf shape curves with different δ values. Hereafter, for convenience, we refer to Equations (4) and (5) as the Mbeta and MLRF equations. For each of the four models, the parameters can be used to quantify leaf shape characteristics. For instance, the ratio of x c to x 2 can reflect whether the maximum leaf width is more approximate to the middle of the lamina, and the ratio of 2y c to x 2 approximates the ratio of lamina width to length. For each of the four models, the parameters can be used to quantify leaf shape characteristics. For instance, the ratio of xc to x2 can reflect whether the maximum leaf width is more approximate to the middle of the lamina, and the ratio of 2yc to x2 approximates the ratio of lamina width to length.

Parameter Estimation for the Models
The coordinates of the leaf edge can be obtained by scanning leaf profile images [1]. The coordinates of the leaf base are not exactly located at (0, 0), and the angle between the leaf length axis (i.e., a straight line through leaf base and leaf apex) and the x-axis is usually not equal to 0 ( Figure 1d). Thus, to obtain the estimates of the model parameters, three additional parameters are required, i.e., x0, y0, and θ, where (x0, y0) represent the coordinates of the leaf base, and θ is the angle between the aforementioned two axes. Thus, for the beta and LRF equations, there are six parameters in total (i.e., three location parameters plus three model parameters); for the Mbeta and MLRF equations, there are seven parameters in total (i.e., three location parameters plus four model parameters). The residual sum of squares (RSS) of the y coordinate was used to measure the goodness of fit:

Parameter Estimation for the Models
The coordinates of the leaf edge can be obtained by scanning leaf profile images [1]. The coordinates of the leaf base are not exactly located at (0, 0), and the angle between the leaf length axis (i.e., a straight line through leaf base and leaf apex) and the x-axis is usually not equal to 0 ( Figure 1d). Thus, to obtain the estimates of the model parameters, three additional parameters are required, i.e., x 0 , y 0 , and θ, where (x 0 , y 0 ) represent the coordinates of the leaf base, and θ is the angle between the aforementioned two axes. Thus, for the beta and LRF equations, there are six parameters in total (i.e., three location parameters plus three model parameters); for the Mbeta and MLRF equations, there are seven parameters in total (i.e., three location parameters plus four model parameters). The residual sum of squares (RSS) of the y coordinate was used to measure the goodness of fit: where i denotes the ith point on the leaf edge; n denotes the number of data points on the leaf edge; and y i with a circumflex denotes the predicted value of y i using a nonlinear To reflect the trade-off between the goodness of fit and the model's complexity, the Akaike information criterion (AIC) was also calculated [23]: where p is the number of model parameters plus an error item. The equation with the lowest AIC is considered as the best one.
To diminish the influence of leaf size on the RSS, we used the adjusted root mean square error (RMSE adj ) (see References [24,25]) to intuitively compare the goodness of fit for different leaves: where A denotes the leaf area. The adjusted RMSE accounts for the proportion of the average absolute deviation in the y coordinate to the radius of an assumed circle whose area is equal to the leaf area, and this is used as a reference indicator during the comparison of models.

Samples of Representative Ovate Leaves
Representative plants of Neocinnamomum (Lauraceae) species with representative ovate leaf shapes were used to determine which among the four candidate equations provided the best fit. The leaves were collected from southern China, where the plants grew naturally (Table 1). Two healthy and mature representative leaves from each species were collected ( Figure 2). To further compare the four candidate equations, we used 72 leaves of Chimonanthus praecox (Calycanthaceae) growing in the Nanjing Forestry University campus collected on 20 October 2017. There is no evidence to show that the geographical location, plant age, leaf position in the canopy, or light conditions can significantly alter the basic leaf shape characteristics of ovate leaves. Thus, we randomly sampled leaves from the plants that we had access to. Under any circumstances, any significant change in leaf shape would have no effect on evaluating which equation best described leaf shape. Indeed that is the whole point of this study.

Leaf Image Processing and Data Acquisition
The Neocinnamomum leaves were photographed in jpg format, and the leaves of C. praecox were scanned using a photo scanner (HP Scanjet 4850; Hewlett-Packard Company, Palo Alto, California) at 400 dpi resolution. A Matlab procedure (version ≥ 2009a) [1] was used to extract planar coordinates for each leaf. The numbers of data points on the leaf edge for the 14 Neocinnamomum specimens ranged from 1000 to 2700 depending on the leaf size; the data points of the leaves of C. praecox were set to be ca. 3000. These numbers of data points are sufficient to form an integrated leaf shape. Each leaf image was transferred to a black-white image in bmp format. We fitted the four candidate equations (i.e., the beta, LRF, Mbeta and MLRF equations) to the leaf shape data using the Nelder-Mead optimization method [26], which was carried out using the "optim" function in R software (version 3.6.1) [27]. The target function is the RSS (see Equation (6) for details).

Leaf Image Processing and Data Acquisition
The Neocinnamomum leaves were photographed in jpg format, and the leaves of C. praecox were scanned using a photo scanner (HP Scanjet 4850; Hewlett-Packard Company, Palo Alto, California) at 400 dpi resolution. A Matlab procedure (version ≥ 2009a) [1] was used to extract planar coordinates for each leaf. The numbers of data points on the leaf edge for the 14 Neocinnamomum specimens ranged from 1000 to 2700 depending on the leaf size; the data points of the leaves of C. praecox were set to be ca. 3000. These numbers of data points are sufficient to form an integrated leaf shape. Each leaf image was transferred to a black-white image in bmp format. We fitted the four candidate equations (i.e., the beta, LRF, Mbeta and MLRF equations) to the leaf shape data using the Nelder-Mead optimization method [26], which was carried out using the "optim" function in R software (version 3.6.1) [27]. The target function is the RSS (see Equation (6) for details).

Results
For the 14 Neocinnamomum leaf specimens, the Mbeta and MLRF equations provided better fits than the beta or LRF equations. Specifically, 11 out of 14 AIC and RMSEadj values of the MLRF equation were smaller than those using the Mbeta equation (Table 2). This result was interpreted to indicate that the MLRF equation is the best model among the other candidate equations. The adjusted RMSE values using the MLRF equation for the 14 leaves ranged from 0.0163 to 0.0424, i.e., the average absolute deviation of the curve fitting was less than 5% of the radius of a circle, assuming that its area is equal to the leaf area. This result further validated the superiority of the MLRF equation for describing ovate leaf shapes. Figure 3 shows the fitted results using the MLRF equation.

Results
For the 14 Neocinnamomum leaf specimens, the Mbeta and MLRF equations provided better fits than the beta or LRF equations. Specifically, 11 out of 14 AIC and RMSE adj values of the MLRF equation were smaller than those using the Mbeta equation (Table 2). This result was interpreted to indicate that the MLRF equation is the best model among the other candidate equations. The adjusted RMSE values using the MLRF equation for the 14 leaves ranged from 0.0163 to 0.0424, i.e., the average absolute deviation of the curve fitting was less than 5% of the radius of a circle, assuming that its area is equal to the leaf area. This result further validated the superiority of the MLRF equation for describing ovate leaf shapes. Figure 3 shows the fitted results using the MLRF equation.
The MLRF equation was also found to be the best one among the four candidate equations based on the 72 individual leaves collected from C. praecox (Figure 4a). Specifically, the AIC values of the Mbeta equation were greater than those of the MLRF equation for 67 out of 72 leaves. In addition, the adjusted RMSE values of the Mbeta equation were greater than those of the MLRF equation for the same number of leaves (Figure 4b). The adjusted RMSE values using the MLRF equation for the 72 leaves ranged from 0.0087 to 0.0385, indicating that the average absolute deviation of the curve fitting was less than 4% of the radius of a circle whose area is hypothesized to equal the leaf area. The MLRF equation fitted the leaf edge data of C. praecox well (see Figure 5 for example). adjusted RMSE values using the MLRF equation for the 72 leaves ranged from 0.0087 to 0.0385, indicating that the average absolute deviation of the curve fitting was less than 4% of the radius of a circle whose area is hypothesized to equal the leaf area. The MLRF equation fitted the leaf edge data of C. praecox well (see Figure 5 for example). Tables S1 and S2 in the online Supplementary Materials tabulate the estimates of the model parameters and the goodness of fit for the 86 leaves (14 + 72).  Figure 3. Comparisons between the scanned leaf edges (gray curves) and the predicted leaf edges (red curves) using the MLRF equation for the seven Neocinnamomum species. For each species, there are two replicates, with the first number of the code representing the species, and the second the replicate. For the species code, "1" represents N. lecomtei; "2" N. complanifructum; "3" N. fargesii; "4" N. delavayi; "5" N. caudatum; "6" N. caudatum var. macrocarpum; "7" N. mekongense. Figure 3. Comparisons between the scanned leaf edges (gray curves) and the predicted leaf edges (red curves) using the MLRF equation for the seven Neocinnamomum species. For each species, there are two replicates, with the first number of the code representing the species, and the second the replicate. For the species code, "1" represents N. lecomtei; "2" N. complanifructum; "3" N. fargesii; "4" N. delavayi; "5" N. caudatum; "6" N. caudatum var. macrocarpum; "7" N. mekongense.

Discussion
Dornbusch et al. [15] proposed an empirical model to describe the leaf shapes of several crop species (i.e., Triticum aestivum, Hordeum vulgare and Zea mays cultivars). The advantage of this model is that it only requires information pertaining to leaf length, leaf  Tables S1 and S2 in the online Supplementary Materials tabulate the estimates of the model parameters and the goodness of fit for the 86 leaves (14 + 72).

Discussion
Dornbusch et al. [15] proposed an empirical model to describe the leaf shapes of several crop species (i.e., Triticum aestivum, Hordeum vulgare and Zea mays cultivars). The advantage of this model is that it only requires information pertaining to leaf length, leaf maximum width, and the ratio of the distance from the leaf apex to the point on the leaf length axis associated with the leaf maximum width. However, this model consists of two step functions (i.e., from the leaf apex to the leaf maximum width and from the leaf maximum width to leaf base), which does not necessarily lead to smoothness at the junction between the two step functions (i.e., the point on the leaf edge curve at which the leaf has its maximum width). Shi et al. [5] proposed a simplified version of generating leaf shapes based on the original Gielis equation [4] using three location parameters and two model parameters. The modified Gielis equation can produce a similar but smoother curve than that of Dornbusch et al. [15], and it has fewer parameters (only two model parameters). Lin et al. [6] tested the validity of the simplified Gielis equation using the leaves of 42 bamboo species, and found that the equation fit them well. However, the simplified Gielis equation cannot produce a concave curve close to the leaf apex, the same limitation as that of the equation of Dornbusch et al. [15].
In contrast, the curve produced by the MLRF equation presented here not only produces an ovate leaf shape similar to that produced by the simplified Gielis equation (see Figure 1c for the case of δ = 0.5), but it also produces a concave curve close to the leaf apex. If an ovate leaf is not strictly bilaterally symmetrical for the parts close to the leaf apex (e.g., when one side is a linear or convex curve, whereas the other side is a concave curve), more parameters can be set for one or the other side (e.g., two δ values can be fitted, one for each side of the leaf). Thus, the MLRF equation can describe a broad variety of ovate leaf shapes. For some linear leaf shapes that have the same structural characteristics as ovate leaf shapes but a lower y c value in the MLRF equation, the MLRF equation is still applicable, because there is no strict distinction between an ovate leaf shape and a linear leaf shape of the other kind. Schrader et al. [3] distinguished between the two leaf forms using the ratio of leaf width to leaf length, such that leaf shapes with ≤1/10 width/length ratio are defined as linear rather than ovate. Although this is an empirical definition, it shows that an ovate leaf shape can be defined as a linear leaf shape with a lower ratio of leaf width to length. Clearly, there are other linear leaf shapes that have structural characteristics that differ from the ovate leaf shape, i.e., linear leaf shapes appear to be similar because the small ratio of leaf width to length has hidden their different structural characteristics. However, these linear forms probably result from other leaf shapes that have been described as acicular, acuminate, cuneate, deltoid, or elliptical shapes. Examining these forms using the approach described here merits further investigation.
An important limitation of the MLRF equation is that its ability to mimic ovate shapes manifests errors at the leaf apex that can be large for some leaves (Figure 3). The actual leaf apex is longer and more skewed compared to the predicted form. This likely results from the evolution of what is called a drip-point in some Lauraceae species that can quickly drain water from the leaf surface (for a discussion of this phenomenon, see ref. [28]). We believe that this process-driven mechanism may be similar to temperature-dependent growth dynamics [29][30][31]. The process of leaf size growth has been demonstrated to follow a sigmoid curve [32]. According to References [8,14], a sigmoid curve results from integrating the beta equation or the LRF equation. Thus, the growth rate of a leaf results in curves that can still be simulated by the MLRF equation.
There are many statistical approaches, especially non-parametric models, such as elliptic Fourier, kernel estimation, local regression, neural network, and spline [33,34], which can depict leaf shape very accurately. However, these approaches have significant disadvantages such as requiring too many parameters in each step of multi-step calculations [33], and biologically important shape parameters fail to appear in these non-parametric approaches. In comparison with these models, the Mbeta and MLRF equations provide clear geometric and biological analogs for the curves they generate. For example, x c represents the point on the leaf length axis associated with half of the leaf maximum width (i.e., y c ), x 2 represents leaf length, and δ can be used to adjust the concavity and convexity of a curve. There are other advantages of using the Mbeta and MLRF equations over other nonlinear temperature-dependent developmental rate equations such as the Briére, Lactin and square root equations [35][36][37][38]. In these models, with the exception of the lower and upper thresholds (i.e., x 1 and x 2 ), there are other parameters that cannot directly provide a link to a curve's shape. However, the parameters in the MLRF equation can be intuitively understood and used to quantify the complexity of leaf shape. It can be used to assess leaf shape differences among different individuals of the same species or among different species that are taxonomically closely or distantly related. In addition, the MLRF equation, which generates bilaterally symmetrical leaf shapes perfectly, can be used to quantify the extent of bilateral asymmetry.

Conclusions
The modified beta and LRF equations (i.e., the Mbeta and MLRF equations) are demonstrated to be more flexible for data fitting than their original versions (the beta and LRF equations), and the MLRF equation has a lower AIC or adjusted RMSE value than the Mbeta equation for most of the data sets studied here. The MLRF equation can generate two types of ovate leaf shapes: for the part of the curve from the leaf apex to a point on the leaf edge not reaching the leaf maximum width (wherein one type is approximately linear and somewhat convex), and the other is concave. The MLRF appears to be more useful than other available leaf morphometric models for describing ovate leaf shapes. For example, the ratios of the average absolute deviation in the y coordinate to the radius of the assumed circle whose area equals leaf area are smaller than 5% for all leaves investigated. In addition, with the exception of δ, the other parameters of the MLRF can be intuitively understood in terms of geometric and biological features. Thus, the work presented has potential value for future studies understanding the influence of environmental factors on leaf shape for plants with ovate leaf shapes.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/sym13081524/s1, Table S1: fitted results of the four nonlinear equations to the leaf edge data of seven Neocinnamomum species, Table S2: