Towards More Realistic Leaf Shapes in Functional-Structural Plant Models

Fluctuating asymmetry in plant leaves is a widely used measure in geometric morphometrics for assessing random deviations from perfect symmetry. In this study, we considered the concept of fluctuating asymmetry to improve the prototype leaf shape of the functional-structural plant model L-Cucumber. The overall objective was to provide a realistic geometric representation of the leaves for the light sensitive plant reactions in the virtual plant model. Based on three-dimensional data from several hundred in situ digitized cucumber leaves comparisons of model leaves and measurements were conducted. Robust Bayesian comparison of groups was used to assess statistical differences between leaf halves while respecting fluctuating asymmetries. Results indicated almost no directional asymmetry in leaves comparing different distances from the prototype while detecting systematic deviations shared by both halves. This information was successfully included in an improved leaf prototype and implemented in the virtual plant model L-Cucumber. Comparative virtual plant simulations revealed a slight improvement in plant internode development against experimental data using the novel leaf shape. Further studies can now focus on analyses of stress on the 3D-deformation of the leaf and the development of a dynamic leaf shape model.


Introduction
The field of plant modeling provides insights to plant responses on changes in environmental parameters where experimental evaluation becomes infeasible.Simulation studies can include effects related to climate change, such as the increase in temperature and atmospheric CO 2 concentration, or focus on changes in controlled greenhouse conditions.While in the past two kinds of plant modeling approaches-architectural models and process-based models-were analyzed separately, they were combined to functional-structural plant models (FSPM, aka virtual plant models) more recently [1][2][3].Using process-based models, whole-plant reactions on environmental conditions are studied, while architectural models focus on establishing algorithms to simulate plant morphological development often described by Lindenmayer-Systems (L-Systems) [4].FSPMs allow studying both, the physiological and structural responses to global and local environmental conditions.Most prominent environmental factors include temperature, nutrient supply and light, e.g., photosynthetically active radiation (PAR) [1,2,[5][6][7].
Functional-structural plant models have already been developed for various crops, for instance for arabidopsis, wheat, cotton or sorghum, as well as for some trees, including fruit trees, such as peaches and apples.[1,2,[8][9][10][11][12][13].In vegetable sciences, FPSMs are limited to a few model organisms.Cucumber is one of those model crops and recently became the first genome-sequenced vegetable [14], as it is not only of interest for plant physiologist, but also implies economical interests for the food industry.Hence, the FSPM for the greenhouse cucumber Cucumis sativus L. is under steady development [5,6,15,16].
FSPM models usually rely on simplified representations of plant organs such as leaves and fruits.For instance, the current three-dimensional model leaf used in the cucumber FPSM (L-Cucumber [6]) consists of four equal-sized equilateral triangles, two per each leaf half, and has a mirror-symmetrical shape.While extensive research on two-dimensional leaf shapes exists, three-dimensional analysis, using e.g., 3D laser scanners or a multi-camera setup, focus mostly on entire plants or characteristics for automated phenotyping and has its limits, when plant organs overlap or are hidden inside the canopy [17][18][19][20][21].To overcome this limitation, digitalization with magnetic field digitizers, where three-dimensional coordinates of specified landmarks are captured using a manual pointer offer a labor intensive alternative.Previous studies on cucumber plant morphology utilizing data acquired with this approach were focused on major plant growth parameters, such as internode development, leaf areas or phyllotaxis angles to be included in a virtual plant model [22,23].Landmark-based morphometric methods and related measures, such as fluctuating asymmetry [24-27], were not yet under consideration in this field of research, while they are common in two-dimensional leaf shape analysis.Resulting in size free shape description they allow studies on leaf shape variation under different environmental conditions or among species [28][29][30][31][32].
The primary aim of this study was to identify a realistic leaf shape model from three-dimensional landmark data and compare it to the original leaf prototype of a cucumber virtual plant model.Most recent studies on virtual plant models focused on Arrhenius type models for plant responses to temperature variation and identified a model for internode development based on canopy light quantity and quality [5,6].Based on the findings on local light fluctuations affecting plant development, this study further aims at estimating the effect of the virtual plant model's leaf shape on the simulation results, as the leaf shape clearly affects the canopy-light-interaction [33,34].Furthermore, more realistic 3D plant canopies were found to be of importance in studies related to remote sensing, where canopy reflectance or other radiometric measures affect the remote sensing quality and interpretability [35].

Virtual Plant Model
The dynamic light-sensitive functional-structural plant model L-Cucumber from Kahlen and Chen [6] with an option for temperature-sensitivity ("MA-T"), hereinafter referred to as L-Cucumber, was used as the reference model for this study.It was set up for virtual plant simulations of greenhouse grown cucumber plants and considers ambient temperature and local light conditions in the production of dry matter and growth of the shoot.The coupled light model provides local data from Quasi-Monte Carlo simulation [36] on photosynthetically active radiation (PAR) and red-to-far-red ratio (R:FR) for leaves and internodes.This information is incorporated in model equations regulating plant development of L-Cucumber [5,6], e.g., the growth of an internode, ER(t) at day t (cm), is modeled following where F av is a temperature response function in the same period as the PAR signal and IL(t) is the internode length at day t [5,6].The final internode length (FIL, (cm)) is a function of both PAR and R:FR [5]: Hence, IL(t) (cm) is a function of PAR, R:FR and temperature as it equals the sum of growth (ER(t)) over the duration of organ growth.
In contrast to the published version of L-Cucumber, we modified the code to use daily environmental light information instead of data from a fixed date in the Quasi-Monte Carlo simulations.The L-Cucumber model utilizes an axis symmetric prototype leaf shape, where the estimated leaf area is equally distributed on four equilateral triangles.The left and right leaf halves are tilted down by 15 • to resemble observations of cucumber leaf shapes.A sketch of the original L-Cucumber model leaf is shown in Figure 1.The equilateral triangle side w (cm) is used in the leaf construction can be calculated from leaf area A L (cm 2 ) as follows: Vertical leaf orientation in L-Cucumber is modeled by an angle between midrib and vertical axis (z) in dependence of leaf age.Hence, leaf shape in this study is always considered relative to the midrib orientation as constructed in the virtual plant model eliminating these dimensions in the shape estimation.
More detailed information on the virtual plant model are provided in [5,6].

Digitized Leaves Dataset
All data on digitized leaves were collected in a greenhouse experiment with Cucumis sativus L., variety 'Aramon' in 2003.The experiment took place in the Institute of Biological Production Systems, Leibniz University Hannover, Germany (lat.52 • 23 35 N, long.9 • 42 09 E).It was set up as a salt stress experiment using a randomized complete block design with four treatment levels, four greenhouses and three measurement plants per treatment per greenhouse, resulting in a total of 48 measurement plants.A more detailed description of the experiment is given in [23]; we only provide the key parameters relevant to this study.
During a growth period of 36 d all plants were digitized up to seven times using a magnetic field digitizer (3D Digitizer Fastrak, Polhemus, Colchester, VT, USA).At the final date, leaves were also removed from the plant and measured with a leaf area meter (LI-3100; LI-COR, Lincoln, NE, USA).Digitalization of leaves followed a scheme of 17 unique landmarks (see Figure 2a).The landmarks 1 to 17 correspond to the standardized sequence for digitizing cucumber leaves also described in [23] (their Figure 1B).For example, landmark 1 represents the leaf base, the connection between leaf and petiole, and landmark 2 is the tip of the leaf.From a symmetrical point of view, seven landmark pairs can be defined as follows: 3-8, 4-9, 5-10, 6-11, 7-12, 14-16 and 15-17 (see Figure 2c).Except for the last two pairs, all landmarks are located on the leaf edge.Figure 2b shows the triangulation scheme applied for leaf area approximation from digitized data, considering a total of 20 triangles, 10 per each leaf half.The morphometric analysis in this study was limited to the 12 measurement plants from the control treatment to exclude saline stress effects on leaf shape.Thus, from 3866 digitized leaves we only considered 1125 digitized measurements for the shape analysis.However, 325 leaf area meter measurements across all treatments (data now shown), were used in validating the leaf area approximation by triangulation.
Figure 3 shows the distribution of estimated leaf area A L (cm 2 ) and corresponding equilateral triangle side w (cm) for 1125 digitized leaves used for shape estimation, with the majority being in the range of 150 cm 2 ≤ A L ≤ 2000 cm 2 or 9.3 cm ≤ w ≤ 34 cm.

Plant Morphological Data
Virtual plant simulations were compared against experimental data from two greenhouse studies conducted in May (E1) and June (E2) 2012 at the Department of Vegetable Crops, Geisenheim University, Germany (lat.49 • 59 N, long.7 • 58 E).Monitored environmental parameters include daily global radiation data and ambient temperature, with daily averages over the simulation time-frames for PAR of 430 ± 127 µmol m −2 s −1 and 388 ± 134 µmol m −2 s −1 and for temperature of 22.0 ± 1.5 • C and 22.8 ± 1.4 • C, respectively.
Plant morphological measurements were focused on internode length.Internode length was measured using a ruler for internodes from the fifth rank to the top of the plant.In both experiments four measurement dates were realized (in E1 at day 132, 137, 144, and 150 of year; in E2 at day 179, 182, 189, and 192 of year).Internodes at lower ranks were excluded as they had already been developed during cultivation, before the plants were treated under the experimental conditions.Furthermore, not to risk plant damage during measurements, the minimal considerable internode lengths was 3 cm.More details on the experimental setup are provided in [6,23].

Morphometric Leaf Shape Analysis
With the aim to improve the prototype leaf shape used in the L-Cucumber simulations, a method was developed that allows comparisons of the original leaf prototype with measurement data while considering the concept of fluctuating asymmetry.In contrast to most leaf shape analysis we were interested in the unaltered 3D-shape and not in a 2D-outline, as this is how the plant interacts with environmental light conditions.Furthermore, following the construction of leaves in the virtual plant model, we estimate the deviations from measurements and model leaf, thus, unlike in Procrustes superimposition [37], we limit the rotations in the alignment process.In addition, as we expect higher variability in landmarks further away from the leaf base, the well-known Pinocchio effect might cause problems, when we would apply a generalized Procrustes analysis [38][39][40].

Preprocessing
With two distinct fixed points from the original leaf prototype construction, the leaf base (P 1 ) and the leaf tip (P 2 ), we begin the transformation in our method with a translation of the leaf base to the origin and align the P 1 P 2 to the x-axis by applying yaw-and pitch rotations only.This allows us to analyze deviations from the original leaf prototype to the measured data and facilitates the implementation of a new shape in the virtual plant model.The alignment is similar to a three-dimensional equivalent of Bookstein's two-point registration approach, or a sliding baseline registration [39,[41][42][43], but with a notable difference in the scaling method.In our approach, scaling is based on the equilateral triangle side w, a measure derived from the approximated leaf area (see Equation ( 3) and Section 2.3.2).We do not scale leaves to unit centroid size (the square root of the sum of squared distances of all landmarks to their centroid) as the latter only is uncorrelated with shape, if landmark variation is relatively small and approximately isotropic [37,44].Centroid size is affected by natural folding and rolling deformations of leaves.Differently deformed leaves with equal initial surface areas do not have the same centroid size, and hence scaling would introduce artificial variance at landmarks unaffected by such deformations.For example, a flat original prototype leaf with a leaf area of A L = 100 cm 2 has a centroid size of 22.92 cm, estimated over all 17 landmarks, whereas a leaf, where both halves are folded downward to the maximum angle of 90 • has a centroid size of approx.16.44 cm.Hence, scaling to unit centroid size would lead to considerably different surfaces areas (sizes) of approx.0.19 cm 2 and approx.0.37 cm 2 , respectively.This would introduce variance in all three dimensions, although only two dimensions were altered.When scaling by equilateral triangle side w, estimated from the leaf area (Equation (3)), the surface area of both leaves will remain equal (here: w ≈ 7.60 cm).A general example of the differences introduced when scaling by centroid size or by w is given in the Supplementary Materials (Figure S1, Table S1).
We also fix landmark 1-the leaf base-to the origin, as it represents the connection between leaf and petiole, hence, there is no natural source of variation here we want to account for.When, in addition, alignment would be based on the centroid, as it is done during a GPA, Pinocchio effects are introduced at the leaf base, as the variance is equally distributed over all landmarks.
We decided to use w in the estimation of relative displacements between model leaf and measurement (see Section 2.3.2), as it was already present in the virtual plant model.Equation (3) shows that the same results could have been achieved when scaling by the square root of the leaf area instead, as w can also be interpreted as a fixed factor multiplied by the square root of the leaf area.Hence, provided that an estimation of the leaf area is available, a similar scaling could be used for other leaves, too.
The decision, not to roll the leaf around the x-axis in the transformation process is based on the premise to maintain the leaf surface orientation relative to the environmental light exposure.An exemplary alignment and scaling is shown in Figure 4, including the corresponding original model leaf shape with the same area, to which it is compared.

Morphometric Measures
To measure the deviation between the prototype leaf of the original L-Cucumber and the digitized leaves, we defined 17 similar landmarks on the model leaf (Figure 2c) and estimated the point-wise displacements in three dimensions as follows: The three dimensions are s x , the displacement along the x-axis, s y , the displacement in distance to the x-axis, and s z , the displacement measured as roll angle around the x-axis.Normalization for a scale-independent shape description was realized by normalizing the distance measurements s x and s y with the prototype's equilateral triangle side w (cf.Equation (3)), while the angle s z already is a relative measure.The reasons for choosing equilateral triangle side w over centroid size were already discussed in Section 2.3.1.Using two distances and an angle rather than distances between coordinates reflects the construction process of the original model leaf.It assures that leaf size is preserved when estimating the average displacement for each dimension and is related to the fact that we do not remove the roll-dimension around the x-axis in the alignment of leaves.
An example of the estimation of the three dimensions of displacement between model and measurement is given in Figure 5.

Symmetry Analysis
For the determination of leaf shape we included a symmetry analysis between paired landmarks, with the focus on accounting for fluctuating asymmetry.This is embedded in a robust Bayesian estimation [45,46] of signed differences between paired landmark shifts.Robust in this context means, that it utilizes the outlier tolerant properties of Student's t-distribution in the Bayesian estimation model in contrast to the Gaussian normal distribution [45].In this way, natural outliers and outliers due to human errors during digitalization are accounted for in the statistical analysis.Robust Bayesian estimation was performed using the BESTmcmc-function from BEST-package (v0.5.0) [46] in R (v3.4.0) [47].The model parameters were kept at the default values, to assess the models robust applicability.For details on this method we refer to Kruschke [45] and the BEST-package manual Kruschke and Meredith [46].
The following procedure was conducted for all landmark pairs and the three displacement dimensions to determine the mean shift including the information on symmetry: 1.
calculate the left-right signed difference between the paired landmark shift for a single dimension (s x , s y or s z ) 2.
conduct robust Bayesian estimation on the differences (n = 1125) to estimate posterior mean (ŝ) and posterior standard deviation ( σ) as measures for directional and fluctuating asymmetry if necessary: adjust average displacements using directional asymmetry information To detect antisymmetry in leaves the posterior distributions of the Bayesian estimation on the signed differences were subjected to visual inspection.For the unpaired landmarks on the midrib, 2 and 13, symmetry evaluation becomes obsolete and only a robust Bayesian estimation for the absolute average displacements is performed (n = 1125).Finally, with the results on average displacement between landmarks from the original L-Cucumber model leaf in all three dimensions, the new leaf shape is determined.

Virtual Plant Simulations
To assess the effect of the new leaf shape within virtual plant simulations, the L-Cucumber code [6] was modified to include the construction of the new model leaf.The leaf construction was implemented following the proposed triangulation with 20 triangles shown in Figure 2b.We exploited the edge length (w) information of the original L-Cucumber leaf prototype for rescaling the distances in order to build the new leaf with the correct surface area.
A set of 100 L-Cucumber simulations were run for both leaf shapes with a virtual stand of 15 (3 × 5) plants.To exclude edge effects from the simulation's cultivation area, data were taken from the central plant only and averaged over the 100 simulations.Furthermore, to resemble natural variation and get an estimate of the overall model performance, initial plant orientation was varied between simulations based on a latin hypercube sampling from plant-wise initial orientations between 0 • and 144 • .This orientation interval corresponds to the phyllotaxis of the plant.The 2/5-phyllotaxis indicated two rotations along the main stem for five consecutive leaves.The sixth (1 + 5) leaf has the same orientation as the first one.The angle 144 • equals 2/5 × 360 • .
To compare the simulation results against the experimental data from the greenhouse experiments E1 and E2 (Section 2.2.2) internode lengths from rank 5 and higher were exported for each time step (day) (see also [5,6]).Experimental data were available for four measurement dates per growth period (E1, E2).The corresponding simulations steps use in the comparison were steps 11, 16, 23, and 29, and 12, 15, 22, and 25, for E1 and E2, respectively.
In addition to the local internode measurements, simulations were compared based on summary data from the final time step including produced dry weight per plant and day (DW, g), leaf area per plant (LAP, cm 2 ), daily stem red-to-far-red ratio from internodes weighted by internode length R:FR Stem (-) and total plant length (L P , cm).
A robust Bayesian estimation [45] of relative differences between paired simulations (same initial plant orientation) was used to determine significant effects of leaf shape on simulation outcome.To compare the data across L-Cucumber simulations for both measurement periods (E1, E2) the results from the original leaf shape simulations were used as the reference.

Average Leaf Shape
Significant differences from the three-dimensional shape analysis between the original prototype leaf from L-Cucumber and the measurements were transferred into a new leaf shape.Figure 6 shows a comparison between the original leaf shape, the new shape and all 1125 measurements.The large spread in measurement data shows the high variability in the three-dimensional shape remaining after alignment.Comparing original and new shape, most distinct deviations were found at the tip and at the outer edges, corresponding to landmarks 2, 6 and 11 (cf.Figure 2a).In relation to the leaf base point-the fixed reference point, which represents the connection between leaf and petiole-we identified the tip to be approx.29 % (s x = 0.29w) further apart from the base then in the original prototype leaf.As leaf area (size) is retained, this elongation of the leaf along the midrib is compensated by a more pronounced heart-shaped outline.This is expressed in the outer edge landmarks, 6 and 11, by being closer to the base point in the new leaf model.Changes in curvature of the leaf are easier to recognize when looking at the three-dimensional representation in Figure 7.While there was only a single roll angle present in the original shape, in the new leaf shape the points in proximity to the leaf tip are closer or even above the xy-plane (Figure 6, side view).This is reflected by strongly positive s z -displacement, the change in roll-angle around the x-axis, estimated for the landmarks 3, 8, 4 and 9 (see Table 1).Distinct negative s z -displacements were found for the landmark pairs 6-11 and 7-12 with approximately −16.6 • to −19.9 • leading to a greater curvature from left to right leaf halves than the one assumed in the original L-Cucumber model leaf.
With s z = 0 and s y = 0.16w landmark 13 was found to have a vertical shift away from the midrib without having a tendency to one half of the leaf.This introduces a curvature along the midrib, where in the original prototype leaf the midrib is a straight line through the landmarks 1, 13 and 2. Shifts estimated at the inner landmark's (14 to 17) are supporting the leaf's three-dimensional curvature, as described before.
Focusing on the s x displacement along the x-axis, represented in the top view of the original and the new leaf shape (Figure 6), it can be noted that landmarks 5-15-1-17-10 still lie on an almost straight line through the origin (P 1 ) (s x ∈ {0.001w, 0.009w, 0.02w}, Table 1), as it was the case in the original L-Cucumber prototype shape (Figure 2c).This supports the assumption of folding and rolling deformations in relation to the midrib.

Leaf Symmetry
We found the new leaf shape to be almost perfectly mirror-symmetric.In Figure 8, where the right leaf half is mirrored onto the left half across the xz-plane, only minor deviations between paired landmarks can be found.This is the result of the robust Bayesian estimation on the differences in paired landmark's shifts (cf.Section 2.3.3).  1) indicate directional asymmetry.The shaded area represents the outline of the leaf.
Figure 9 shows the estimated posterior means and their 95 % HDI.No significant deviation from average symmetry was found for the landmark pairs 1, 2, 3 and 6, as indicated by the 95 % HDI overlapping zero in all three dimensions (s x , s y , s z ).Directional asymmetries were detected for the pairs 4, 5 and 7 in the s x and s y dimension.The estimated differences in means are in the range of 0.008w to 0.02w.This corresponds to a maximum deviation of 0.68 cm, considering leaf areas from 50 cm 2 to 2000 cm 2 (5 cm < w < 34 cm) (cf. Figure 3).Previous studies found measurement errors for similar magnetic field digitizers, calculated as the root-mean-square deviation (RMSD), to be in the range of 0.15 cm to 1.0 cm.Less deviation was found in laboratory trials compared to field experiments [23,48,49].Hence, we cannot fully disregard this differences in symmetry on the basis of measurement error, when interpreted as the region of practical equivalence (ROPE) [50].We therefore keep these results on directional asymmetry for the further study, although we do not expect significant effects on light interaction in the virtual plant simulations.The extent of fluctuating asymmetry is shown in Figure 10 as the estimated posterior standard deviations (ŝ x , ŝy , ŝz ) and the corresponding 95 % HDIs of the robust Bayesian estimation on differences in shifts of paired landmarks.For all landmark pairs fluctuating asymmetry is significantly different from zero, with values approximately ten times the extent of the maximum values for directional asymmetry (cf. Figure 9).This is not in contradiction with the findings on directional asymmetry-directional asymmetry and fluctuating asymmetry can both be present [27].It shows that there is no dimension that is free from fluctuating asymmetry.This is not surprising, as we are studying the three dimensional shape of the leaves, and did not use an alignment method that minimizes variability between corresponding landmarks, as for example in a GPA.Hence, the fluctuations in this study are indicators for the large variability in three-dimensional leaf deformation, i.e., leaf shape, due to folding, rolling, tilting or dangling of leaf parts.It should therefore not to be confused with the fluctuating asymmetry measured in flattened 2D leaf shapes.Hence, in contrast to analysis of developmental instabilities using 2D leaf shapes, the extent in fluctuating asymmetry in our study vastly exceeds the expected measurement error [23,51].This different perspective offers a chance to study 3D responses of leaf deformation on environmental (stress) conditions, although the extreme variability necessitates large sample sizes.A possibility to reduce variability, that we deliberately kept in the study, in order not to effect leaf orientation in relation to environmental light, would be to use an alignment more similar to Procrustes superimposition.Hence, such an analysis on 3D deformations would focus on slightly different objectives based on alternative assumptions.However, as we discussed in Section 2.3.1, scaling should still be based on measures related to the square root of leaf surface area rather than on centroid size.
In summary, we could successfully reconstruct a more realistic leaf shape from the experimental data by considering the concept of fluctuating asymmetry within the process with a robust, outlier-tolerant approach.Minor directional asymmetries were detected.Whether these have a natural cause or are a measurement error could not be verified in this context and needs further investigation.Fluctuating asymmetry was present in all landmark pairs and reflects a high variability in three-dimensional leaf shape, while no signs of antisymmetry were found.

Effects of Leaf Shape on Virtual Plant Simulations
In the virtual plant model internode growth is effected by both, light quantity and light quality-more precisely by the R:FR ratio-at the internodes' light receptors (see Equation ( 1)) [5,6].
Both, the original L-Cucumber leaf and the new prototype leaf virtual plant simulations perform well in comparison to experimental data on internode lengths (Figures 11 and 12) for the two different growth periods, E1/May and E2/June.Error estimation by root mean squared deviation (RMSD) per date indicate minor effects of leaf shape on internode development.Overall, we found slightly positive effects of the new leaf shape, especially on account of the June experiment (E2).Comparing the summary parameters of the final simulation step, similar patterns were found in the results representing different environmental conditions E1 and E2 (see Figure 13).The parameter directly connected to internode growth R:FR Stem is significantly effected by the new leaf shape in both simulations scenarios, based on the 95 % HDI around the posterior mean not overlapping zero and estimated mean differences of approx.5 % to 8 %.We found higher R:FR Stem in the simulations with the new leaf shape associated with reduced internode growth, which is consistent with lower total plant length L P (see Figure 13).We also observed a reduction in DW by over 5 % using the new prototype shape.As this parameter is linked to a decrease in photosynthetic activity derived from the amount of PAR reaching the leaves' surfaces [5,22], we can conclude, that the new leaf shape induces more inter-or intra-plant shading in comparison to the simple original L-Cucumber shape.In a simulation study Sarlikioti et al. [34] found higher canopy photosynthesis rates for more elongated leaves compared to short-wide leaves.This might promote the hypothesis of increased intra-plant shading as a cause for reduced photosynthetic activity of the new leaf shape simulations, possibly due to the redistribution of model leaf area towards the leaf tip.Hence, effects of planting density might become more pronounced utilizing the new leaf shape in virtual plant studies [52].However, finding possible causes for this shading such as the increase in leaf curvature or the extent of the leaf tip is beyond the scope of this study.Furthermore, as an inter-simulation measure this supports our initial expectation regarding effects of leaf shape variation on simulation results, but is not a measure of leaf shape quality.
Comparing the computational expenses of original and new leaf shape simulations, we found an increase in computational time from 13 min to 31 min for ten parallel simulations on an Intel R Xeon R CPU E5-2660v3 with 64 GB Ram.This is mainly caused by an increased effort in Quasi-Monte Carlo light modeling, where, instead of the original four triangles, 20 triangles are evaluated for the new shape.Although computational times almost tripled, they are still on the scale of minutes and hence irrelevant in comparison to real life plant development cycles.
In summary, we found positive effects of the new leaf shape when comparing simulation results against experimental data on internode growth.This finding suggests that the effort for identifying more realistic leaf shapes is beneficial for simulation studies.Significant differences in other growth parameters between simulations with original and new prototype leaves shows the sensitivity of the functional-structural plant model on environmental effects and calls for studies considering further model extensions, e.g., on leaf movement (phototropism) [15].Relative deviation of summary parameters from new leaf shape simulations at the final simulation step in reference to the L-Cucumber simulations using the original shape for both experimental conditions (E1, E2).Posterior means and 95 % HDI estimated by robust Bayesian estimation of relative differences from 100 simulations varying in initial plant orientation.

Conclusions
We found significant effects on growth parameters related to canopy light interaction using a more realistic leaf shape model in the virtual plant model L-Cucumber.The proposed method to determine the leaf shape from three-dimensional measurements reflects leaf shape and orientation in relation to environmental light interaction.Assuming a mirror-symmetric leaf shape, the concept of fluctuating asymmetry was adapted in the robust estimation of landmark differences.It revealed only minor directional asymmetries, but a high variability in all dimensions, related to the various possible leaf deformations in the three-dimensional space, such as rolling, folding or tilting.Since knowing local light fluctuations affects plant development in reality and in the simulations, we conclude that it is advisable to identify more realistic leaf shapes for virtual plant studies.Further research will therefore focus on the identification of three-dimensional leaf shape development and its implementation as a dynamic model in L-Cucumber.Aspects of phototropism could also be of significance to simulation results and offer potential for additional studies.The effects of high variability observed in the digitized leaves could be addressed by utilizing stochastic models to achieve non-uniform leaf shapes in the virtual plant simulations.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-8994/10/7/278/s1, Figure S1: Different folding angles (0 • , 45 • , 90 • ) of a simplified leaf (isosceles triangle) with four landmarks and their centroid (arithmetic mean of all landmarks points in the shape), Table S1: Dimensionless landmark distances to the centroid (dC) for different folding angles (0 • , 45 • , 90 • ) (see Supplementary Figure 1) of a simplified leaf (isosceles triangle, A = 0.5 ) for comparing centroid size (CS) scaling with the proposed scaling factor w = A L √ 3 , which is related to the square root of the surface area.Centroid size is defined as the square root of the sum of squared distances of all landmarks to their centroid.

Figure 1 .
Figure 1.Top and front view of the original prototype leaf shape used in L-Cucumber.Four equilateral triangles.Midrib defines symmetry axis.Leaf halves are rotated (rolled left/right) around the midrib to form a 15 • downward angle.

Figure 2 .
Figure 2. (a) Seventeen unique landmarks in the digitalization process of cucumber leaves; (b) Triangulation scheme for leaf area calculation; (c) 17 similar landmarks assigned to the original L-Cucumber prototype leaf.Dashed arrows connect symmetrical landmark pairs.(Own representation based on Kahlen and Stützel [23].)

Figure 4 .
Figure 4. Exemplary preprocessing transformation of a measurement leaf to align with the corresponding prototype leaf shape of the original L-Cucumber model.Alignment of P 1 P 2 by translation to origin, and yaw-and pitch rotation.

5 l 5 ≈ l 5 +α 5 ≈ α 5 +Figure 5 .
Figure 5. Example of the point-wise displacements in three dimensions between a measurement point 5 (blue) and its prototype equivalent 5 (red) (a) s x -displacement along the x-axis; (b) s y -displacement in distance to the x-axis; (c) s z -displacement in roll angle around the x-axis.
highest probability density interval (HDI) of ŝ includes zero ⇒ no directional asymmetry (b) if 95 % highest probability density interval (HDI) of ŝ excludes zero ⇒ directional asymmetry 3. conduct robust Bayesian estimation on the pooled left-right data (n = 2250) to estimate overall posterior mean (absolute average displacement) 4.

Figure 6 .
Figure 6.Size-free comparison of new leaf shape, original L-Cucumber shape and all 1125 measurements by their outlines.Additionally, closed and open circles indicate the landmark positions in the new and original shape.Axis limits represent maximum spread of measurement data.

Figure 7 .
Figure 7. Three-dimensional impression of equal-sized original (left) versus new (right) leaf shape.Color representation is based on pairs of corresponding triangles from the left and right halve.

Figure 8 .
Figure 8. Representation of new leaf shape's symmetry by mirroring the right (R) leaf half onto the left (L) half across the xz-plane.Non-overlapping paired landmarks (cf.Table1) indicate directional asymmetry.The shaded area represents the outline of the leaf.

Figure 11 .Figure 12 .
Figure 11.Internode length from experiment E1 for rank 5 to maximal rank 20 against L-Cucumber simulation results at four timepoints for original and new leaf shape.(A) Rank-wise average and standard deviation of 5 experimental measurements or 100 simulations.(B) Root-mean-square deviation per day.

Figure 13 .
Figure 13.Relative deviation of summary parameters from new leaf shape simulations at the final simulation step in reference to the L-Cucumber simulations using the original shape for both experimental conditions (E1, E2).Posterior means and 95 % HDI estimated by robust Bayesian estimation of relative differences from 100 simulations varying in initial plant orientation.

Table 1 .
Overall posterior means (n = 2250) of point-wise three-dimensional displacement (s x , s y , s z ) for each landmark with consideration of differences in pairs from robust Bayesian estimation (no zero-overlap of the 95 % HDI) (see Figure9).A single value for a landmark pair indicates no directional asymmetry.