Cervical Spine Finite Element Models for Healthy Subjects: Development and Validation

The cervical spine protects the spinal cord and nerve roots and provides structural support and flexibility of movement for the head and neck complex. Injuries to the cervical spine can be fatal [1]. Finite element (FE) modeling can help the scientific community understand injury mechanisms and dysfunctions of the head and neck complex [2,3].

A virtual representation of a physical object or process to its real-time digital counterpart is known as a digital twin [4]. With drastically improved computational power, simulation techniques now can solve more complicated problems by the models of different designs and accurately predict the behavior of the systems [5,6]. Recent developments in numerical simulation, artificial intelligence, and signal processing brought the application of digital twins in human biomechanics closer to reality. Real-time monitoring and analysis with digital twins of spine has become one of the most cutting-edge technologies in the field of biomechanics and medical science [7]. The first step to make a digital twin is to create a physics-based digital human model. In spine research, a wide variety of researchers use computed tomography (CT)-scan-based finite element models for this step [7].

There has been a substantial amount of work in modeling the cervical spine using finite element techniques in the past few decades [8–11]. Very early models used lumped parameter head–neck modeling to predict the response of the human head–neck system [12]. The geometries of these models were defined by simple rigid masses to be vertebrae and beam and spring elements that were used to define ligaments, intervertebral discs, and muscles [13–15]. These models did not have the biomechanically relevant geometric boundaries of the cervical spine parts and predicted results deviated from actual conditions [16]. Later models were more focused on finite element techniques and geometric boundaries of the cervical vertebrae [9,17–20]. A suitable source of these geometries was CT scans. Yoganandan et al. [9,18–20] developed a model of the lower cervical spine with C4–C6 vertebrae using CT scans. The model was validated under axial compression against experimental data. Goel and Clausen [21] developed an FE model of the C5–C6 spinal unit and used the model to predict intervertebral disc pressure, ligament tensions, and facet joint forces. This model adopted nonlinear properties for the ligaments and was validated under various loading modes (axial compression, axial rotation, and lateral bending). None of the above-mentioned models represented the full cervical spine as they had only one or two spinal units. These models lacked the physical aspects of a full or multilevel cervical spine. The necessity of a multilevel cervical spine model was defined in Ref. [16] and several full cervical spine models were developed [22–24]. Stemper et al. [23] built a full cervical spine model including the skull, C1–T1 vertebrae, intervertebral discs, and relevant ligaments. They studied the material sensitivity of motor vehicle inertial impact by applying loads to the first thoracic vertebra (T1). Meyer et al. [24] developed a complete head and neck finite element model and validated it against in-house experimental data.

Recent FE models of the cervical spine were more detailed and biomechanically accurate [10,11,25–39]. Zhang et al. [31] developed a nonlinear C0–C7 head and neck model and compared the range of motion (ROM) in different loading conditions against the experimental results in Panjabi [2]. Wheeldon et al. [29] validated a finite element model of the lower cervical spine (C4–C7) to study flexion/extension, left and right lateral bending, and left and right axial rotation. The model adopted a three-zone disc annulus to represent different material properties at the anterior, lateral, and posterior annulus. Kallemeyn et al. [25] developed a finite element model that included the C2–C7 vertebrae. The model was validated using in-house experimental data for a load range of ±1 N m. Panzer et al. [10] studied the cervical spine response in a frontal crash using a detailed finite element model of the head and neck system. The model was validated for flexion/extension and studied the cervical spine response for frontal impact for 8G to 22G severity. Erbulut et al. [11] developed a finite element model of C2–T1 vertebrae to study the effect of soft tissues on stability. Preclinical evaluation with patient-specific modeling is now possible with these image-based finite element models. With the advancement in the field of constraint-based computer-aided design (CAD) [40], artificial intervertebral discs can now be tested simultaneously with the finite element models of the cervical spine. Disc arthroplasty (artificial disc replacement) is required for patients to restore the mobility of the spine by replacing damaged disc parts. An early assessment with a finite element model and computer-designed disc models can be of importance to evaluate pre- and post-operative states [41,42]. More recent models target a specific population. For example, Finley et al. [39] developed a model specific for pediatric subjects. Models develop by Herron et al. [37] are specific to female subjects.

Most models in the literature developed finite element models based on 50th percentile male subjects. Lack in the variety of female cervical spine finite element models affects capturing gender-specific differences in biomechanical response [34]. The probability of chronic neck pain in female patients is twice as high than male patients [43]. Being smaller in size and 20–30% weaker in flexion–extension [44], it is not recommended to scale a male cervical spine model to represent a female subject [34]. Östh et al. [34] created a cervical spine model from CT scans of a 26-year-old female subject. The model specifically aimed to validate a female finite element cervical spine model with biomechanically relevant ligaments. Herron et al. [37] developed three individual cervical spine models from three different female subjects. Though the models were developed with the same material properties, results showed variation in the range of motion and disc pressure between subjects.

The finite element models described above, in general, adopted geometries from one subject and implemented material properties from the literature that may not match the anthropometry, sex, or age of a specific subject. Though these models were consistent with experimental results, one can argue about the compatibility of the model to predict the intervertebral motion of other subjects due to intersubject variation. Dreischarf et al. [45] compared eight different finite element models of the lumbar spine and showed that the combined effect of the pooled median of individual model results improves the prediction of spinal motion. Xu et al. [46,47] implemented this idea and built five finite element models from five different healthy subjects to show an improved prediction method.

This research aims to investigate this intersubject variability for the cervical spine by creating five finite element models from five different female subjects based on their CT scans. These models include C2–C7 vertebrae, intervertebral discs, endplates, and five major ligaments of the cervical spine. Mesh convergence and material sensitivity tests are performed to ensure the efficiency of the finite element models. The simulation results from these five models are validated against experimental results as well as against finite element studies in the literature for flexion/extension, lateral bending, axial rotation, and intervertebral disc pressure [3,6,8,25,31,48–52]. Intersubject variability among these models is studied to improve the prediction of different subjects.

The geometries of the neck finite element models are developed based on CT scans of five female subjects (Table 1). These scans provide the basic shape of six cervical vertebrae (C2—C7). The three-dimensional vertebral geometries are created by importing the two-dimensional CT scan images to 3D Slicer [53]. The segmentation is done manually by drawing contours at the boundary of each slice of the CT scan. Contours from all three directions create a vertebral 3D surface geometry. Slice thickness plays an important role during segmentation. Data used in this research have an average thickness of 0.3 mm between each slice. The higher the thickness between each slice, the lower the resolution of the geometry will be. The converted 3D geometries consist of irregularities like spikes and holes. Geometric surfaces of the vertebrae go through a smoothing process using Blender^® (Blender Foundation, Amsterdam, Netherlands). This procedure removes all spikes and unwanted holes on the surface of the vertebral geometry. Meshing is performed in IA-FEMesh (University of Iowa, Iowa City, IA) [54] using a multiblock mesh technique [26]. Three-dimensional blocks surrounding each vertebra are used to determine the number of divisions for hexahedral eight-node meshes. Mesh size is varied to keep the geometry as original as possible because the low-resolution mesh can end up providing a distorted geometry. The vertebral surface, endplates, and intervertebral discs are modeled with the same number of concentric elements at the contact surfaces. The steps involved to build the vertebral model are shown in Fig. 1.

The vertebral bodies of the model are made up of the cortical bone, the cancellous core, and the posterior bony structure. The cortical bone is modeled as 2D shell elements that surround the cancellous bone whereas the cancellous and posterior bone are modeled as 3D eight-node elements [55] using ls-dyna prepost (Livermore Software Technology Corporation, Livermore, CA). The cortical shells had a constant thickness of 0.5 mm [34]. The intervertebral discs are modeled between two vertebrae after meshing using ls-dyna prepost. The neighboring vertebrae surfaces of each segment are merged to create a complete disc model. An intervertebral disc consists of a nucleus at the center, endplates, an annulus ground that are all other elements of the disc that do not contain the nucleus or endplates and an annulus fibrosus. These discs are modeled with concentric layers of elements where the nucleus remains centered. The nucleus is around 40% of the entire disc volume [11]. The annulus fibrosus is modeled as concentric 2D shell elements between annulus ground element layers. An endplate of approximately 1 mm thickness is integrated between each disc and vertebra with 3D eight-node elements [18,56]. Figure 2 shows the intervertebral disc components modeled.

Ligaments in the cervical spine provide stability and flexibility of motion [57]. Five major ligaments (anterior longitudinal ligament: ALL, posterior longitudinal ligament: PLL, ligamentum flavum: LF, interspinous ligament: ISL, and capsular ligament: CL) are incorporated in the model as found across literature [3,10,11,25,48]. Two-node cable elements are used to model the ligaments that only allow a tensile load. The fixed cross-sectional area of each element of ligaments is taken from the literature [25]. The ALL and PLL are at the anterior and posterior of the vertebral body. The ISL connects the spinous processes. The LF is connected to the internal side of the posterior bones. Figure 3 shows all five cervical spine finite element models with all components.

The material properties are mostly adapted from the experiment conducted by Kallemeyn et al. [25] and compared with the computational study by Erbulut et al. [11]. Vertebrae parts (cortical, cancellous, and posterior bone) and endplates are modeled as an isotropic linear elastic material and the annulus fibrosus as concentric 2D sheet elements between annulus grounds. The nucleus is considered incompressible fluid elements [25]. Though vertebrae and intervertebral disc parts can be modeled with linear elastic materials, ligaments show highly nonlinear behavior [10,25]. Nonlinear stress versus strain curves reported by Kallemeyn et al. [25] are adopted to model the ligaments with an appropriate cross-sectional area as shown in Fig. 4. Detailed material properties are listed in Table 2.

The models developed are set to be validated in three directions of motions: flexion/extension, lateral bending, and axial rotation. A pure moment is applied at the top surface of the odontoid process of the C2 vertebra ranging from 0 to 2 N m for flexion/extension. For lateral bending and axial rotation, the load is fixed to 1 N m. The nodes at the bottom surface of the C7 vertebra are fixed in all degrees of freedom. Note that no compressive follower loads are applied due to its low contribution to intervertebral rotation [50].

The computational analysis is carried out through the commercial software ls-dyna (Livermore Software Technology Corporation, Livermore, CA) using the High-Performance Computing Center (HPCC) (server: Quanah, benchmark performance: 485 TFLOPS, total cores: 16,812) at Texas Tech University.

A mesh convergence test is necessary to justify the selected number of elements in a finite element model. This process ensures an accurate result and also a reduction in simulation time [46]. To reduce the complexity of the process, only one model (subject 1) is being tested here for mesh convergence. Three mesh resolutions of the full model are created (coarse, medium, and fine) for mesh sensitivity analysis. The fine mesh has approximately 2.60 times the number of elements compared with the medium mesh shown in Table 3. Figure 5 shows the models of three mesh resolutions. The three mesh resolutions are simulated under flexion with a moment of 1 N m. Intervertebral rotation is taken as the target parameter for mesh convergence analysis [58,59]. The medium mesh is considered converged if the percentage difference of the specific parameter (in this case, intervertebral rotation) is less than 10% from the fine mesh [58,60–63].

The five finite element models’ mean and standard deviation are compared with experimental results and validated finite element models in the literature [3,6,8,25,31,48–52]. In the validation, the models are tested for different loading conditions with pure moments applied at the top surface of the C2 vertebra shown in Table 4: flexion/extension, lateral bending, and axial rotation. Intervertebral rotation and intervertebral disc pressure are considered parameters for validation.

The percentage differences in angular displacement at each vertebral level between coarse and medium mesh and between medium and fine mesh are shown in Fig. 6 for subject 1. The differences between the coarse and medium mesh are higher than those between medium and fine mesh. Level C2–C3 is found to be most sensitive to mesh resolution with a maximum difference of 29.23%. As for the differences between the medium and fine mesh, all the vertebral levels are found to be below 10%. Therefore, medium mesh resolution is considered as converged and the same resolution is applied to the rest of the models [58,60–63].

The models are validated under different loading conditions for flexion/extension, lateral bending, and axial rotation against the literature. The load is varied between 0 and 2 N m with 0.5 N m increments for flexion and extension. At each load, the intervertebral rotation is recorded for all the vertebral levels (C2–C3, C3–C4, C4–C5, C5–C6, and C6–C7) [3,10,48,52]. As for lateral bending and axial rotation, only a fixed load of ±1 N m is applied to keep it consistent with the literature [8,25,31,49]. The full range of motion (left and right) is recorded for all the vertebral levels. Throughout the following sections, “current study” represents the average property between all five subjects.

The intervertebral rotations under pure moment are plotted in Figs. 7 and 8 along with selected experimental and finite element studies in the literature [3,10,48,52]. The results show asymmetry in flexion (+ve) and extension (−ve) which is consistent with the literature. With the increase in pure moment, the angle of rotation increases, thus increasing the ROM. The combined results for flexion/extension are within the experimental corridor in Fig. 7 and match well with other finite element model results in the literature in Fig. 8.

Figure 9 shows the ROM for the left and right lateral bending at each vertebral level under a pure moment of ±1 N m. Here the ROM represents the sum of left and right lateral bending angular displacements. Both C2–C3 and C3–C4 show a comparatively smaller ROM than other intervertebral levels. The current model is overall stiffer in lateral bending compared to the results of Panjabi et al. [8] and Kallemeyn et al. [25]. The simulated ROM mostly agrees with the results of Zhang et al. [31] and Traynelis et al. [49].

Figure 10 shows the ROM at different intervertebral levels compared to the results in the literature [8,25,31,49]. The models predict the smallest ROM at the C3–C4 level and the largest at the C6–C7 level. The results generally agree with the results in the literature except the studies in Kallemeyn et al. [25] which are more flexible than most studies in the literature.

Intervertebral discs are important for carrying loads and absorbing shocks. Any difference in intradiscal pressure affects the nutrition flow of the disc, thus affecting its shock-absorbing ability [51]. In this section, simulated intradiscal pressures of the models are compared to existing experimental results [51] and one finite element simulation result in the literature [50]. The pressure of five models is represented as mean and standard deviation in Fig. 11. Though our models predict all five levels of cervical spine (C2–C3, C3–C4, C4–C5, C5–C6, C6–C7), the experiment in the mentioned literature contains only C3–C4 and C5–C6. Model-predicted results agree more on extension than flexion with the experiment. Overall intradiscal pressures vary greatly between subjects, as seen in Fig. 11, with very high standard deviations.

Figures 12–14 show the ROM for flexion/extension, lateral bending, and axial rotation for all five subjects. For flexion/extension, the ROM is under 2 N m load and for lateral bending and axial rotation, the load is 1 N m only.

Mesh convergence tests are necessary to ensure that the results do not depend on the mesh resolution and that the simulation has a feasible computation time. Though this study considers intervertebral rotation as a target parameter for mesh convergence, other properties like disc pressure, strain energy, and facet joint force can also be metrics of mesh sensitivity [58]. The model converged at around 100k similar to the finite element study by Panzer et al. [10].

The models are validated in flexion/extension against three published results among which two are experimental studies and one is a finite element simulation study. The intervertebral rotation angle at each vertebral level shows asymmetry around the center which agrees with experimental results by Wheeldon et al. [48] and Nightingale et al. [3,52]. As shown in Fig. 7, the mean intervertebral rotations are within the experimental corridor for extension. The flexion is stiffer compared to the experiments for C2–C3, C3–C4, C4–C5, and C5–C6, which agrees with the finite element study of the cervical spine by Panzer et al. [10] shown in Fig. 8. The reason behind this inconsistency may originate from the facet joint orientations [11,48] as these joints are more involved for extension than flexion. Facet joints provide restriction between vertebrae during extension and geometries of these facet joints vary between subjects. One important factor here to consider is that all the models developed in this study are from female subjects whereas experimental data are dominated by male subjects. The subjects for the experiment by Nightingale et al. [3,52] are all male (16 male subjects) and the experiment by Wheeldon et al. [48] has a majority of male subjects (five male, two female). Previous study results show that female cervical curvature is significantly lower compared to a male’s [64]. This variation leads to different neutral positions between males and females and therefore can affect bending stiffness between subjects [10].

For lateral bending validation, the smallest mean for ROM is found at the C2–C3 level which agrees with results in Zhang et al. [28]. For all vertebral levels, the results are mostly close to those in Zhang et al. [31]. For axial rotation, the ROM closely follows the results of Panjabi et al. [8] for C3–C4, C4–C5, and C5–C6. The smallest range of motion is found at the C3–C4 level and the largest at the C2–C3 level similar to Zhang et al. [31]. Panjabi et al. [8] showed the smallest ROM at the C6–C7 level and the largest at the C4–C5 level.

With limited experimental data, validation for intradiscal pressure was limited to C3–C4 and C5–C6 only. Pospiech et al. [51] performed experiments based on cadaveric specimens with cables connected to vertebrae to simulate muscles. Though models in this study do not contain muscles, the results agree well for extension. The differences in pressure for flexion may originate from the fact that during the experiment, the specimen was preloaded with a 10 N axial preload to keep the levels together. The upper body load is much lower for cervical spine than lumbar spine with cervical spine finite element models in the literature usually avoiding preloading. Cai et al.’s [50] results show a follower load and have very little effect on ROM and high effect on disc pressure for the cervical spine. As this study is mostly interested in the kinematics of cervical motions under a very small, pure moment (up to 2 N m), we decided to avoid preloading. This does not guarantee that all loading conditions (e.g., high-g acceleration) can be predicted without a preload. Intradiscal pressure in this study varied a lot between subjects. Niemeyer et al. [65] stated variation in the spine geometry strongly influences intradiscal pressure. In this study, an average of 0.22 MPa with a standard deviation of 0.19 MPa was obtained for extension under 1 N m load. This high variability agrees to the fact that vertebral geometry is influencing the disc pressure. Little and Adam [66] found a similar high variation in disc pressure in their study for lumbar spine.

Figure 12 shows the ROM for flexion/extension for all five subjects under a 2 N m load. Even though ROM varies between each subject, combined results with mean and standard deviations agree well with experimental results in Fig. 7. For example, subject 4 shows a lower ROM compared to others. Different geometries of vertebrae and ligaments’ insertion points influence the result. Similar trends can be noticed for lateral bending and axial rotation in Figs. 13 and 14. It is a matter of concern that the models are being compared to experimental data of different subjects (different ages, weight, sex, and/or stature). It is extremely difficult to match the subjects between simulation and experiment. The variation between subjects mostly originates from the differences in the geometry and dimension. Range of motion may get lower due to a smaller disc thickness. Individual subject specificity varies their geometry as well as material properties. Even though in this research all the material properties are derived from the literature, in reality, material properties are not the same in all subjects. There will always be uncertainty using finite element models of the cervical spine to apply to a real subject if both the model and the real subject are not similar. In the literature, researchers normally developed a generic cervical model to represent all subjects in spine biomechanics analysis. In Figs. 12–14, it is obvious that intersubject variations exist and it is critical to consider this variation in spine biomechanical analysis. Therefore, the major contribution in this study is to consider the intersubject variation by developing a series of cervical spine models.

One of the major reasons for the deviation from previous studies is due to the approximated boundaries of the intervertebral discs. Models in the literature used MRI scans to approximate the boundary of the intervertebral discs [25]. In this study, the discs are modeled by linearly connecting endplates without MRI scan data. Furthermore, this study adapted material properties from the literature [11,25] which were collected through controlled laboratory experiments for a specific subject. Some of material properties depend on the subject’s age, gender, level of hydration, and composition level [18]. A wide variation in the literature can be noticed in linear and nonlinear material properties. A preliminary material sensitivity test with this model can be available in the Supplemental Materials on the ASME Digital Collection. Sensitivity analysis on material properties shows a 25% change in the elastic modulus in flexion can lead to up to a 25% change in desired ROM. To save simulation time, all the materials in the vertebrae or discs are changed simultaneously in the material sensitivity study. Overall, vertebrae materials are less sensitive. Changes in disc material properties tend to show a larger effect on the output. A detailed sensitivity analysis of the individual material properties of the cervical spine can show the effect on the output and may be a useful tool to methodically determine the appropriate properties of desired materials [18]. Moreover, we believe changes in material properties will affect the output for lateral bending and axial rotation also. A comprehensive sensitivity analysis is ongoing.

In this study, five detailed finite element models of the cervical spine were developed as a first step to model digital twin for human cervical spine. The model included C2–C7 vertebrae, intervertebral discs, endplates, and nonlinear ligaments. A mesh convergence test was performed for verification of the model. The ROM was taken as a target property for convergence. After verification, the models in terms of mean and standard deviation were validated against existing well-established experimental results and simulation study results. Different loading conditions including flexion/extension, left and right lateral bending, and left and right axial rotation were applied for validation. The predicted intradiscal pressures were also compared with existing literature. With appropriate modification, it is expected these models can be applied to simulate real-life situations like whiplash in car accidents, disc degenerative disease, or high-g acceleration during flight maneuvers.

This research was partially supported by the National Science Foundation (NSF) (Grant No. CBET 1703093).

There are no conflicts of interest.

No data, models, or code were generated or used for this paper.