Authors: Samuel D. Schaefer, Benjamin M. Davies, Virginia F.J. Newcombe, Michael P.F. Sutcliffe
Categories: Short Communication, Biomechanics, Cervical spondylotic myelopathy, Computational case control study, Degenerative cervical myelopathy, Finite element analysis, Oscillation
Source: Brain & Spine
Degenerative Cervical Myelopathy [DCM] is a slow-motion spinal cord injury. Compression and dynamic compression have been considered disease hallmarks. However, this is likely an oversimplification, as compression is more commonly incidental and has only modest correlation to disease severity. MRI studies have recently suggested spinal cord oscillation could play a role.
To determine if spinal cord oscillation could contribute to spinal cord injury in degenerative cervical myelopathy.
A computational model of an oscillating spinal cord was developed from imaging of a healthy volunteer. Using finite element analysis, the observed implications of stress and strain, were measured in the context of a simulated disc herniation. The significance was bench marked by comparison to a more recognised dynamic injury mechanism; a flexion extension model of dynamic compression.
Spinal cord oscillation altered both compressive and shear strain on the spinal cord. Following initial compression, compressive strain moves from within the spinal cord to the spinal cord surface, whilst shear strain is magnified by 0.1–0.2, depending on the amplitude of oscillation. These orders of magnitude are equivalent to a dynamic compression model.
Spinal cord oscillation could significantly contribute to spinal cord damage across DCM. Its repeated occurrence with every heartbeat, draws parallels to the concept of fatigue damage, which could reconcile differing theories on the origins of DCM. This remains hypothetical at this stage, and further investigations are required.
Keywords: Finite element analysis, Degenerative cervical myelopathy, Cervical spondylotic myelopathy, Biomechanics, Oscillation, Computational case control study
Degenerative Cervical Myelopathy [DCM] is a form of spinal cord injury caused by narrowing of the spinal canal through degenerative and/or congenital changes (Davies et al., 2022a, 2022b). DCM is estimated to affect 1 in 50 adults and today is associated with significant disability, including loss of dexterity, imbalance and pain (Davies et al., 2022c). A recent estimate from England, United Kingdom suggested a cost to society of £0.7bn (Davies et al., 2022d). This is conservative, as it does not account for the widespread underdiagnosis (Grodzinski et al., 2022).
Despite this unmet need, the aetiology of DCM is poorly understood (Davies et al., 2022b; Badhiwala et al., 2020). At an elementary level we observe a chronic interaction of degenerative changes with the spinal cord, that leads to a progressive spinal cord injury and disability. This can be checked and often partially reversed through surgery to release the spinal cord. This has led to the notion that DCM is a disease of ‘spinal cord compression.’ However, compression has failed to provide a single unifying explanation for what is observed in these patients (Davies et al., 2022b; Witiw et al., 2017). For example, compression on MRI is most commonly incidental (Smith et al., 2020) and the amount of compression correlates poorly with disability or response to treatment (Martin et al., 2022). The spinal cord is a viscoelastic material and can tolerate some compression, and compression is not the only loading mechanism. It is more likely therefore that pathological loading arises from the interaction of a range of forces, modified by an individual's intrinsic vulnerability to sustain spinal cord injury and time (Davies et al., 2022b).
Longitudinal imaging studies of the degenerative spine, including DCM, highlight that there is seldom significant change in the relationship of degenerative changes with the spinal cord (Adamova et al., 2017; Bednarik et al., 2004, 2008; Nouri et al., 2022). Whilst injury could solely be explained by an overwhelming of the mechanisms that resist or retard injury (repair, or functional reserve capacity) (Davies et al., 2022b) it has also led to the proposal of a dynamic process – repetitive microtrauma (Badhiwala et al., 2020; Nouri et al., 2015). So far, this has largely been discussed in terms of the movement of the sub-axial spine, leading to ‘dynamic compression’ or ‘stretch’ loading (Nouri et al., 2015; Henderson et al., 2005). However recent work has refocused attention on spinal cord oscillation (Hupp et al., 2021).
Phase contrast imaging studies have demonstrated that the spinal cord oscillates with the cardiac cycle and this oscillation frequency and magnitude may change with cervical stenosis (Hupp et al., 2021; Wolf et al., 2018; Mikulis et al., 1994). Conceptually this motion pattern could lead to a further loading mechanism on the spinal cord affected by cervical stenosis (for example as the cord oscillates across the surface of a disc prolapse) and contribute to tissue injury (Davies et al., 2022b).
The objective of this study was therefore to explore whether oscillation could contribute to spinal cord loading, and therefore the aetiology of DCM, using finite element analysis [FEA].
This study used a simplified model to explore the estimated stress and strain on an oscillating spinal cord across a range of parameters. This is referred to as the oscillation model. Pre-clinical experiments from animal models have estimated strain thresholds of 0.1–0.2 for neural tissue injury (Russell et al., 2012). However to further benchmark the significance of the observed estimates in stress and strain to DCM (i.e., tissue injury), a separate model of ‘dynamic’ compression, more widely accepted to compound injury in DCM (Badhiwala et al., 2020), under similar conditions was first created. This is referred to as the flexion model. To ensure movement was consistent with physiological conditions, the flexion model included spinal vertebrae and ligaments, from C3 to T1. The spinal cord in each instance was modelled from the MRI of a healthy 22-year-old male. The study had obtained the necessary ethical approvals (Health Research Authority, IMAGE-DCM,20/EE/0037), and the participant had provided informed consent.
The spinal cord was extracted from the MRI using Spinal Cord Toolbox (Leener et al., 2016), which included segmentation of the grey and white matter. Experiments have shown these structures have different material properties (Davies, 2022). 3D Slicer (Kikinis et al., 2013) was used to convert these to separate mesh objects. Vertebrae meshes were obtained from the BodyParts3D database (Mitsuhashi et al., 2009). The mesh objects were then modified using Solidworks (Dassault Systèmes, 2018) and Abaqus (Dassault Systèmes, 2018) to obtain a 3D finite element half-model of the lower cervical spine (C3-T1).
Each vertebra was assumed to behave as a rigid body (Davies, 2022). The T1 vertebra was taken as a fixed reference. Discs were approximated as a uniform isotropic material. Ligaments were defined as axial connectors between adjacent vertebrae and given non-linear elastic properties according to Wheeldon et al. (2008) (Wheeldon et al., 2008).
To model DCM, a single level disc herniation was simulated. This is a common isolated finding in DCM (Nouri et al., 2016a). A portion of the disc was separated and pushed out into the spinal canal, deforming the healthy neural anatomy to create the diseased model. For the neural anatomy, the area of contact with the disc herniation experiences the most deformation, so element size was decreased in that area to model the behaviour accurately. Material properties were defined using a search of the literature and shortlisting process (Davies, 2022). Further detail of the flexion model is given in Supplementary Data A.
Simulations were performed using the Abaqus/Standard solver. General contact was enforced with a hard normal behaviour and a tangential friction behaviour with coefficient of friction 0.3. Flexion was achieved by applying a gradually increasing pure moment about the positive X axis to the top (C3) vertebra, up to 1.4 Nm, while fixing the bottom (T1) vertebra in place. A symmetric boundary condition was applied to all faces, edges and reference points on the midplane. These choices, alongside the mesh parameters, are summarised in Table 1.
The oscillation model used the same spinal cord, but for computational simplicity the disc prolapse was represented in isolation, as a static spherical cap. A compliant plate was located on the opposite side of the cord to allow compression of the cord.
The ends of the cord were held in place but allowed to rotate. The rear face of the compliant plate was fixed in place. General contact was enforced with hard normal behaviour and penalty loss tangential behaviour with a variable friction coefficient. The compliant plate had frictionless interaction with the cord. A symmetric boundary condition was applied to all faces, edges and reference points on the midplane. Further detail is given in Supplementary Data B.
The disc was moved via a reference point coupled to its outer face. To model disc prolapse, the reference point was pushed horizontally into the cord. To model spinal cord oscillations, the reference point was moved vertically according to decaying sinusoidal motion based on empirical measurements of the spinal cord (Mikulis et al., 1994).
Friction coefficient, spinal cord compression, and peak vertical displacement were variable parameters. Spinal cord compression and peak vertical displacement based on the aforementioned literature in DCM, are hypothesized to relate to DCM. Peak vertical displacement describes the oscillation amplitude. A fixed rate of decay was used in all cases. The benchmark criteria selected were a friction coefficient 0.15, compression ratio 3.6%, and peak vertical displacement 0.76 mm. The three parameters were studied independently by varying each parameter from this benchmark case.
The friction coefficient is an unknown for these surfaces, so was varied from the benchmark (0.15) down to 0 (i.e. frictionless behaviour), to determine whether overestimating the frictional stresses had significant impact on the results. A higher compression ratio increased the computation required to complete the simulation, so a maximum compression of 6.2% was used. This is relatively low compared to compression observed in severe cases of DCM (Tempest-Mitchell et al., 2019; Nouri et al., 2016b). Oscillation amplitude has been shown to increase with decreasing spinal cord diameter (Hupp et al., 2019). In DCM, this effect is seen around the location of canal stenosis, with a maximum observed peak displacement of 1.86 mm in DCM patients (Hupp et al., 2021; Tanaka et al., 1997).
Fields of strain and absolute stress from a section of the spinal cord near contact with the disc herniation were analysed in Python. For the node experiencing a maximum over all frames, the value at that node was plotted for all frames. Experimental work and previous FEA studies have found both principal strain and shear strain to be key indicators of tissue damage in brain tissue (Russell et al., 2012; Sutcliffe and Pan, 2021), so components of these were analysed. Elements far from the disc herniation were ignored as the coarser mesh size produced spurious local maxima (see Fig. 1).
Fig. 1 Flow chart illustrating construction of finite element models. Upper: Construction of finite element model for flexion study. Lower: Construction of finite element model for cord oscillation study.
Fig. 2 shows the minimum principal strain on the spinal cord for a section across the median plane (A) and the anterior face (B) for the start, intermediate, and end stages of flexion. Examining the direction of the principal axes showed this strain to be largely perpendicular to the cord surface, so gives compressive strain.
Fig. 2 Flexion results. A: Contour plots of minimum principal strain on the spinal cord across the median plane. B: Contour plots of minimum principal strain on the anterior face of the spinal cord. C: Peak compressive stress (left) and peak compressive strain (right) in the spinal cord, plotted against stage of flexion (given as the rotation of the C3 vertebra).
Three areas of elevated strain show the points of contact with the spinal two areas with a peak compressive strain around 0.1 at the vertebral arches, and a large area with a peak compressive strain of around 0.15 at the disc herniation.
The minimum principal stress and strain distributions approximately followed each other. The peak compressive stress and strain (negative minimum principal) are shown in Fig. 2C. Values are plotted against the stage of flexion, given as the rotation of the C3 vertebra.
Again, minimum principal strain is compressive strain on the cord. The disc herniation causes compressive strain on the spinal cord at the point of contact. The initial peak compressive strain is below the surface of the cord, at a value of 0.19. When oscillation starts, the peak compressive strain moves to the surface. The peak compressive strain, shown as the large red areas of Fig. 3A, moves with the point of contact between the disc and cord. The peak compressive strain increases slightly at the extremes of the oscillation.
Fig. 3 Oscillation results. A: Contour plots of minimum principal strain, which is largely compressive strain. Peak compressive strain labelled. Dashed lines indicate motion of the area of high strain. B: Contour plots of YZ shear strain before oscillation (left), at the peak of cord upstroke (centre), and at the peak of cord downstroke (right). C: Results for the three variables analysed in the oscillation study, showing shear strain oscillations at the location of peak shear strain. Horizontal bar signifies initial shear strain due to transverse compression; vertical arrows show range of shear strain experienced at that point.
For the benchmark case, contour maps of XY shear strain along the midplane are shown in Fig. 3B (where Y and Z are axes in the inferior-superior and anterior-posterior directions, respectively). Initially, the shear strain due to transverse compression has peaks above and below the disc, under the surface of the white matter at the interface with the grey matter. On the cord upstroke, the positive shear strain decreases, while the negative shear strain increases. On the cord downstroke, the reverse is seen, so that over a full cycle the shear strain oscillates about the initial value due to compression only.
These oscillations are summarised by the diagrams in Fig. 3C, showing the peak shear strain for each of the three studies. The horizontal bar signifies the initial shear strain due to transverse compression, while the vertical arrows show the range of shear strain experienced at that point.
Increasing the extent of transverse cord compression increased both the initial shear strain and the amplitude of oscillating shear strain. Increasing the peak vertical cord displacement caused slight increases in the maximum shear strain and large decreases during the cord upstroke. Increasing friction had negligible effect at low values, with shear oscillations still present when friction was removed entirely.
Using a computation model and FEA, we have illustrated that spinal cord oscillation can cause dynamic loading of the spinal cord. Notably the measured stresses and/or strains were comparable to a dynamic compression model under similar conditions (a mechanism more widely accepted to contribute to tissue injury in DCM) and thresholds for tissue injury established in preclinical models.
Whilst this only represents a hypothetical experiment and its generalisation should be cautious, these findings appear credible. Firstly, the model was built using anatomy extracted from human imaging, and the behaviour of each anatomical element carefully selected following a literature review (Davies, 2022). The model development therefore follows methods employed for existing FEA studies using DCM patients. Second, in both models, a compression strain at the level of the disc herniation was simulated in keeping with the theoretical basis of DCM (Badhiwala et al., 2020), it was simply magnified by dynamic loading. This included flexion and extension of the cervical spine in keeping with broader opinion, as well as oscillation. Specifically, during spinal cord oscillations, the location of peak loading moved along the cord, causing local areas to experience an oscillating compressive strain with an amplitude around 0.05. Significant shear strain oscillations were also experienced with an amplitude of 0.1–0.2, depending on characteristics of the oscillation. This is the preclinical threshold hypothesized to be sufficient for tissue injury (Russell et al., 2012). Finally, this would align with the observation that static and dynamic compression does not completely explain DCM (Martin et al., 2022), and the emerging imaging evidence suggesting a significance of cord oscillation in DCM, with changes observed (compared to asymptomatic individuals) and with some parameters demonstrating a relationship to clinical measures of disease severity (Hupp et al., 2021; Wolf et al., 2019).
This creates interesting implications for DCM theory. Whilst dynamic compression is clearly possible, conceptually how pervasive this is to all DCM has been more difficult to reconcile. For example, with cervical spondylosis either arthritis itself, or symptoms can reduce movement of the sub-axial spine (Binder, 2007) and the superiority of motion restricting surgery for DCM has not been clearly demonstrated (Yang et al., 2021; Ghogawala et al., 2021). Further, although spondylolisthesis has been associated with greater disease severity and poor surgical outcomes, it represents only a small proportion of DCM (Gondar et al., 2020; Nouri et al., 2020) with fewer still demonstrating segmental instability (Jiang et al., 2011). Conversely, spinal cord oscillation occurs with every heartbeat. It will therefore occur in all patients, with far greater frequency.
This would therefore support the concept that dynamic loading is a critical mechanism in DCM, but broaden the definition of ‘dynamic’ from just dynamic compression during flexion and extension. Further, the high frequency occurrence of spinal cord oscillation prompts interesting parallels to other concepts in engineering, in particular fatigue damage. Fatigue damage is the failure of a material after many cycles of loading, through the gradual propagation of cracks. It is more commonly considered in mechanical engineering with for example metals. However, the high-amplitude cycles of compressive and shear strains occurring with cervical stenosis could be relevant; assuming a heartrate of 80 bpm, one month of oscillations equates to 10^6^ cycles of shearing at significant amplitude. Such in-plane shearing is one mechanism of crack growth, suggesting crack propagation as a potential damage mechanism.
The vulnerability of the spinal cord to fatigue damage will depend on its material properties. Although little is known about the mechanical failure of spinal cord white matter, it may be analogous to unidirectional fibre composites due to the structure of aligned fibres (fascicles) embedded in an extracellular matrix (ECM). These unidirectional fibre composites are particularly vulnerable to damage in compression and in-plane shearing. The failure criteria in both cases are that of the matrix. Young, healthy ECM is a ductile material with high fracture toughness, so can absorb a lot of energy before crack propagation (Nyman and Makowski, 2012). However these properties change with age. In humans, this has largely been considered in the context of osteoporotic fractures (Martin, 1992) but also observed in neural tissue (Ryu et al., 2021). Aging is considered an important part of DCM, alongside other vulnerability factors.
This study demonstrates the potential for spinal cord oscillation to significantly contribute to spinal cord damage in DCM and extends the paradigm of ‘dynamic’ injury. It introduces a well-recognised engineering concept of fatigue damage to the field. Together, these may help to reconcile many differing mechanistic concepts. Whilst credible for the reasons outlined, these remain hypothetical at this stage, and further investigations are required.
BMD is supported by the National Institute for Health Research (NIHR) [https://www.nihr.ac.uk/] Clinical Doctoral Research Fellowship. The views expressed in this publication are those of the authors and not necessarily those of the NHS, the National Institute for Health Research or the Department of Health and Social Care. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
This research aligns with the AO Spine RECODE-DCM James Lind Alliance top research priorities, selected by people living and working with DCM. This includes “Biological Basis”, and to a lesser extent “Individualizing Surgery” and “Imaging and Electrophysiology”. For further information on how this process was conducted, why these questions were prioritised, and global updates on currently aligned research, please visit aospine.org/recode.
BMD is supported by the National Institute for Health Research (NIHR) [https://www.nihr.ac.uk/] Clinical Doctoral Research Fellowship. The views expressed in this publication are those of the authors and not necessarily those of the NHS, the National Institute for Health Research or the Department of Health and Social Care. This project was supported by an award from the Rostree Foundation with the Storygate Trust (A2844) The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Handling Editor: Prof F Kandziora
The following are the Supplementary data to this