Abstract
Spine degeneration is a normal aging process. It may lead to stenotic spines that may have implications for pain and quality of life. The diagnosis is based on clinical symptomatology and imaging. Magnetic resonance images often reveal the nature and degree of stenosis of the spine. Stenosis is concerning to clinicians and patients because of the decreased space in the spinal canal and potential for elevated risk of cord and/or osteoligamentous spinal column injuries. Numerous finite element models of the cervical spine have been developed to study the biomechanics of the osteoligamentous column such as range of motion and vertebral stress; however, spinal cord modeling is often ignored. The objective of this study was to determine the external column and internal cord and disc responses of stenotic spines using finite element modeling. A validated model of the subaxial spinal column was used. The osteoligamentous column was modified to include the spinal cord. Mild, moderate, and severe degrees of stenosis commonly identified in civilian populations were simulated at C5–C6. The column-cord model was subjected to postero-anterior acceleration at T1. The range of motion, disc pressure, and cord stress–strain were obtained at the index and superior and inferior adjacent levels of the stenosis. The external metric representing the segmental motion was insensitive while the intrinsic disc and cord variables were more sensitive, and the index level was more affected by stenosis. These findings may influence surgical planning and patient education in personalized medicine.
Introduction
From a functional perspective, the spinal cord is situated within a complex columnar structural arrangement of the spine. It lies posterior to the vertebral body and anterior to the lamina-spinous process [1]. The cord initiates from the mesencephalon in the lower brain region and continues inferiorly to occupy the three mobile regions of the vertebral column: cervical to thoracic to lower lumbar spine. While allowing motions between various intervertebral joints, the three regions of the spinal column maintain the integrity of the cord and preserve its functionality [2]. The osteoligamentous column protects the cord.
From a structural perspective, the human cervical spinal column transmits the in vivo load from the head to the torso and pelvis via a system of bones and joints [2,3]. At the superior end, C1 and C2 act as a transition structure from the base of the skull to inferior intervertebral joints from C2 to T1. Spinal ligaments from the base of the skull to T1 provide continuity to joints. While posteriorly located facet joints are present from the skull to T1, discs are present only from C2–C3 to C7–T1 levels in the cervical spine. The complex anatomical arrangement of multiple discs, facet joints, and ligaments at each level has a lordotic curvature to the column to balance the in vivo head mass, as its center of gravity lies posterior to the spine. Segmented joints coupled with lordotic curvature allow physiologic motions such as flexion and extension, respond to disease states by generally limiting the motions, and sustain other mechanical loads such as vibration and accelerative loading whiplash environments.
From a constitutive material perspective, the cancellous and cortical components of the vertebrae are stiffer than the other components of the spinal column [4]. The five major spinal ligaments in the subaxial spine and annulus fibers and nucleus pulpous of discs are considered softer materials with their own constitutive properties. The spinal cord is considered a soft tissue while it is heterogeneous with horn-like internal gray matter and an outside white matter. It is important to study the column and cord as a unit because of their coupling in form and function. Column response affects cord response and vice versa.
From a biomechanical perspective, the osteoligamentous column and cord sustain deformations under physiological and traumatic loading [2,3]. Although ligaments resist the loading uniaxially, vertebrae and discs of the column and cord resist the external load in a three-dimensional mode [5–8]. The responses of the column (disc and facet joints connected by ligaments and vertebrae) and cord are interrelated. To determine the response of the cervical spine to mechanical loading, clinical studies with normal subjects and patients, and experiments with postmortem human surrogates, animals, and physical models are conducted. Depending on the insult, human cadaver studies can be designed to apply loading and record kinematics and resisted loads at either end of the preparation. Local responses such as segmental range of motion are used to replicate clinical scenarios; however, to understand the intrinsic responses, computational modeling is necessary. This study uses finite element modeling as the tool to quantify responses under dynamic loading.
From a modeling perspective, finite element models of the cervical spine have been developed to investigate the biomechanics of the osteoligamentous spine [9–13]. Models have varied from isolated disc segments to functional units to multilevel motion segments to subaxial column and head-neck. The early purpose was to determine the biomechanical responses under physiological loading conditions such as flexion, extension, lateral bending, and axial rotation of the intact model, then the model was modified to simulate degenerated spines and different types of surgeries and develop medical devices for civilian environments. Results from human cadavers were used for validation. With the introduction of artificial discs as an adjunct to anterior cervical discectomy and fusion, finite element models have proliferated [14]. While these applications are not all inclusive, most models have focused on the responses of the osteoligamentous spine, with little regard to spinal cord simulations [15–17]. Although some isolated cord models exist, the interrelationship between the cord and column is not discussed. They are important for the understanding of the structural and functional biomechanics of the cervical spine. This study uses the column-cord as one structure to study its biomechanical responses via finite element modeling.
From an aging perspective, it is known that spine degeneration advances with increasing age [18]. Increasing levels of degeneration affect spine motions and may lead to stenosis. Increasing levels of stenosis is a major determinant of cord-related dysfunction. In the clinical literature, it is referred to as cervical myelopathy. Patients seek counsel on loads such as whiplash that can be endured without adverse effects. The objective of this study was to determine the response of the column and cord under postero-anterior acceleration loading. This was achieved by quantifying the external segmental range of motion, internal disc pressure, and cord stress and strain metrics.
Materials and Methods
Different grades of stenosis at the C5–C6 segment were simulated using an in-house developed and validated three-dimensional finite element model of the osteoligamentous subaxial cervical spinal column (Fig. 1) [19]. The model simulated the vertebral uncinate anatomy and consisted of the cortical shell and cancellous core of the bodies in the anterior column of the spine. From the dorsal to ventral region, the posterior elements of the vertebrae included the spinous process, bilateral laminae, bilateral pedicles, and bilateral mass, i.e., the facet pillars. The intervertebral disc included the endplates, annulus fibrosus, and nucleus pulposus which were sandwiched between the endplates and connected the superior and inferior vertebral bodies. The anterior and posterior longitudinal ligaments formed the outward anatomical layers of the vertebral bodies, while the capsular ligaments bounded the bilateral facet joints. The ligamentum flavum was attached to the laminae at each segment. The spinous ligaments connected the spinous process to lamina in the dorsal region of the column.
The cortical and cancellous bones of the spinal vertebrae were simulated as linear isotropic materials. The thickness of the cortical bone was 0.5 mm, and it encased the trabecular bone circumferentially. The thickness of the endplate was 0.2 mm, and it encased the disc along the axial direction. The ground substance of the disc was simulated using the hyper-elastic foam with Hill strain energy function, and the nucleus was simulated as a fluid element. Nonlinear ligament-specific stress–strain relationships were used to simulate individual spinal ligaments. The cord was simulated along with the head and musculature. The spinal cord components consisted of white and gray matter, pia mater, dura mater, denticulate ligaments, and cerebrospinal fluid. Material properties used in the spine and cord components were obtained from previously published literature studies [20–30]. Details of the model are given in Table 1. The model was subjected to postero-anterior acceleration loading at the first thoracic vertebra at a peak impact velocity of 6.5 km/h. Stenosis was simulated at the C5–C6 level, and they were classified as mild, moderate, and severe, with spinal canal diameters of 14 mm, 10 mm, and 6 mm, respectively. They were simulated by posteriorly protruding the mid disc, compressing the spinal canal, and altering the shape of the spinal cord.
Details of the finite element model
| Component | Element type | Constitutive model | Parameters |
|---|---|---|---|
| Cortical bone | Quadrilateral shell | Linear elastic | E = 16.8 GPa, ν = 0.3 |
| Trabecular bone | Hexahedral solid | Linear elastic | E = 0.442 GPa, ν = 0.3 |
| Endplates | Quadrilateral shell | Linear elastic | E = 5.6 GPa, ν = 0.3 |
| Facet cartilage | Quadrilateral shell | Linear elastic | E = 0.01 GPa, ν = 0.3 |
| Annulus ground substance | Hexahedral solid | Hill foam | C1 = 0.000115 GPa, C2 = 0.002101 GPa, C3 = −0.000893 GPa, b1 = 4, b2 = −1, b3 = −2 |
| Annulus fibrosus | Quadrilateral membrane | Orthotropic nonlinear | Fiber angle: 45–60 deg |
| Nucleus pulposus | Hexahedral solid | Fluid | K =1.720 GPa |
| Ligaments | Quadrilateral membrane | Nonlinear | Stress–strain curves |
| Spinal cord | Hexahedral solid | Viscoelastic | ρ = 1.04 × 103 kg/m3, C1 = 0.5345, C2 = 1.0665, C3 = 1.0113, G1 = 0.8927, G2 = 0.8926, G3 = 0.8917, β1 = −0.0137, β2 = 0.00775, β3 = 0.035 |
| Pia mater | Quadrilateral shell | Linear elastic | ρ = 1.13 × 103 kg/m3, E = 2.3 MPa, ν = 0.49 |
| Dura mater | Quadrilateral shell | Linear elastic | ρ = 1.130 × 103 kg/m3, E = 80 MPa, ν = 0.49 |
| Denticulate ligaments | Quadrilateral shell | Linear elastic | ρ = 1.040 × 103 kg/m3, E = 0.0058 GPa, ν = 0.45 |
| Spinal fluid | Hexahedral solid | Viscoelastic | ρ = 1.040 × 103 kg/m3, K = 2.19 GPa, G0 = 5 × 10−7, G1 = 1 × 10−7 |
| Component | Element type | Constitutive model | Parameters |
|---|---|---|---|
| Cortical bone | Quadrilateral shell | Linear elastic | E = 16.8 GPa, ν = 0.3 |
| Trabecular bone | Hexahedral solid | Linear elastic | E = 0.442 GPa, ν = 0.3 |
| Endplates | Quadrilateral shell | Linear elastic | E = 5.6 GPa, ν = 0.3 |
| Facet cartilage | Quadrilateral shell | Linear elastic | E = 0.01 GPa, ν = 0.3 |
| Annulus ground substance | Hexahedral solid | Hill foam | C1 = 0.000115 GPa, C2 = 0.002101 GPa, C3 = −0.000893 GPa, b1 = 4, b2 = −1, b3 = −2 |
| Annulus fibrosus | Quadrilateral membrane | Orthotropic nonlinear | Fiber angle: 45–60 deg |
| Nucleus pulposus | Hexahedral solid | Fluid | K =1.720 GPa |
| Ligaments | Quadrilateral membrane | Nonlinear | Stress–strain curves |
| Spinal cord | Hexahedral solid | Viscoelastic | ρ = 1.04 × 103 kg/m3, C1 = 0.5345, C2 = 1.0665, C3 = 1.0113, G1 = 0.8927, G2 = 0.8926, G3 = 0.8917, β1 = −0.0137, β2 = 0.00775, β3 = 0.035 |
| Pia mater | Quadrilateral shell | Linear elastic | ρ = 1.13 × 103 kg/m3, E = 2.3 MPa, ν = 0.49 |
| Dura mater | Quadrilateral shell | Linear elastic | ρ = 1.130 × 103 kg/m3, E = 80 MPa, ν = 0.49 |
| Denticulate ligaments | Quadrilateral shell | Linear elastic | ρ = 1.040 × 103 kg/m3, E = 0.0058 GPa, ν = 0.45 |
| Spinal fluid | Hexahedral solid | Viscoelastic | ρ = 1.040 × 103 kg/m3, K = 2.19 GPa, G0 = 5 × 10−7, G1 = 1 × 10−7 |
Metric refers to range of motion of the intervertebral joints, disc pressure, spinal cord stress, or spinal cord strain.
Finite element models of the severe (left column: top sagittal view and bottom axial view), moderate (middle column), and mild stenotic (right column) spines. Stenosis is at the C5–C6 level.
Finite element models of the severe (left column: top sagittal view and bottom axial view), moderate (middle column), and mild stenotic (right column) spines. Stenosis is at the C5–C6 level.
Results
At the superior level, for the mild stenotic spine, the range of motion, intradiscal pressure, spinal cord stress, and strain were 7.51 deg, 725 kPa, 17.8 kPa, and 0.4%, respectively. For the moderately stenotic spine, they were 7.57 deg, 708 kPa, 16 kPa, and 0.4%, respectively. For the severely stenotic spine, they were 7.64 deg, 760 kPa, 15.6 kPa, and 0.5%, respectively. At the index level, for the mild stenotic spine, the range of motion, intradiscal pressure, spinal cord stress and strain were 7.96 deg, 1025 kPa, 17.5 kPa, and 0.1%, respectively. For the moderately stenotic spine, they were 7.96 deg, 1152 kPa, 32.0 kPa, and 0.2%, respectively. For the severely stenotic spine, data were 8.08 deg, 1251 kPa, 42.3 kPa, and 1.8%, respectively. At the inferior level, for the mild stenotic spine, the range of motion, intradiscal pressure, spinal cord stress and strain were 7.99 deg, 139 kPa, 30.2 kPa, and 3.4%, respectively. For the moderately stenotic spine, they were 7.95 deg, 136 kPa, 31.8 kPa, and 4.5%, respectively. For the severely stenotic spine, they were 7.89 deg, 123 kPa, 30.8 kPa, and 4%, respectively. Figures 2–4 show the intradiscal pressure, spinal cord stress, and strain data for the three spines at the three levels.
With respect to the mild stenotic spine, the range of motion at the superior and index levels increased by up to 1%, and at the inferior level, it decreased by the same magnitude. With respect to the moderate stenotic spine, the range of motion showed the same trend: increased by up to 1.5% at the superior and index levels and decreased by 0.75% at the inferior level (Fig. 5). With respect to the mild stenotic spine, intradiscal pressures at the superior, and index levels increased up to 22%, whereas at the inferior level it decreased by up to 12%. With respect to the moderate stenotic spine, the intradiscal pressure showed the same trend: increased up to 9% at the superior and index levels and decreased by 10% at the inferior level (Fig. 6). With respect to the mild stenotic spine, spinal cord stress at the index level increased by up to 1.4 times, and at the superior level it decreased by up to 11%. With respect to the moderate stenotic spine, the spinal cord stress increased by 32% at the index level and decreased by up to 3% at the adjacent levels (Fig. 7).
Change in the range of motion of moderate and severe stenotic spines normalized with respect to mild spine, and severe stenotic spine normalized with respect to spine with moderate stenosis
Change in the range of motion of moderate and severe stenotic spines normalized with respect to mild spine, and severe stenotic spine normalized with respect to spine with moderate stenosis
Change in intradiscal pressure of moderate and severe stenotic spines normalized with respect to mild spine, and severe stenotic spine normalized with respect to spine with moderate stenosis
Change in intradiscal pressure of moderate and severe stenotic spines normalized with respect to mild spine, and severe stenotic spine normalized with respect to spine with moderate stenosis
Change in spinal cord stress of moderate and severe stenotic spines normalized with respect to mild spine, and severe stenotic spine normalized with respect to spine with moderate stenosis
Change in spinal cord stress of moderate and severe stenotic spines normalized with respect to mild spine, and severe stenotic spine normalized with respect to spine with moderate stenosis
Discussion
The objective of this study was to quantify the response of the column and cord for stenotic spines under postero-anterior acceleration using a finite element model. The original model was validated with human cadaver tests [12,31]. A mesh convergence study was done in addition to element quality checks such as the Jacobian and aspect ratio. Although the cord stress was not validated due to lack of published data, disc pressures and facet loads were validated [32,33]. These initial efforts provided the confidence to analyze stenotic spine cord and column responses. The range of motion at all the three levels from this study matched well with literature studies, and their magnitudes were below injury levels under physiologic loading [31,34–38]. A similar observation can be made of the intradiscal pressure [39]. Experimental studies for cord stress and strain during whiplash loading are sparse in the literature. While not the subject matter of this study, to determine the optimal metric that defines human tolerance to injury, i.e., cord stress, cord strain, or disc pressure, additional experimental and modeling studies are necessary.
The mild, moderate, and severe degrees of stenosis were simulated based on classifications reported in clinical literature: three, four, and five degrees of stenosis have been identified based on the appearance and extent of the disc material in the spinal canal and impingement to the cord [40–43]. Although a more granular definition covering the four- or five-degree classification system could have been used, as an initial step, a three-degree system was considered for determining the effects of varying stages of cervical disease on external and internal biomechanical variables. It should be noted that a clear distinction and consensus between different levels of stenotic classifications have not been reached based on numerical data on the impingement of the intervertebral disc material into the canal-cord space. The presently used 14 mm, 10 mm, and 6 mm of space were considered to represent the three degrees of stenotic canals, and they are generally accepted as normative values by clinicians, including the neurosurgeon author of this paper.
Based on clinical considerations, normalizations in this study compared changes in spines with moderate and severe stenosis with respect to the spine with mild stenosis and compared changes in severe stenotic spines with respect to the spine with moderate stenosis. Conservative or nonconservative treatment is often decided by the patient in consultations with the health care provider (surgeon), and decisions are based on the current status of disc impingement, patient's symptoms (neck pain), and expected future advancement of the disease and effect on the quality of life. Therefore, data were normalized to the “current” degree, i.e., when stenosis is in a mild stage, what could be the potential prognosis from future moderate and severe changes, and if the patient's anatomy and symptoms fall within the moderate category, from a biomechanical perspective, what could be the potential prognosis from future severe stenotic changes. While the study used a generic model to examine these issues, personalized finite element models of the spine are more applicable to specific patients. Because the present model was developed using a mapping block approach, it can be morphed to the specific anatomy of the patient to quantitatively determine these biomechanical variables and serve as decision-making tools for the patient and surgeon [19].
Motion at the index and superior and inferior levels represent the external response of the stenotic spine, and it is a clinical measure. The disc pressure represents the load-sharing by the anterior column of the spine [44]. Stress and strain metrics at three levels represent the intrinsic responses of the spine, focused on the cord. While motion can be directly measured using X-rays, disc pressure, and stress and strain cannot be easily measured. The rationale for the selection of the three levels was based on the consideration of whiplash loading effect on the immediate superior and inferior adjacent segments, and these levels are examined most frequently for treatment options. For example, decisions are often based on adjacent segment disease, and this is true for spinal fixations including artificial disc replacements [15–17]. Thus, the selection of the level and metrics have clinical and biomechanical implications.
From a clinical perspective, the C5–C6 level is the most frequently involved level in patients with spinal disorders such as cervical myelopathy. Superior and inferior levels were chosen because clinicians are interested in how the disease progresses over time, i.e., adjacent segment disease issues. The external metric of range of motion is commonly used to assess the mobility of the cervical spine, and flexion–extension X-rays are used for this purpose in a clinic. From a biomechanical perspective, the metrics of range of motion, disc pressure, and cord stress serve as external and internal responses of the spine to stenosis.
Because all four variables were obtained at all three levels, it is possible to assess their sensitivity with degree of stenosis. The relatively nonsignificant change in motion at any of the levels and for any degree of stenosis indicates that the presence of local stenosis does not appreciably alter segmental motion. Both cord stress and strain were found to be dependent on the degree of stenosis and spinal level. The greatest cord strains at the inferior level compared to its two superior levels indicate that the stretch migrates caudally for all degrees of stenosis. It should be noted, however, that the strain magnitudes were low. For both disc pressure and stress, in general, comparative/normalized changes were pronounced at the index level, C5–C6, while the actual magnitude of changes was lower for disc pressure than cord stress. Spine with severe stenosis sustained the greatest change in both pressure and stress measures. The inferior spinal level responded with decreased disc pressure while the superior level generally showed increased pressures (except moderate with respect to mild) suggesting that the inferior level is shielded by the increased disc pressure at the index level. These results suggest that an inferior disc may not advance to a degenerative state because of stenosis before the index level disc. From a cord stress perspective, the magnitudes of changes were considerably low at the adjacent levels compared to the increase at the index stenotic level. These findings suggest that cord stress and disc pressure are the two most sensitive parameters compared to the external range of motion, when defining the effects of whiplash on the stenotic spine.
This study used a generic model to determine the responses of stenotic cervical spinal cord and column. This approach delineated the role of stenosis without the confounding effects of other variables, commonly present in patients. Because the generic model outcomes are not patient-specific, conclusions drawn in this study are only a first step and are applicable broadly to stenosis cases. In most patients with degenerative diseases such as cervical myelopathy, stenosis will be associated with other geometrical or anatomical variations, and they can be quantified with medical imaging. Commonly obtained magnetic resonance images in patients can be used to incorporate the actual anatomical geometry. To better describe the relationship between a patient's stenotic anatomy and cord-column biomechanics, it would be necessary to develop patient-specific models. Findings from patient-specific models can then be used in personalized medicine for surgical decision making and patient counseling.
The range of motion is an integrated parameter as it represents the motion of the entire intervertebral joint with contributions from the disc and bilateral facet joints. Cord and disc mechanical parameters per se are not directly included in this external biomechanical metric. In contrast, cord stress and disc pressure are intrinsic metrics, focused on specific spine components. From a clinical standpoint, cord stress has important implications for risk of spinal cord injury and surgeries or interventions that reduce cord stress during whiplash will likely reduce the risk of spinal cord injury.
The original osteoligamentous finite element model was initially validated with human cadaver experiments in terms of range of motion at every level and in the time domain [12,31]. A mesh convergence study was also done in addition to element quality checks such as the Jacobian and aspect ratio. Although the cord stress was not validated due to lack of published data, disc pressures and facet loads were validated [32,33]. These initial efforts provided the confidence to analyze stenotic spine cord and column responses.
While this study simulated mild, moderate, and severe degrees of stenosis by altering the impingement of the disc material into the spinal canal, no other changes were made to the osteoligamentous spinal column. It is likely that severe stenotic changes also elicit changes to the column such as decreased disc height at the index level. The addition of such alterations to the severe stenotic spine may confound the results. Instead, focusing on a single change, as defined by the geometric impingement of the disc material, allowed the present modeling study to determine the pure effects of varying stenosis on column-cord mechanics. Because all three stenotic models were developed using a validated finite element model, the results are considered robust. Validations were done for segmental motions at all three levels, and disc pressures were also done. The spinal cord was not simulated in the originally validated finite element model. As indicated, it is difficult to validate against stress and strain of the cord, especially with stenotic spines, because of the intrinsic nature of the cord anatomy. Additional studies are needed to validate stenotic-specific finite element models.
The presence of asymptomatic spinal cord compression is seen in over 35% of patients over the age of 65 years. The risk of spinal cord injury after minor trauma in this patient population cannot be accurately quantified presently; however, cord compression predisposes to spinal cord injury. Using spine and spinal cord finite element modeling with simulation of whiplash injury, we can quantify the cord stress/strain and determine if these values are higher than known thresholds for spinal cord injury. Extending this approach to patient-specific finite element modeling, we can provide patients with a personalized risk assessment of spinal cord injury after whiplash loading. Additional parametric studies are necessary to determine the impact of spinal range of motion, alignment, and ankylosis on effects of whiplash on spinal cord biomechanics. Together, this approach will be a valuable addition to the clinical management of patients with cervical stenosis and spinal cord compression.
This study modeled stenosis by altering the impingement of the disc material into the spinal canal, while other changes associated with stenosis are likely present. The use of the single parameter variation allowed us to delineate the effects of stenosis without confounders that are generally present in patients. Addition of co-existing factors such as osteophytes that may change the disc-vertebral body relationship, osteopenia or osteoporosis that may soften the trabecular bone, and facet arthrosis that may compromise the joint anatomy, may play an additional, perhaps, cumulative role on the biomechanics of the stenotic spine. Inclusion of such associated alteration is best studied by patient-specific finite element models, and this study has laid a foundation to pursue those efforts.
Conclusions
Motions at the index and adjacent levels do not change significantly with stenosis, suggesting its insensitivity as a clinical measure. Intradiscal pressure is affected at the index level compared to adjacent levels suggesting that disc disease may initiate at this level before migrating to other levels. Although cord strains were low, they were greater at the inferior to the index and superior levels, suggesting that the inferior cord stretches more in stenotic spines. Cord stresses followed disc pressures, and these parameters appear to be most affected in stenotic spines. External responses are not always fully reflective of intrinsic responses in stenotic spines.
Acknowledgment
Drs. Yoganandan and Vedantam are employees of the VA Medical Center. It was supported by NASS and AO Spine. Opinions, interpretations, and conclusions are solely from the authors.
Funding Data
National Center for Advancing Translational Sciences and National Institutes of Health (Award No. UL1TR001436; Funder ID: 10.13039/100000002).
U.S. Army Medical Research and Materiel Command and Office of the Assistant Secretary of Defense for Health Affairs, through the Broad Agency Announcement (Award No. W81XWH-16-1-0010; Funder ID: 10.13039/100000182).
Department of Veterans Affairs Medical Center. It was supported by NASS and AO Spine.
Data Availability Statement
The authors attest that all data for this study are included in the paper.
