Development and Stochastic Validation of a Parameterized Model of Maize Stalk Flexure and Buckling

Maize stalk lodging is the structural failure of the stalk prior to harvest and is a major problem for maize (corn) producers and plant breeders. To address this problem, it is critical to understand precisely how geometric and material parameters of the maize stalk influence stalk strength. Computational models could be a powerful tool in such investigations, but current methods of creating computational models are costly, time-consuming, and most importantly, do not provide parameterized control of the maize stalk parameters. The purpose of this study was to develop and validate a parameterized three-dimensional model of the maize stalk. The parameterized model provides independent control over all aspects of the maize stalk geometry and material properties. The model accurately captures the shape of actual maize stalks and is predictive of maize stalk stiffness and strength. The model was validated using stochastic sampling of material properties to account for uncertainty in the values and influence of mechanical tissue properties. Results indicated that buckling is influenced by material properties to a greater extent that flexural stiffness. Finally, we demonstrate that this model can be used to create an unlimited number of synthetic stalks from within the parameter space. This model will enable the future implementation of parameter sweep studies, sensitivity analysis and optimization studies, and can be used to create computational models of maize stalks with any desired combination of geometric and material properties.


INTRODUCTION
Maize stalk lodging (wind-induced failure of the stalk) represents a major source of loss for farmers (Flint-Garcia et al., 2003;Duvick, 2005). In late-season lodging, maize stalks fail due to localized (Brazier) buckling which manifests as a characteristic "creasing" of the stem just above a node (Robertson et al., 2015) The term "buckling" refers to a sudden change in shape of a structural component under loading, often resulting in collapse. Brazier buckling involves localized shape change and should not be confused with the more well-known Euler buckling, which applies to large-scale buckling of columns. The creasing of a drinking straw that occurs when the straw is gently bent is the most commonplace example of Brazier buckling.
Stalk bending strength refers to the maximum bending load that can be applied to a stalk prior to the onset of buckling. Stalk bending strength is a major factor in determining lodging resistance and has been shown in both empirical and computational studies to be closely related to stalk morphology (Gomez et al., 2017;Robertson et al., 2017Robertson et al., , 2022. In this paper, "strength" always refers to ultimate bending strength. Tissue strength is also of interest in predicting ultimate strength but tissue strength is not addressed in this study. It has been proposed that relatively minor changes in stalk morphology might lead to significant increases in stalk strength (Von Forell et al., 2015). To improve maize stalk strength through breeding and/or biotechnology, two key advances are needed: (1) the ability to rapidly measure the morphology of the maize stalk, and (2) insight into which morphological features are most closely related to maize stalk strength. This study addresses the second need by developing a method for manipulating the morphology of maize stalks.
The finite element method has previously been used to create models of maize stalks (Stubbs, 2019;Stubbs et al., 2019Stubbs et al., , 2020Stubbs et al., , 2022Von Forell et al., 2015). In the broader biomechanics literature, the finite element method is a powerful and pervasive tool. In human biomechanics, the finite element method has been used to study a wide range of topics, including the human knee, foot, pelvis, liver, brain, heart arteries, and many others (Andersen et al., 2021;Besnault et al., 1998;Erdemir et al., 2006). There are two primary methods used to create specimen-specific models. First, image-based modeling is commonly used to create models directly from imaging data such as CT or MRI scans. Such models typically provide very high geometric fidelity and are thus very close approximations to the specimens from which they were derived. These models typically require a significant amount of manual segmentation and modeling effort by a trained technician (Patil & Deore, 2013). As a result, these models are costly and timeconsuming to create. Moreover, once created, the geometry of these models are essentially static and cannot be manipulated to investigate the influence of geometry (Cook et al., 2014).
To overcome these limitations, parameterized modeling schemes have been developed. Parameterized models typically utilize geometric features and/or landmarks of the system to be modeled (Danelson et al., 2009;Maurel et al., 1997;Sigal et al., 2010). Parameterized models provide good approximations of the geometry of the desired specimen and provide control over model shape through the adjustment of model parameters. This approach dramatically expands the utility of these models by enabling parametric studies, sensitivity studies, and optimization studies. However, the parameterization scheme typically requires geometric approximations, so parameterized models do not allow the same level of geometric detail as image-based models. Ottesen et al., (2022) used machine learning techniques to develop a twodimensional parameterization of the maize cross-section. The Ottesen model was based on a simple ellipse as the exterior boundary of the stalk and an offset ellipse as the boundary between rind and pith tissue. Finer geometric features such as the ear groove and various types of asymmetries in the model were captured using principal components derived from actual cross sections. Ottesen et al. reported that any desired level of geometric fidelity could be obtained by adding additional principal components.
The model was validated under a variety of loading cases, including bending, torsion, axial tension/compression, and transverse compression. Notably, the ellipse alone (no principal components) provided high predictive accuracy in all cases and required just three geometric parameters: major diameter, minor diameter, and rind thickness. The primary limitation of the Ottesen model is that it did not account for axial variation of the maize stalk.
A fully parameterized, three-dimensional model of the maize stalk would allow control over individual geometric features. Such a model could be used to determine which geometric features of the maize stalk are most closely related to maize stalk strength. It could also be used to determine the optimal stalk architecture that balances biomass with both strength and flexibility. The purpose of this study was to develop a three-dimensional parameterized model of the maize stalk. We validated this model in three ways. First, by assessing its geometric similarity to original stalks. Second, by determining the ability of this model to accurately predict flexural stiffness, and third, by assessing the ability of this model to predict the strength of maize stalks. In contrast with prior studies, stochastic (randomized) variation in material properties was used to account for the influence of material uncertainties on the validation process. The resulting model opens many avenues for advanced computational analyses of maize stalk strength.

Overview
This study is based upon two important data sources from previous studies: (1) CT scans of maize stalks, and (2) corresponding physical 3-point bending tests (Robertson et al., 2016. Previous research has confirmed that 3-point bending tests can be used to asses stalk bending strength because these tests induce the same type of buckling failure observed in natural lodging (Robertson et al., 2014;Robertson et al., 2015a;Robertson, et al., 2015b) Figure 1 shows the physical testing arrangements along with sample CT scan images. This figure also depicts the creation of CT-based finite-element models , the use of CT scan data to create 2D parameterized models (Ottesen et al., 2022) as well as the 3D parameterized models developed in this study.

Empirical Data
Three-point bending tests and the CT-scanning process were described in previous studies (Robertson et al., 2016. Briefly summarized, approximately 900 maize stalks were scanned and then tested in three-point bending. During bending tests, stalks were gradually loaded until localized buckling occurred. This testing approach provided empirical measurement of two mechanical features of each stalk: flexural stiffness (also known as flexural rigidity), and ultimate bending strength. CT-scans provided detailed three-dimensional data on each maize stalk at a resolution of 78 um/voxel. This resolution typically results in over 200 pixels across the minor diameter of the stalk and over 100k pixels within each cross-section. Representative CT crosssections are shown in Figure 1.
The CT data for each stalk consisted of approximately 1500 cross-sectional images. Because adjacent images were nearly identical to each other, a non-uniform sampling scheme was used. In this scheme, 46 cross-sectional images were extracted from each stalk using higher sampling density near the node, and lower sampling density in the internodes. The 46 sampling points can be seen as light gray vertical lines that intersect the longitudinal CT cross-section image in Figure 1.

Ellipse-based axial sweep
Each cross-section of each stalk was approximated as an ellipse as described previously (Ottesen et al., 2022). The interior ellipse was offset from the exterior boundary by the rind thickness. The rind thickness was determined as the average distance between the interior and exterior boundaries of the maize stalk cross-section . This process produced profiles for the major diameter a, minor diameter b, and rind thickness t as functions of axial location z for 900 individual stalks.
For each individual stalk, these three profiles define a three-dimensional shape that closely approximates the original stalk geometry. Figure 2 provides an illustration of the swept-ellipse geometry.

Parameter Profiles
The profiles, a(z), b(z), and t(z), exhibited distinct features that were present in virtually all stalks. Figure 3 depicts the data distributions for each parameter profile, with dark shading indicating high data density and light shading indicating low data density. The solid white line in the center of each chart represents the average profile.

Profile Parameterization
Machine learning techniques and landmark identification were used to parameterize the essential features of each parameter profile. This was accomplished by iterating between principal component analysis (also known as empirical eigenfunction analysis) and landmark identification. Early attempts applied principal component analysis to each profile separately, or to all three profiles simultaneously. Although this approach does provide an immediate parameterization based on variance decomposition (Jackson, 2003), this method does not enable independent control of physical features and was therefore abandoned. For similar reasons, polynomial-based regression approaches were excluded because this provides no physical connection between coefficient values and physical features of the corn stalk.
Landmarks were selected based on consistent and distinctive features which were found to be present in virtually all stalk profiles. Most landmarks were peaks or valleys in the respective profiles. Each landmark was defined by two coordinate values. The first coordinate was the axial location z. The z-axis origin was defined as the central peak in the rind thickness profile. This feature was chosen as a reference because it was (a) the most distinct of all landmarks, and (b) identifies the center of each node (the dividing line between stalk segments). The second coordinate was the perpendicular distance of each profile (major diameter, minor diameter, and rind thickness) from the z-axis.  Custom software was developed to automatically identify these landmarks for each of the 980 x 3 = 2940 profiles used in the study. Principal component analysis was then used to capture and capture and parameterize the transition patterns between landmarks. A local, normalized coordinate system was defined for each neighboring pair of landmark points. The original and subsequent coordinate systems are shown in Figure 3. The transverse values retained their original units but the transverse value of the first landmark was subtracted so that the coordinate value of the first landmark was (0,0). The local coordinate system is shown in Figure 3B.
Next, the landmark points were connected in the local coordinate system by a straight line. The final "subtracted" coordinate system captures the vertical distance between the line connecting the landmarks and individual points on the profile. The subtracted coordinate system is shown in Figure 5C. typically captured approximately 96% of the variation in a transition pattern. This approach allows any desired geometric fidelity to be captured by simply adding additional eigenfunction terms.
To summarize, the parameterized model of maize stalk morphology consists of an axially varying elliptical cross-section. The variation in ellipse shape is defined by three parameter profiles: major diameter, minor diameter, and rind thickness ( Figure 2).
Each of these profiles is defined by six landmarks (Figure 4). Straight lines connect adjacent landmarks and transition patterns between landmarks are captured by linearly scaled empirical eigenfunctions ( Figure 5). Any number of empirical eigenfunctions can be used in this scheme, which allows for complete reconstruction of the axial patterns of ellipse profiles when all eigenfunctions are utilized. We found that the first empirical eigenfunction captured the vast majority of geometric variation, so only one eigenfunction was used in this study. With a single empirical eigenfunction, the shape of the parameterized 3D model is defined by 51 parameters, each of which can be independently adjusted/controlled. These parameters are shown in Figure 4 and listed Table S1 of the supplementary data that accompanies this paper.

Generating Parameterized Models
The parameterized model provides the ability to fit the parameterized model to actual maize specimens. The parameterized model was fit to each of the original 900 maize stalks in the database. Once the fitting process was completed, a subset of corresponding finite-element models were created.
It is also possible to use the parameterized model to generate an arbitrary number of novel "synthetic" stalk models. This was accomplished using the normal copula method (McClarren et al., 2018). This method produces a synthetic population in which the mean values, standard deviations, correlations, and parameter distributions shapes of the synthetic population match those of the original data set. This approach preserves key relationships between parameters. Synthetic models can also be generated to meet any arbitrary set of pre-determined statistical attributes. The methods for model fitting, modification, and generating synthetic stalks are described in the following sections.

Creating Specimen-Specific Parameterized Models
Several steps were required to fit the parameterized model to an individual stalk specimen. First, ellipses were used to capture each cross-section of the original stalk geometry. The ellipse fitting resulted in axial profiles for the major diameter, minor diameter, and rind thickness similar to those shown in Figure 4. The landmarks were then identified. The coordinates of landmark points provided 36 of the 51 parameters (see Table S1).
Next, transitions between each pair of landmark points were captured using empirical eigenfunctions. Each transition pattern between landmark points was converted to local and subtracted coordinate systems as described in Equations 1 -4.
This produced a vector of {y sub }values for each transition pattern. The following equation was then used to approximate {y sub }:

Creating a Population of Synthetic Models
The Gaussian copula method was used to create a population of synthetic models that had the similar statistical characteristics as the original set of 900 stalks. First, the parameterized model was fit to all 900 stalks in the original data set. This produced a 900x51 matrix of parameter values. This matrix was then used to calculate the 51x51 correlation matrix, R as well as empirical cumulative density functions for each parameter.
The Gaussian copula is formed by performing Choelesky factorization on the R matrix to obtain a lower diagonal matrix L, sampling a set of standard normal random values, and then mapping those values into new coordinates according to the relationships between variables captured by R, and probability density characteristics captured by the cumulative density functions. This process is described in more detail in (McClarren et al., 2018).
This process created a new population of N stalks. The value N can be set to any desired value. The quality of Gaussian copula was assessed by creating a new synthetic set of 900 stalks. The mean and standard deviation for both the original and synthetic sets were compared. The correlation matrix of the synthetic set was also compared termby-term with the original correlation matrix. Finally, the empirical cumulative density functions of the synthetic stalk set were compared term-by-term to the original empirical density functions.

Finite Element Models
The finite element method was used to create parameterized models of maize stalks. Models were generated using Abaqus/CAE 2022 by defining the appropriate material properties, mesh, and boundary conditions, as described in the following sections.

Mesh
The mesh of each parameterized finite element model was broken into three sections: top end, node region, and bottom end. The node region spans 10 mm with the fixed boundary condition in the center. The top and bottom ends were meshed using a hexahedral mesh due to the simplicity/regularity of the geometry in this region. The nodal region was meshed using a tetrahedral mesh to account for the more complex internal geometry of this region (see Figure 6). The rind was given high mesh density while the pith was given low mesh density. This was done as most stresses in bending occur in the rind. An extensive mesh convergence study was performed to ensure that this mesh produced mesh-converged results.

Boundary Conditions
The boundary conditions for the parameterized finite element 3D model match the shear and moment loads that were present upon the modeled segment during physical 3-point bend testing. This approach has been previously used and validated . Because the finite element models only account for a portion of the stalk, as shown in Figure 1, the relevant shear and moment loads were applied to the end faces of the finite element model. The loading anvil at the center of the stalk was held stationary with no degrees of freedom allowed at the point of contact between the anvil and the stalk to prevent rigid body rotation/translation. The expressions for shear and moment values are provided in Table 1 and the diagram in Table 1 shows these applied loads. Using the measured test values, the forces and moments for each finite element model was calculated using the equations in Table 1. In the finite element model, the central anvil of the finite element models was assumed to be fixed in six degrees of freedom.
distance between left-hand support and applied load in the 3-point bending test B distance between right-hand support and applied load in the 3-point bending test M fail maximum bending moment applied during physical 3-point bending tests

Mechanical Tissue Properties
Both pith and rind tissues are characterized by fibers which run longitudinally through the stalk. Plant tissues typically experience low levels of strain. As a result, both rind and pith tissues were modeled as linearly elastic transversely isotropic, as has been done in several previous studies (Stubbs et al., 2018Von Forell et al., 2015). This material model requires five independent material properties for each tissue type, consisting of the modulus of elasticity in the fiber direction ( E 3 ), the modulus of elasticity in the transverse plane ( E 1 and E 2 ), the shear modulus in the fiber direction ( G 13 and G 23 ), the shear modulus in the transverse plane ( G 12 ) and Poisson's ratio (ν).
Currently, only E 1 , E 2 , and E 3 in the rind as well as E 1 and E 2 in the pith have been

Stochastic Sampling of Material Property Values
Because not all the properties of maize tissue have been measured experimentally, some estimation is always required to obtain a full set of material properties. In previous studies, these estimates have had a limited scope: typically a single estimate was used for all other tissues (Stubbs et al., 2018Von Forell et al., 2015) An unnaturally low level of biological variation in models can obscure the influence of material properties (Cook et al., 2014). In contrast, a direct treatment of variation and uncertainty provides much more reliable insights (Nelson et al., 2019) In this study, realistic ranges were estimated for all material properties. Values for individual models were randomly sampled to account for the fact that these material properties are subject to significant levels of uncertainty. While many mechanical properties of maize tissues have not been measured, the relevant properties have been measured for many varieties of wood, both hardwood and softwood (Green et al., 1999).
The longitudinal modulus of elasticity across many wood species has been reported by Green et al. as ranging from 4 Gpa to 16 GPa. This range corresponds almost exactly to the range reported by Al-Zube et al, 2018 for maize rind tissue: 3 -17 GPa. Estimates for ranges of the missing maize material properties were based upon the ranges and ratios of the same material properties reported for wood (Green et al., 1999). To aid in this process, we also referred to thermodynamic constraints which exist when using a transversely isotropic material model. These constraints are (Lempriere, 1968): Here ν 12 represents passive deformation in the 1 direction due to applied tension in the 2 direction.
The most important material property of maize stalks is the longitudinal modulus of the rind (Ottesen et al., 2022). In this study, the longitudinal modulus of the rind was measured directly for each specimen using the bending test method described in Al Zube et al 2018. The remaining properties were randomly sampled from the ranges provided in Table 2 below. Normal distributions were used where possible. If insufficient information was available, uniform distributions were used.

Finite Element Analyses
Two analyses were performed for each model: flexural stiffness and linear buckling analysis. The flexural stiffness of the model was obtained using a static structural analysis. The force/deformation relationship from simulations was used to calculate an equivalent flexural stiffness which was then compared to the measured value of flexural stiffness. Critical buckling loads were obtained using a linear buckling analysis. This type of analysis solves for instabilities in the structure and can be used for traditional column buckling or localized buckling. In this study, linear buckling analysis provided the loat at which localized buckling instabilities occurred. This load represents and estimate of the ultimate bending strength of the stalk. It should be noted that linear buckling analysis always overpredicts ultimate bending strength, but it is computationally much more efficient than dynamic explicit analyses (Bolton, 2020;Xu et al., 2013). Both flexure and linear buckling methods have been previously validated for this purpose .

Cross-Fold Validation of Parameterization Method
When using machine learning methods, a portion of the data is used to train the model, with the remaining data held in reserve for validation of the model. The geometric maize stalk data used in this study was obtained from five different hybrid varieties. As these varieties are genetically distinct, the hybrid type was used as the basis for partitioning between training and validation sets. Rather than arbitrarily choosing certain hybrids for training and others for validation, we used the more robust crossfold validation approach to investigate all possible training/validation combinations (Kolodiazhnyi, 2020). In each training/validation set, the training set was composed of four hybrids with the remaining hybrid used for validation. Cross-validation consisted of repeating this process five times (one for each hybrid was used as the validation set).
Model validation was quantified by computing the R 2 value between validation set parameter profiles and the corresponding "fit" obtained when using the model derived from the training data. This was performed for each stalk in the validation set and the profiles obtained after fitting the stalk with the model.

Validation of Finite Element Model Predictions
Finite element models were validated by comparing their predictions to the corresponding empirical test data (flexural stiffness and ultimate strength). Each parameterized specimen specific model was created by fitting the parameterized ellipse model to the geometry of the corresponding CT scan. The specimen-specific longitudinal modulus of the rind (E 3 ) was also obtained from experimental test data. To account for uncertainty in the remaining material properties, ten versions of each specimen-specific model were created. Each version utilized the same geometry and rind modulus, but different random values of the remaining material properties. This approach helps ensure that the validation results are not an artifact of a narrowly-specified set of material properties.

Geometric cross-validation
The parameterized model geometry was compared to the actual stalk geometry by comparing actual profiles for major diameter, minor diameter, and rind thickness to their parameterized model counterparts. For all 5 hybrids and across all three profiles, the R 2 values between original profiles and the parameterized profiles averaged 0.96 with the lowest value across all hybrids and profiles of 0.91. This indicates that the parameterization approach that has been developed is sufficiently flexible to accurately capture a range of maize stalk geometries. It should be noted here that each hybrid was grown at a wide range of planting densities and that planting density is known to have a strong effect on stalk morphology (Jun et al., 2017;Sher et al., 2017). The crossvalidation statistics across validation sets is provided in Table 3. Table 3: Cross-validation R 2 values between the CT-based profile and corresponding parameterized models. Sample sizes varied between sets. The smallest sample size was n =188.

Model creation time
The parameterized model automates the process of model creation, allowing models to be created programmatically instead of manually. The previous method required CT scanning and manual processing of the CT data to create each finiteelement model. Including scanning time and manual manipulation of data, this process required approximately 45 minutes per model/specimen.
The parameterized model is capable of building models at a rate of 40 seconds/model. This is 70x faster than previously possible and no manual model manipulation is required. Within a few minutes, the computer code that controls the model building process can be programmed to create hundreds or thousands of models.
Model solve time is unchanged with respect to prior models. Using a Dell Precision 5820 Tower Workstation, each flexural analysis required 30 minutes of solving time and each buckling analysis required 90 minutes. Excluding mesh convergence, the results described below represent the creation and analysis of 1160 unique finite-element models.

Predictive accuracy of finite-element models
Direct comparisons between experimental results and simulated behavior were performed to assess the validity of the parameterized model. A set of 59 stalks were randomly selected from the original set of 900 tested stalks. For each of these stalk geometries, a sample of 10 replicates were created, each having a unique set of plausible material properties. This resulted in 590 unique finite-element models. Across the 590 simulation results, predicted flexural stiffness compared very favorably to the measured flexural stiffness, with an R 2 value of 0.97 and a best fit line of slope 1.11. The average response across all replicates was calculated for each stalk geometry and also compared to measured flexural stiffness values. Both data sets are shown graphically in Figure 7.
Random variations in material properties had a relatively minor influence on the predicted values of flexural stiffness. This effect was quantified using variance decomposition. Only 1% of the total variance was due to variation in material properties.
This suggests that flexural stiffness is primarily dependent upon geometry and the longitudinal modulus of the rind (which was not varied because it was measured for each specimen). A similar comparison was performed for the prediction of ultimate bending strength, which was assessed using linear buckling analysis. It is well-known that linear buckling analysis tends to overpredict actual buckling loads (Bolton, 2020;Xu et al., 2013). In a previous study using CT-based maize stalk geometries, linear buckling analysis was correlated with ultimate strength, but overpredicted it by a factor of approximately 10x . In this study, linear buckling results were also well-correlated with ultimate strength but a more moderate over-prediction of about 5x was obtained. Results for the full set of simulated stalks as well as replicate means are shown graphically in Figure 8. With material variation included, the R 2 value between simulations and test results was 0.56. With material variation excluded (replicate averaging), the R 2 value between simulations and test results increased to 0.73. Variation in material properties had a very noticeable influence on predicted ultimate strength. Using variance decomposition, it was found that 20% of the variation in predicted ultimate strength was caused by material property variation. This explains the broader spread of data in Figure 8 as compared with Figure 7. The much stronger influence of material properties suggests that material properties have a significant effect on ultimate strength. confidence intervals on the group mean for each stalk geometry.

Assessment of the Synthetic Population
The synthetic population of 900 stalks exhibited nearly identical statistical characteristics as the original population, thus indicating that the Gaussian copula method can be used to create unique parameterized maize stalk models. Comparisons between the statistical characteristics of the original data set and the synthetic data set are provided in Figure 9 below. As shown in the upper left-hand panel of Figure 9, the mean parameter values of the original set were virtually identical to the mean parameter values of the synthetic data set. The same was true for the standard deviations and parameter correlation values. Finally, term-by-term comparisons of the original and synthetic cumulative density functions were also very similar (lower right panel of Figure 9).

Using the parameterized model in studies on maize biomechanics
The parameterized model can be used in two main ways. First, the model can be used to provide a close approximation to the geometry of a real maize stalk. This is The specimen-specific "fitting" approach requires detailed geometric information on the stalk to be modeled. Such data could be obtained from a CT scan or from photographs. Once the parameterized model has been used to match a physical specimen, the behavior of the specimen specific model can be simulated. In addition, the model can be modified to assess sensitivity of the specimen to hypothetical changes in geometry and material properties. Alternatively, optimization can be performed, using the specimen-specific model as a starting point.
The second approach is to create synthetic stalk populations without the need for specimen-specific data. This is done by specifying the set of 51 parameter values. In this approach, the chosen values must be mutually compatible, and profiles should be representative of actual maize stalks. However, the data that was used to create the parameterized model can be used as a starting point.

Using this approach to model other grain species
All grains rely upon a similar geometric architecture. Thus, it may be possible to directly apply the parameterized model developed in this study to other grain species such as sorghum, wheat, oats, and rice. More research would be needed to determine the applicability of this model to other species. Alternatively, the methods described in this study could certainly be used to develop a parameterization that is unique for each grain species.

Model Limitations and Future improvements
The model described in this study captures gross anatomical features, but does not model smaller features such as internal vascular bundles, internal transitions between internodes, and the tendency of the stalk cross-section to change dramatically when primary or secondary ears are developed. Additional research will be needed in these areas.
This model relies upon the assumption that each cross-section of the stalk is purely elliptical. In our previous research, we showed that this assumption provides a favorable balance between model complexity and predictive power. In some loading conditions, the presence of the ear groove may be necessary to accurately predict stalk mechanics. If the ear groove is needed, additional modeling efforts would be required to add a parameterized ear groove to the model. We previously described a parameterized cross-sectional model of the maize stalk that includes the ear groove and other minor geometric factors. That model could be combined with the model in this paper to simulate the ear groove with relatively high fidelity.
This model assumes that the maize stalk can be accurately modeled using two distinct tissue types: rind and pith. In reality, the boundary between rind tissue and pith tissue is gradual rather than distinct. More research would be needed to capture the functionally graded aspect of maize stalks, but the exterior shape of the gross anatomy provided by the model outlined in this paper could be used as a starting point for such a model.

CONCLUSION
This study provides a major improvement over the time-consuming process previously required to create specimen-specific models of maize stalks. As evidenced by validation data, this new model accurately captures the behavior and trends observed in empirical tests of maize stalks. Most importantly, this new model allows individual physical features of the maize stalk to be independently controlled, thus opening the door to future studies focused on geometric sensitivity analysis, optimization studies, and the creation of synthetic populations of maize stalks for future study.

Author Contributions:
NWL identified landmarks, created the landmark identification algorithm, and created the process used to convert profile paths to solid models. MO performed the initial principal component analysis (PCA), integrated PCA results with landmarks to produce parameterized profile paths. JC created the code used to automate the creation of finite element models and performed the finite element analyses. RH integrated and revised the codes created by MO and NWL so that they could be used in an automated manner.
DC supervised each stage of the project. All authors were involved in the writing of the manuscript.

Availability of Data and Models:
We encourage other researchers to use the parameterized maize stalk model. While the model is not currently available in an online format, the model and data used in this study is available from the corresponding author upon reasonable request.