Development and validation of a statistical shape modeling-based finite element model of the cervical spine under low-level multiple direction loading conditions

Introduction

Cervical spine injuries are of significant concern in all trauma injuries, particularly given the potential for spinal cord injury in unstable injuries, with an estimated 42% of all cervical spine injuries being unstable (Milby et al., 2008). Recent military conflicts (i.e., Iraq and Afghanistan) have demonstrated the substantial risk of spinal injuries for the modern warfighter due to the high occurrence of traumatic injuries in combat. Explosive mechanisms resulted in inertial injuries due to direct blast injury and/or subsequent impacts in 75–78% of combat casualties (Belmont et al., 2012) and 28–39% of these injuries were suffered in the head/neck (Wade et al., 2007; Owens et al., 2008; Belmont et al., 2012). Although muscle strain is the most common injury in the neck, warfighters have experienced vertebral compression fractures and fracture of the spinous process in the lower cervical vertebrae, as well as interspinous ligament injuries in the lower cervical spine (Anderson, 1988; Schall, 1989; Coakwell et al., 2004).

Mechanisms of bony or ligamentous injury in the cervical spine and elsewhere have been investigated using finite element modeling, a common tool of structural analysts. Finite element models of the spine or spinal motion segments often employ generalized anatomical geometry or subject-specific geometry where the morphology of spinal components (e.g., vertebrae, intervertebral disks, and ligamentous structures) is described based on a specific set of imaging data, in either an explicit or idealized manner (Kallemeyn et al., 2010). However, such generic or individualized models do not allow for investigation of the full range of effects of variation in vertebral and spinal segment orientation and geometry and the influence of geometrical factors on the mechanical behavior of the spine (Laville et al., 2009).

Statistical shape modeling methods have been used to describe variability in the morphology of a population of anatomical structures in terms of a random field representation (Cootes et al., 1994; Lorenz and Krahnstover, 2000; Kaus et al., 2003; Rueckert et al., 2003). Current applications of statistical shape modeling (SSM) include automated image segmentation, image or object registration, object recognition, and disease diagnosis (Rueckert et al., 2003; Benameur et al., 2005; Dornaika and Ahlberg, 2006; Ferrarini et al., 2006; Rao et al., 2006; Koikkalainen et al., 2007). Statistical shape models capture the variability of biological structures by projecting a high dimensional representation of the structure onto a lower dimensional subspace of possible shapes constructed from a population of training shapes. A modeling approach combining SSM with finite element modeling allows investigation of loading behavior in specific individuals, as well as over the full range of morphological variability described within the model set.

The objectives of this study were to develop and implement methods based on SSM to describe the multivariate morphology and geometry within 3D imaging data for a set of cervical spines, to extend this model to a parametric finite element model of the cervical spine, and to quantitatively validate the performance of the parametric model against non-destructive experimental data. This study is the first step in this research program with the overall goal of investigating the effects of geometry and morphometry variation on the complex risks associated with cervical injury in scenarios that might be experienced in the occupational exposure of the modern warfighter, with the intent to provide the basis for mitigating these injuries in further work.

Methods

Image Processing

Five cadaver specimens representative of the 50th percentile male warfighter (based on weight) were obtained and the cervical (C3–T1) spine of each specimen was scanned using a computed tomography (CT) system (Aquilion 64, Toshiba – Medical Systems, Tokyo, Japan). CT data were filtered using a sequence of median and anisotropic diffusion filtering to reduce data noise. Filtered data were semi-automatically segmented to extract vertebra data from the CT image data (Figure 1) (Seg3D, The Center for Integrative Biomedical Computing, University of Utah, Salt Lake City, UT, USA). Watertight triangulated surfaces were generated to describe the outer boundary of each vertebra (e.g., five cervical spines × seven vertebrae) by computing the isosurface geometry for the segmented data region and smoothing the resulting surfaces to remove any stair-stepping effects due to out of plane image resolution (Figure 1) (MATLAB R2012a, The Mathworks, Inc., Natick, MA, USA). Vertebral surfaces were resampled, resulting in approximately 4,000 faces for each triangulated surface.

FIGURE 1

Figure 1. Image processing pathway from CT data to segmented data to vertebral surface.

Development of Individual Spine Models

All vertebral surfaces were positioned based on visual observation to yield nominal anatomic orientation of the full cervical spines and to correct for any positioning errors present during CT scanning (Scheer et al., 2013). The five surfaces at each vertebral level were registered to each other using an arbitrarily selected vertebral surface as the template. Vertices from the template surfaces were mapped onto the remaining four surfaces at each vertebral level and repositioned using a coherence point drift algorithm such that all vertices were positioned at corresponding anatomic locations between all vertebral surfaces at the same cervical level (Figure 2) (Myronenko and Song, 2010). Thus, the resulting vertebral surfaces were defined by the same surface mesh definition due to vertex correspondence across the set of vertebrae at each vertebral level. Average vertebrae were determined by averaging vertex positions for all vertebrae at each vertebral level.

FIGURE 2

Figure 2. A set of corresponding vertices is identified using blue, yellow, red, magenta, and orange dots on three different vertebral surfaces.

Volumetric tetrahedral meshes consisting of 25,000–35,000 elements (5,500–6,800 nodes) were defined for the average vertebral surfaces and refined to improve mesh quality (Tetgen, Weierstrass Institute for Applied Analysis and Stochastics, Berlin, Germany; Stellar, University of California, Berkeley, CA, USA). At each vertebral level, the average vertebral mesh was elastically warped to match each individual vertebra using displacement vectors calculated between corresponding surface vertices on the average vertebral mesh and each individual vertebral surface, resulting in a set of five corresponding vertebral mesh models for each cervical level (ANSYS v11.0, ANSYS, Inc., Canonsburg, PA, USA). Individual corresponding vertebral meshes were transformed back to the nominal anatomic position of the appropriate vertebra within each cervical spine.

In order to define intervertebral disk models, vertices bounding the approximate location of the intervertebral disk were selected on the adjacent vertebral endplates of the average cervical spine model (Figure 3). Due to mesh correspondence, vertices selected on the average motion segment models defined disk boundaries for all individual motion segments. Splines were fit through each set of selected surface vertices, resulting in a set of closed curves lying on the respective endplate surface and defining disk boundaries on adjacent endplates of the individual cervical spine models (Figure 4). Ellipses with dimensions proportional to the disk boundary curves were defined at the centroid of adjacent endplate curves and projected on to adjacent surfaces in order to define the interface between disk nucleus and annulus (Figure 4). Surfaces defining outer boundaries of the disk annulus and nucleus were generated between appropriate splines on adjacent vertebral endplates and disk endplate surfaces further bounded the intervertebral disk space. A hexahedral mesh of the intervertebral disks was defined within the bounded disk space for each individual spines (Truegrid, XYZ Scientific Applications, Inc., Livermore, CA, USA) (Figure 5).

FIGURE 3

Figure 3. Surface node sets were selected on adjacent vertebrae to specify the boundary of the intervertebral disk.

FIGURE 4

Figure 4. Splines were fit to the sets of nodes defining intervertebral disk boundaries on adjacent vertebrae and the resulting splines were projected to the vertebral endplates. Similarly, an ellipse was defined to represent the interface of the disk nucleus and annulus and projected to the adjacent vertebral endplates. Surfaces were defined by connecting the resulting closed curves and adjacent vertebral endplates further bounded the intervertebral disk space.

FIGURE 5

Figure 5. Typical intervertebral disk model with separate annulus and nucleus models.

Facet joint cartilage was modeled by projecting the set of surface triangles that define facet surfaces of adjacent vertebrae outward along the vertex normals to form a single layer of wedge elements. Facet cartilage elements were defined for each facet surface using a constant thickness, which was iteratively determined to maximize joint contact without facet surface interference. Ligaments (e.g., anterior longitudinal ligaments, posterior longitudinal ligaments, interspinous ligaments, ligamentum flavum, and intertransverse ligaments) and facet joint capsules were modeled using discrete spring elements to connect selected nodes on adjacent vertebrae. The facet regions and nodes representing ligament attachment sites were determined from the average cervical spine mesh. Mesh correspondence allowed the same sets of surface triangles and nodes to be used on each of the individual spine models.

Individual spines were translated such that the centroid of the central vertebra (e.g., C4) was located at the origin of the Cartesian space in order to register individual spines while maintaining individual spine intersegmental spacing and curvature. Each individual volumetric model of the spine included C3–T1 vertebrae, intervertebral disks, facet joints, and relevant ligamentous structures (Figure 6).

FIGURE 6

Figure 6. Average cervical spine model note: vertebral bodies are shown in gray, intervertebral disks are shown in blue, facet elements are shown in yellow, and ligaments are shown in red.

Development of Statistical Shape Models of the Cervical Spine

A SSM was generated to describe and investigate geometric variability in the set of five cervical spines. Joint point distribution models were constructed from all individualmeshes. The volumetric mesh for each individual was described by a shape parameter vector as:

${\mathit{\text{p}}}_{\mathit{\text{i}}}={\left(v_{1x},v_{1y},v_{1z},\dots ,v_{jx},v_{jy},v_{jz}\right)}^T$ (1)

where v_{j( xyz)} are the three-dimensional coordinates of the nodes in the volumetric spine model, j = 1, … , J = 54,960 nodes in the volumetric mesh, and i = 1, … , n = 5 denote each individual spine in the set.

The mean shape of all components (e.g., vertebrae, disks, facets, and ligament attachment points) in the set of cervical vertebrae was defined as:

$\overline{\mathit{\text{p}}}=\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\mathit{\text{p}}}_i$ (2)

and the correlation between individual models in the set was given by the empirical covariance matrix:

$\mathit{\text{S}}=\frac{1}{n}{\displaystyle \sum _{i=1}^n}\left({\mathit{\text{p}}}_i-\overline{\mathit{\text{p}}}\right){\left({\mathit{\text{p}}}_i-\overline{\mathit{\text{p}}}\right)}^T$ (3)

A principal components analysis (PCA) of the covariance matrix, S, results in a set of k = n − 1 eigenvalues (λ_k) and eigenvectors (q _k), which are the principal directions spanning a shape space centered at the mean, $\overline{\mathit{\text{p}}}$. The proportion of the total variance described along each eigenvector is equal to its corresponding eigenvalue divided by the sum of all eigenvalues; eigenvectors corresponding to the largest eigenvalues describe the majority of the variance. Thus, the finite element mesh for each cervical spine in the set were described in terms of the average model and a weighted linear combination of uncorrelated principal shape modes as:

${\mathit{\text{p}}}_v=\overline{\mathit{\text{p}}}+{\displaystyle \sum _{j=1}^m}{c}_j\sqrt{{\mathrm{\lambda }}_j}{\mathit{\text{q}}}_j$ (4)

where p_v is a vector containing coordinates for all nodes in the FE model, m is the number of eigenvalues, λ_j, and deviation from the average spine, $\overline{\mathit{\text{p}}}$, was determined as the sum of the products of a set of scalar weighting factors, c_j, and SSM standard deviations, $\sqrt{{\mathrm{\lambda }}_j},$ along the q_j (eigenvector) directions (Bredbenner et al., 2010; Nicolella and Bredbenner, 2012).

Accordingly, the highly correlated 3D spine geometry variables are reduced into a relatively small set of uncorrelated and independent composite morphological traits. All variability within the original set of spine models (originally described by over 164,880 variables) is now described by the weighting factors for four principal components for each cervical spine. Principal components are new descriptive variables that, by definition, are linear combinations of the original descriptive variables and, furthermore, all geometry information in the original models is retained in the new model descriptions.

In order to investigate the variability in vertebral morphology, intersegmental orientation, and overall curvature in the cervical spine models, a series of variation models were created and compared to the average models. Principal component weighting factors in Eq. 4 were modified to generate models describing the difference of 1.0 standard deviation of each principal component (i.e., shape mode) from the average model.

Development of a Parametric Finite Element Model of the Cervical Spine

The statistical shape model of the cervical spine (including vertebrae, disks, and ligaments) was generated in a form directly applicable to finite element analysis and geometry variation in the finite element model was explicitly described by principal component weighting factors. The average spine segment model was created. Additionally, the effects of spine geometry variation were investigated by modifying the average model with the geometry traits carried by the principal components (Eq. 4). Principal component weighting factors were defined as random variables with a mean, standard deviation, and distribution shape (Table 1). Vertebrae were modeled as rigid bodies, as this investigation did not consider vertebral fractures.

TABLE 1

Table 1. Random variable definitions.

The effects of variation and uncertainty in intervertebral disk and soft-tissue material properties were also investigated by considering appropriate material parameters as random variables. Each disk was modeled with a separate annulus and nucleus. Soft-tissue material properties were modeled based on experimental data found in the literature. Material behavior of the annulus was modeled with a transversely isotropic hyperelastic model with viscosity and material parameters were determined from experimental data provided by Lucas et al. (2006). Lucas et al. collected experimental data for the annulus using a ramp and hold loading protocol, with a high rate ramp (52 mm/s) to 0.88 mm and a 10-s hold period, and also under sinusoidal loading conditions (2 Hz with 0.65 mm peak-to-peak displacement). Material properties used in the present model were determined using the relaxation data and validated successfully against the dynamic loading data. The bulk modulus of the annulus was defined as a random variable (Table 1). A Prony series was defined to approximate the viscoelastic relaxation behavior of the annulus (Table 2). Material behavior of the nucleus was modeled using a fluid material model and material properties were determined from the literature (Table 3) (Teo and Ng, 2001; Nicolella et al., 2006). Material behavior for the ligamentous and capsular structures was defined using experimental force-displacement data that were collected under quasi-static loading conditions (Yoganandan et al., 2000) (Figure 7). As multiple spring elements were used to model anterior longitudinal ligaments, posterior longitudinal ligaments, interspinous ligaments, ligamentum flava, and joint capsules, experimental soft-tissue force-displacement data were scaled to account for distribution of the soft-tissue response over multiple discrete elements. Soft-tissue scale factors were treated as random variables (Table 1).

TABLE 2

Table 2. Viscoelastic properties for the intervertebral disk annulus.

TABLE 3

Table 3. Material properties for the intervertebral disk nucleus.

FIGURE 7

Figure 7. Experimental load-displacement data for soft-tissue material behavior.

The effects of variation in spine geometry and soft-tissue material properties on cervical spine kinematic response were determined by sampling the appropriate variable distributions 100 times using a Latin Hypercube approach within a probabilistic framework (NESSUS v8.0, Southwest Research Institute, San Antonio, TX, USA). Additionally, the contributions of each principal component to model loading response were investigated by creating the average model and models created by modifying the average model with the geometry traits described by +1 standard deviation of each principal component (Eq. 4), where principal components were considered individually. Disk and soft-tissue material properties were defined using mean values. In all model simulations, T1 was fixed and pure flexion-extension, left–right axial rotation, or left–right lateral bending moments of ±2.0 Nm were applied to C3 and models were solved using LS-DYNA v. 971 (LSTC, Livermore, CA, USA).

Validation of the Parametric Finite Element Model of the Cervical Spine

A hierarchical probabilistic verification and validation approach was used to quantify the performance of the parametric cervical spine model (ASME, 2006; Nicolella et al., 2006). Briefly, in previous work, material model parameters for the various ligaments and intervertebral disk components were determined and validated and the kinematic response of each motion segment (e.g., C3–C4, C4–C5, and C6–C7) was validated against independent experimental data without any alteration to the material model parameters. In the present work, the performance of the full C3–T1 parametric spine model was validated against independent experimental response data collected using cervical (C2–T1) specimens obtained from seven “normal” young human cadaver specimens (five males and two females; aged 33.4 ± 11.7 years with a range of 20–51 years) provided by Wheeldon et al. (2006). Ligamentous soft tissue was left intact in the experimental specimens. Wheeldon et al. fixed the T1 vertebra to a six-axis load cell, preconditioned the spine segment, and applied quasi-static pure moments of 0.33, 0.5, 1.0, 1.5, and 2.0 Nm to C2, resulting in cervical spine flexion-extension. Reaction load and vertebral rotation data were recorded throughout the experimental testing. Wheeldon et al. (unpublished data) used an identical loading protocol to load the cervical segments in left–right axial rotation and left–right lateral bending.

In order to evaluate the predictive performance of the parametric cervical spine model under non-destructive loading conditions, the mean and one-standard deviation load response envelopes were qualitatively compared between the model predictions and the experimental data for each loading model. In order to quantify model performance, we implemented a general performance metric that characterizes the disagreement between the model variation and relevant experimental data variation (Ferson et al., 2008; Francis et al., 2012). This quantitative metric provides a generalized approach to validation, rather than focusing on comparison of the mean predicted and experimental behaviors (Ferson et al., 2008). Empirical cumulative distribution functions (CDFs) were fit to the experimental data obtained from seven cervical spine specimens at each of five applied moment values (0.33, 0.5, 1.0, 1.5, and 2.0 Nm) over the loading range for each loading mode (e.g., flexion-extension, axial rotation, and lateral bending). Empirical CDFs were also fit to the predicted loading response for each loading mode at identical applied moment values, where response data were determined for each of the models generated by 100 Latin Hypercube samples of the variable space. Area metric values were determined as the total area difference between the CDFs for the experimental results and the model predictions at each of the five applied moment values for each loading mode and have units of degrees, as with the loading response. However, metric values are directly associated with the loading level; therefore, metric values were normalized by the mean experimental rotational displacement at each loading point. The resulting dimensionless values were used to measure the agreement between the response distributions over the full range of loading for each loading mode.

Results

The first principal component explains 63.8% of the total variability in cervical spine geometry, the first and second together explains 78.1%, and the first three principal components explain 90.3% of the geometric variability in the five cervical spine segments, with the remainder of the variability described by the fourth principal component.

Differences in morphology of individual vertebrae, intersegmental orientation, and overall spinal curvature are evident in qualitative comparisons between the average and variation models of the cervical spine (Figures 8–11). Variations in intersegmental orientation in each of the three major axes are clearly visible and, in some cases, obscure variations in specific vertebral morphology, although differences in posterior process shape and length is most evident. The cumulative effects of intersegmental variation over the cervical spine lead to overall variation of the full cervical spine segment.

FIGURE 8

Figure 8. Variation in geometry and alignment described by principal component (PC) 1.

FIGURE 9

Figure 9. Variation in geometry and alignment described by principal component (PC) 2.

FIGURE 10

Figure 10. Variation in geometry and alignment described by principal component (PC) 3.

FIGURE 11

Figure 11. Variation in geometry and alignment described by principal component (PC) 4.

As expected, investigation of the effects of spinal morphology variation on predicted loading response resulted in quantifiable variation in the loading responses between the average model and models created by combining the average model and the geometry variation described by +1.0 standard deviation of each principal component (Figures 12–14). Response variations were determined as percentage change with respect to the mean model. In the case of flexion-extension (Figure 12), principal component 4 had the largest effect on the flexion response (17.3%) and principal component 3 had the largest effect on extension (23.2%). Principal component 4 had the largest effect on axial rotation (Figure 13), with 19.5 and 25.5% variation for right and left rotation, respectively. In the case of lateral bending (Figure 14), principal component 3 had the largest effect on right bending (22.4%) and principal component 2 had the largest effect on left bending (27.1%).

FIGURE 12

Figure 12. Effects of geometry variation on flexion-extension.

FIGURE 13

Figure 13. Effects of geometry variation on axial rotation.

FIGURE 14

Figure 14. Effects of geometry variation on lateral bending.

Qualitatively, the mean and variation in the probabilistic results closely match those of the experimental results. In general (Figures 15–17), the variability of the predicted response is greater than that of the experimental response; however, the mean predicted response lies within the +1.0 standard deviation envelope for experimental data in all cases. The parametric model also predicts greater variation than the experimental results as the applied moment increases for each loading mode.

FIGURE 15

Figure 15. Kinematic response and normalized CDF area metric for flexion-extension.

FIGURE 16

Figure 16. Kinematic response and normalized CDF area metric for left-right axial rotation.

FIGURE 17

Figure 17. Kinematic response and normalized CDF area metric for left-right lateral bending.

Disagreement between model and experimental CDFs were quantified using the normalized CDF area metric as between 0.22 and 0.45 for each time point in extension (Figure 15), between 0.03 and 0.19 for right axial rotation (Figure 16), and ranges between 0.13 and 0.24 for all time points for all other loading modes (Figures 15–17).

Discussion

This study demonstrated that significant variability was present in the geometry of a small group of cervical spine segments, both in terms of vertebral morphometry and in the intervertebral orientation and overall spinal curvature. Furthermore, SSM is capable of efficiently describing variability in the complex vertebral morphometry, intersegmental orientation, and overall spinal curvature and demonstrates the complex relationship between predicted response and variation in model input parameters. Although we have investigated conditions within the non-destructive loading range of the cervical spine, we suggest that, with additional validation under destructive loading conditions, the present implementation of parametric SSM-based finite element analysis methods is applicable to investigations of cervical spine injury mechanisms. These loading scenarios include those that might be encountered during high-speed impact in a moving vehicle, as the result of an explosive blast, or the ejection of a fighter pilot. More importantly, the parametric high-fidelity description of variability in spinal morphology, along with the ability to vary relevant material properties within a probabilistic framework allows the investigation of the effects of body size, position, mass distribution, and other relevant factors on injury risk prediction.

Intersegmental orientation and spinal curvature is highly variable between individuals and has a substantial role in individual disposition to neck and back injuries, including muscular fatigue or soft-tissue injuries, as well as more serious injuries related to hyperflexion and hyperextension (Coakwell et al., 2004; Frechede et al., 2006). Variation in intersegmental orientation and overall spinal curvature are modeled implicitly in the current SSM approach; however, explicitly modeling vertebral orientation and spinal curvature variables would allow the investigation of the role and interaction of vertebral morphology, intersegmental orientation, and spinal curvature variables on the likelihood of injury under various loading conditions.

The small sample size of spines is a limitation of the present study and it remains to be seen whether the range of morphological differences observed within the study sample are sufficient to describe a larger population of warfighters. In this study, we choose to focus on the development of a model of the 50th percentile male as a starting point, with the intent to develop methodology somewhat representative of the loading response of an “average” male warfighter. The geometric and material property variation described in the cervical spine model may well be different than that within the experimental set of spines; however, probabilistic methods allowed investigation of the effects of variability in morphology and material properties on kinematic response. Based on qualitative and quantitative comparison of the model simulations to experimental loading response data that is currently available, we suggest that the model performs adequately under relatively low-level loading conditions in multiple loading directions. In ongoing work, the cervical spine model is being exercised and validated against destructive loading modeling blast conditions for the same set of spines used to construct the model. It remains to be determined whether differences in vertebral morphology and spine geometry and the ability to describe these differences using a small set of principal shape modes will be capable of future risk classification under the widely varying scenarios that may be encountered by a warfighter. We also note ligament material properties were defined based on quasi-static loading conditions and this may present a limitation in evaluation of the model under destructive, high rate loading conditions.

In earlier work involving a parametric cervical spine model based on idealized vertebral geometry, we have demonstrated that both morphological differences between small female, large female, small male, and large male groups (determined by body mass) and material properties of ligaments and the intervertebral disks have a substantial effect on the predicted kinematic response under loading modes (e.g., flexion-extension, axial rotation, and lateral bending) with identical applied moments (Nicolella et al., 2006). Other groups using idealized parametric models of models with subject-specific geometry have also found that geometry and orientation and material models employed in describing soft-tissue behavior in both cervical spine motion segments lead to differences in the loading response (Kumaresan et al., 1999; del Palomar et al., 2008; Laville et al., 2009; Kallemeyn et al., 2010). In other work, finite element models of the cervical spine were deformed to model gross global changes in curvature (e.g., lordosis, straight, and kyphosis models) based on specified Cobb angle (Frechede et al., 2006). Distributions of strains, forces, and moments along the cervical spine were found to be dependent on loading condition (e.g., vehicular rear-end, frontal, lateral, and oblique impact) and global cervical curvature affected the magnitude of strains, forces, and moments experienced by the spine in all loading directions. However, we are unaware of probabilistic investigations where SSM methods were combined with finite element analyses to systematically investigate the effects of the continuous normal geometric variation in a sample of spines and material property uncertainty on loading response in the cervical spine. We suggest that the ability to continuously vary vertebral morphology, intersegmental orientation, and spinal curvature along with soft-tissue material properties within realistic parameter spaces strongly enhances the current investigative tools employed to understand and mitigate the risk of fracture and soft-tissue injuries in military combat scenarios.

We note that in validation of the present model under low-loading conditions in multiple directions, uncertainty and variability in both experimental results and model inputs (and, therefore, model results) are considered stochastically through the determination of a empirical distributions to describe experimental results and the use of probabilistic modeling methods to determine computation results. Therefore, we have quantitatively compared the statistical distributions resulting from both simulation predictions and experimental observations (Liu et al., 2011). The CDF area metric employed here to quantify disagreement between model predictions and experimental data evaluates the ability of the model to not only predict the mean, but the amount of variation in the experimental response. Additionally, the metric prevents incorrect “validation” of model predictions based on experimental data with large corridors of uncertainty. At each time point for each loading mode, the CDF area metric is a quantitative measure of agreement between model predictions and experimental data with reliance on expert opinion to evaluate the accuracy of the predictions. This is relatively straightforward to understand since the metric can be displayed graphically and has an intuitive meaning. Furthermore, unlike other computational model validation approaches, the area metric incorporates variability in both the experiment and model. It has a minimum value of 0, indicating the experimental and model CDF are the same, and has an infinite maximum value with increasing values indicating an increasingly poor match. We note that the quantitative error metric provides relative error (rather than absolute error) within the range of data analyzed and is, of course, dependent on whether the experimental data used for comparison is representative of the conditions modeled. As such, the area-based metric is limited in the case where the experimental data are insufficient (Liu et al., 2011). We have not exercised the model to produce injury in this study, and therefore, have not fully explored the kinematic and dynamic range of the model and Wheldon et al. have not thoroughly exercised the experimental specimens (Wheeldon et al., 2006; unpublished data). Accordingly, we suggest that some concern regarding insufficiency of experimental data may be reduced.

A potential use of a quantitative metric, such as the one used in the current analysis, is to investigate the effect of the use of alternative sources of data for material constitutive modeling, the incorporation of additional validation data when available, and the comparison of alternative or competing models and modeling approaches. A limitation of the area validation metric, however, is that an established standard for what constitutes an “acceptable” validation threshold has not been established. Therefore, this and other validation metrics should be used within the context of the overall objectives and goals of the use of the computational model. Qualitative comparison of the mean and standard deviation envelopes of the experimental and simulated loading responses suggest that the amount of error realized in the simulations are acceptable for the applied loading conditions. The quantitative metric values in each loading mode suggest that error is greater near the so-called “slack region” close to an applied moment of 0 Nm, where the model is not explicitly accounting for the ligamentous laxity that exists within cadaveric specimens. We note that in subsequent investigations utilizing this parametric cervical spine model with injurious loading conditions, we will be less concerned with areas of laxity during low-level loading and more concerned with the peak regions of the loading range; however, it remains to be seen whether a similar pattern of error will be present over a larger loading range.

In conclusion, SSM provides a means of explicitly describing complete spinal morphology and geometry and allows the complex spatial variation in intersegmental orientation to be statistically investigated between individuals. SSM methods coupled with finite element analyses within a probabilistic framework, along with the ability to quantitatively validate the overall model performance, provides innovative and important steps toward describing differences in vertebral morphology, intersegmental orientation, and spinal position and curvature, as well as important variation in material properties, which may directly lead to understanding and mitigating the risks of soft-tissue and hard tissue injury in the spine and other musculoskeletal structures.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to express their appreciation to the Department of Defense Congressionally Directed Medical Research Programs (CDMRP) under contract W81XVVH-09-2-0168. We would also acknowledge the support of the Naval Air Warfare Center Aircraft Division (NAWCAD) (Dr. Barry Shender, program manager) for development of a portion of the methods utilized in this project. The authors would also like to acknowledge the assistance of Richard L. Debellevue in segmenting QCT data.

References

Anderson, H. T. (1988). Neck injury sustained during exposure to high-G forces in the F16B. Aviat. Space Environ. Med. 59, 356–358.

Pubmed Abstract | Pubmed Full Text | Google Scholar

ASME. (2006). Guide for Verification and Validation in Computational Solid Mechanics. New York, NY: American Society of Mechanical Engineers.

Google Scholar

Belmont, P. J. Jr., McCriskin, B. J., Sieg, R. N., Burks, R., and Schoenfeld, A. J. (2012). Combat wounds in Iraq and Afghanistan from 2005 to 2009. J. Trauma Acute Care Surg. 73, 3–12. doi:10.1097/TA.0b013e318250bfb4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Benameur, S., Mignotte, M., Labelle, H., and De Guise, J. A. (2005). A hierarchical statistical modeling approach for the unsupervised 3-D biplanar reconstruction of the scoliotic spine. IEEE Trans. Biomed. Eng. 52, 2041–2057. doi:10.1109/TBME.2005.857665

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bredbenner, T. L., Eliason, T. D., Potter, R. S., Mason, R. L., Havill, L. M., and Nicolella, D. P. (2010). Statistical shape modeling describes variation in tibia and femur surface geometry between control and incidence groups from the osteoarthritis initiative database. J. Biomech. 43, 1780–1786. doi:10.1016/j.jbiomech.2010.02.015

CrossRef Full Text | Google Scholar

Coakwell, M. R., Bloswick, D. S., and Moser, R. Jr. (2004). High-risk head and neck movements at high G and interventions to reduce associated neck injury. Aviat. Space Environ. Med. 75, 68–80.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Cootes, T. F., Hill, A., Taylor, C. J., and Haslam, J. (1994). Use of active shape models for locating structure in medical images. Image Vis. Comput. 12, 355–365. doi:10.1016/0262-8856(94)90060-4

CrossRef Full Text | Google Scholar

del Palomar, A. P., Calvo, B., and Doblare, M. (2008). An accurate finite element model of the cervical spine under quasi-static loading. J. Biomech. 41, 523–531. doi:10.1016/j.jbiomech.2007.10.012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Dornaika, F., and Ahlberg, J. (2006). Fitting 3D face models for tracking and active appearance model training. Image Vis. Comput. 24, 1010–1024. doi:10.1016/j.imavis.2006.02.025

CrossRef Full Text | Google Scholar

Ferrarini, L., Palm, W. M., Olofsen, H., van Buchem, M. A., Reiber, J. H., and dmiraal-Behloul, F. A. (2006). Shape differences of the brain ventricles in Alzheimer’s disease. Neuroimage 32, 1060–1069. doi:10.1016/j.neuroimage.2006.05.048

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ferson, S., Oberkampf, W. L., and Ginzburg, L. (2008). Model validation and predictive capability for the thermal challenge problem. Comput. Methods Appl. Mech. Engin. 197, 2408–2430. doi:10.1016/j.cma.2007.07.030

CrossRef Full Text | Google Scholar

Francis, W. L., Eliason, T. D., Thacker, B. H., Paskoff, G. R., Shender, B. S., and Nicolella, D. P. (2012). Implementation and validation of probabilistic models of the anterior longitudinal ligament and posterior longitudinal ligament of the cervical spine. Comput. Methods Biomech. Biomed. Engin. 17, 905–916. doi:10.1080/10255842.2012.726353

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Frechede, B., Bertholon, N., Saillant, G., Lavaste, F., and Skalli, W. (2006). Finite element model of the human neck during omni-directional impacts. Part II: relation between cervical curvature and risk of injury. Comput. Methods Biomech. Biomed. Engin. 9, 379–386. doi:10.1080/10255840600980940

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kallemeyn, N., Gandhi, A., Kode, S., Shivanna, K., Smucker, J., and Grosland, N. (2010). Validation of a C2-C7 cervical spine model using specimen-specific flexibility data. Med. Eng. Phys. 32, 482–489. doi:10.1016/j.medengphy.2010.03.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kaus, M. R., Pekar, V., Lorenz, C., Truyen, R., Lobregt, S., and Weese, J. (2003). Automated 3-D PDM construction from segmented images using deformable models. IEEE Trans. Med. Imaging 22, 1005–1013. doi:10.1109/TMI.2003.815864

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Koikkalainen, J., Hirvonen, J., Nyman, M., Lotjonen, J., Hietala, J., and Ruotsalainen, U. (2007). Shape variability of the human striatum – Effects of age and gender. Neuroimage 34, 85–93. doi:10.1016/j.neuroimage.2006.08.039

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kumaresan, S., Yoganandan, N., Pintar, F. A., and Maiman, D. J. (1999). Finite element modeling of the cervical spine: role of intervertebral disc under axial and eccentric loads. Med. Eng. Phys. 21, 689–700. doi:10.1016/S1350-4533(00)00002-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Laville, A., Laporte, S., and Skalli, W. (2009). Parametric and subject-specific finite element modeling of the lower cervical spine. Influence of geometrical parameters on the motion patterns. J. Biomech. 42, 1409–1415. doi:10.1016/j.jbiomech.2009.04.007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Liu, Y., Chen, W., Arendt, P., and Huang, H.-Z. (2011). Toward a better understanding of model validation metrics. J. Mech. Des. 133, 71005. doi:10.1115/1.4004223

CrossRef Full Text | Google Scholar

Lorenz, C., and Krahnstover, N. (2000). Generation of point-based 3D statistical shape models for anatomical objects. Comput. Vis. Image Underst. 77, 175–191. doi:10.1006/cviu.1999.0814

CrossRef Full Text | Google Scholar

Lucas, S., Bass, C., Salzar, R., Shender, B. S., and Paskoff, G. (2006). High-rate viscoelastic properties of human cervical spinal intervertebral discs. Annual Meeting of the American Society of Biomechanics, American Society of Biomechanics, Blacksburg, VA.

Google Scholar

Milby, A. H., Halpern, C. H., Guo, W., and Stein, S. C. (2008). Prevalence of cervical spinal injury in trauma. Neurosurg. Focus 25, E10. doi:10.3171/FOC.2008.25.11.E10

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Myronenko, A., and Song, X. (2010). Point-set registration: coherent point drift. IEEE Trans. Pattern Anal. Mach. Intell. 32, 2262–2275. doi:10.1109/TPAMI.2010.46

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Nicolella, D. P., and Bredbenner, T. L. (2012). Development of a parametric finite element model of the proximal femur using statistical shape and density modelling. Comput. Methods Biomech. Biomed. Engin. 15, 101–110. doi:10.1080/10255842.2010.515984

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Nicolella, D. P., Francis, W. L., Bonivtch, A. R., Thacker, B. H., Paskoff, G. R., and Shender, B. S. (2006). Development, Verification, and Validation of a Parametric Cervical Spine Injury Prediction Model, 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. Newport, RI: American Institute of Aeronautics and Astronautics.

Google Scholar

Owens, B. D., Kragh, J. F. Jr., Wenke, J. C., Macaitis, J., Wade, C. E., and Holcomb, J. B. (2008). Combat wounds in operation Iraqi freedom and operation enduring freedom. J. Trauma 64, 295–299. doi:10.1097/TA.0b013e318163b875

CrossRef Full Text | Google Scholar

Rao, A., Babalola, K., and Rueckert, D. (2006). Canonical Correlation Analysis of Sub-Cortical Brain Structures Using Non-Rigid Registration, Biomedical Image Registration. Berlin: Springer, 66–74.

Google Scholar

Rueckert, D., Frangi, A. F., and Schnabel, J. A. (2003). Automatic construction of 3-D statistical deformation models of the brain using nonrigid registration. IEEE Trans. Med. Imaging 22, 1014–1025. doi:10.1109/TMI.2003.815865

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Schall, D. G. (1989). Non-ejection cervical spine injuries due to +Gz in high performance aircraft. Aviat. Space Environ. Med. 60, 445–456.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Scheer, J. K., Tang, J. A., Smith, J. S., Acosta, F. L. Jr., Protopsaltis, T. S., Blondel, B., et al. (2013). Study, cervical spine alignment, sagittal deformity, and clinical implications: a review. J. Neurosurg. Spine 19, 141–159. doi:10.3171/2013.4.SPINE12838

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Teo, E. C., and Ng, H. W. (2001). Evaluation of the role of ligaments, facets and disc nucleus in lower cervical spine under compression and sagittal moments using finite element method. Med. Eng. Phys. 23, 155–164. doi:10.1016/S1350-4533(01)00036-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wade, A. L., Dye, J. L., Mohrle, C. R., and Galarneau, M. R. (2007). Head, face, and neck injuries during Operation Iraqi Freedom II: results from the US Navy-Marine Corps Combat Trauma Registry. J. Trauma 63, 836–840. doi:10.1097/01.ta.0000251453.54663.66

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wheeldon, J. A., Pintar, F. A., Knowles, S., and Yoganandan, N. (2006). Experimental flexion/extension data corridors for validation of finite element models of the young, normal cervical spine. J. Biomech. 39, 375–380. doi:10.1016/j.jbiomech.2004.11.014

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Yoganandan, N., Kumaresan, S., and Pintar, F. (2000). Geometric and mechanical properties of human cervical spine ligaments. J. Biomech. Eng. 122, 623–629. doi:10.1115/1.1322034

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar