Maximum Principal Strain Correlates with Spinal Cord Tissue Damage in Contusion and Dislocation Injuries in the Rat Cervical Spine

Abstract

The heterogeneity of the primary mechanical mechanism of spinal cord injury (SCI) is not currently used to tailor treatment strategies because the effects of these distinct patterns of acute mechanical damage on long-term neuropathology have not been fully investigated. A computational model of SCI enables the dynamic analysis of mechanical forces and deformations within the spinal cord tissue that would otherwise not be visible from histological tissue sections. We created a dynamic, three-dimensional finite element (FE) model of the rat cervical spine and simulated contusion and dislocation SCI mechanisms. We investigated the relationship between maximum principal strain and tissue damage, and compared primary injury patterns between mechanisms. The model incorporated the spinal cord white and gray matter, the dura mater, cerebrospinal fluid, spinal ligaments, intervertebral discs, a rigid indenter and vertebrae, and failure criteria for ligaments and vertebral endplates. High-speed (∼ 1 m/sec) contusion and dislocation injuries were simulated between vertebral levels C3 and C6 to match previous animal experiments, and average peak maximum principal strains were calculated for several regions at the injury epicenter and at 1-mm intervals from +5 mm rostral to −5 mm caudal to the lesion. Average peak principal strains were compared to tissue damage measured previously in the same regions via axonal permeability to 10-kD fluorescein-dextran. Linear regression of tissue damage against peak maximum principal strain for pooled data within all white matter regions yielded similar and significant (p<0.0001) correlations for both contusion (R² =0.86) and dislocation (R² =0.52). The model enhances our understanding of the differences in injury patterns between SCI mechanisms, and provides further evidence for the link between principal strain and tissue damage.

Introduction

T raumatic spinal cord injury (SCI) often results in a debilitating condition, with estimated incidence rates of 1800 each year in Canada, and 12,000 per year in the U.S. (Farry and Baxter, 2010; National Spinal Cord Injury Statistical Center, 2009). A variety of treatments for traumatic SCI have been tested in recent decades, but none have proved widely effective for improving neurological outcomes in humans (Hyun and Kim, 2010). The heterogeneity of the SCI population is one possible reason for the lack of effective treatments, in that we do not fully understand the effects of important variables such as age, injury level, severity, and mechanism (Tator, 2006). Further research into these variables is necessary to guide substantial breakthroughs in targeted therapy development. One particular aspect that has received little attention until recently is the possibility of important differences in injury patterns created by the mechanism of primary injury, such as a spinal cord contusion from a vertebral burst fracture or a fracture-dislocation, the two most prevalent clinical mechanisms. Such differences could have implications for differential treatment of patient groups, and can most thoroughly be investigated with a combination of experimental and computational approaches.

In vitro cell culture experiments have shown that the biological and mechanical mechanisms of traumatic neuronal injury are influenced by mechanical loading patterns (Geddes-Klein et al., 2006; Cullen and LaPlaca, 2006). On a larger scale, dynamic injuries in the in vivo rat model have shown varying patterns of tissue damage for different injury mechanisms (Choo et al., 2007; Clarke et al., 2008; Fiford et al., 2004). Furthermore, finite element models of the human spine have demonstrated distinct stress and strain patterns within the spinal cord that depend on the biomechanical mechanism of injury (Greaves et al., 2008; Li and Dai, 2009).

While experimental methods have revealed important details regarding tissue injury thresholds and patterns, they are difficult to fully interpret and apply to clinical injuries, as they do not yield information on internal spinal cord deformations during SCI. Finite element models are ideal for investigating these deformations, but human models to date have been difficult to validate due to a lack of in vivo loading data, and have used linear elastic quasi-static simulations which may not capture the full nature of high-speed cord injuries. One group addressed this with an experimentally-calibrated dynamic finite element model of weight-drop contusion in the rat thoracic spine, but did not investigate other injury mechanisms (Maikos et al., 2008b).

The objective of this study was to develop and validate a dynamic finite element model of the rat cervical spine, and use it to compare internal spinal cord deformations during contusion and dislocation injury mechanisms with previously observed tissue damage. The indicator of neuronal tissue damage used for this comparison was cellular permeability to fluorescein-dextran, as membrane permeability has widely been linked to neuronal pathology (Galbraith et al., 1993; Shi and Blight, 1996; Stys, 2004; Sullivan et al., 2005), and has been used to quantify regional patterns of damage in the spinal cord (Choo et al., 2007). We chose maximum principal strain as the primary measure of cord deformation because it has been shown to be a good predictor of neural damage for several animal models (Bain and Meaney, 2000; Maikos et al., 2008b; Shreiber et al., 1997; Zhu et al., 2006). The maximum principal strain in soft tissues can be interpreted as the greatest tensile distortion of the tissue element. Note that two regions of tissue may experience the same maximum principal strain, but the orientation of the distortion may differ. Improved understanding of the strain distribution in the cord during two distinct mechanisms of SCI will aid interpretation of tissue damage patterns, and may inspire new strategies to treat or prevent injury.

Methods

Model creation

The first step in development of the rat cervical spine model was acquisition of spinal geometry to be used in creating a finite element (FE) mesh. High-resolution magnetic resonance imaging (MRI) was performed with a 7-T animal scanner at the University of British Columbia MRI Research Centre. Two scans, oriented perpendicular to the upper cervical cord (C1–C3), and lower cervical cord (C5–T1), were obtained with 156×156-micron in-plane and 1-mm through-plane resolution. These scans were interpolated to an isotropic 156-micron resolution, zero-padded, and registered to each other using Analyze (AnalyzeDirect, Overland Park, KS), yielding a fused image of pixel dimensions 256×256×240 covering the full range of the cervical spine (Fig. 1).

FIG. 1.

Geometry of the rat cervical spine was extracted from 7-T magnetic resonance images (MRI) using ITK-SNAP, a semi-automated volume segmentation tool using 3-D snake evolution (Yushkevich et al., 2006). A screenshot of the ITK-SNAP interface is shown after segmentation of the white and gray spinal cord, intervertebral discs, and C1–T2 vertebrae. Clockwise from the top left: Axial, sagittal, and frontal views of the MRI data are displayed, with a 3-D cursor that links all three. A 1-cm-long scale bar is shown alongside the spine volumes at bottom left. Color image is available online at www.liebertonline.com/neu

Image segmentation was performed using the semi-automated volume segmentation software ITK-SNAP, which is based on 3-D snake evolution (Yushkevich et al., 2006). Geometry for the white and gray matter of the spinal cord, intervertebral discs, and C1–T2 vertebrae, were exported from ITK-SNAP into Rapidform (INUS Technology, Seoul, Korea), and fit with analytical surfaces using non-uniform rational B-splines (Piegl and Tiller, 1997). These surfaces were meshed initially in HyperMesh (Altair Engineering, Troy, MI), using hexahedral solid elements for the white and gray cord, and tetrahedral elements for the discs and vertebrae. The meshes were imported into PAM-CRASH (ESI Group, Paris, France), an explicit FE software program suitable for impact simulations, for further development. The dura mater could not be reliably identified in the acquired MRI scans, and was instead created by expanding the surface of the cord based on MR images outlining the cerebrospinal fluid (CSF) by Franconi and associates (2000), then assigned a thickness of 90 μm, and meshed with two layers of hexahedral solid elements. Finally, spinal ligaments were created by manually defining two-dimensional bar elements according to anatomical descriptions similarly to the method used by Greaves and colleagues (2008) in a human cervical spine model.

A four-vertebra section of the cervical model corresponding to the local region of the experiments by Choo and associates (2007; C3–C6) was isolated for use in the current study to reduce computational complexity (Fig. 2). Maikos and colleagues (2008b) found that changing model length had minimal impact on simulation results, due to the localized nature of the weight-drop injury, therefore this truncation is not expected to influence simulations of contusion and dislocation injuries for which the axial extent of injury is concentrated within±3 mm from the injury epicenter (Choo et al., 2007), a region spanning roughly half of the C3–C6 model section length. Preliminary simulations demonstrated that increasing element size by up to a factor of 8 times did not change strain patterns, resulting in a fine mesh being chosen with ∼ 96,000 white matter and ∼ 74,000 gray matter elements. Simulations were performed on a quad-core 3.0 GHz desktop computer.

FIG. 2.

(A) Superior view of the rat cervical spine model, with white and gray matter, dura mater, and cerebrospinal fluid (CSF) constituents identified. (B) Lateral view of the four-vertebra cervical model, outlining intervertebral (IV) discs and spinal ligaments. Color image is available online at www.liebertonline.com/neu

Spinal cord and dura

Accurate spinal cord material properties are crucial to yield reliable tissue deformation during FE simulation. While the difficulty of measuring and modeling soft tissue properties has so far precluded gathering enough data for a consensus on cord material properties, material tests of human and animal spinal cord have all demonstrated hyperelastic and viscoelastic behavior (Bilston and Thibault, 1996; Cheng et al., 2008; Fiford and Bilston, 2005; Sparrey and Keaveny, 2011). We chose to use the same hyperviscoelastic Ogden and Prony material properties for the cord and dura presented by Maikos and associates (2008b), for consistency and to further test these material models in different injury conditions (Table 1). Cord hyperelastic properties were based on material tests of rat spinal cord by Fiford and Bilston (2005), combined with viscoelastic properties of brain tissue from Mendis and colleagues (1995), and then calibrated by Maikos and co-workers (2008b) to fit their own weight-drop experimental behavior at impact velocities of 0.489–0.690 m/sec. Properties for the dura were derived by mechanical testing of rat dura mater (Maikos et al., 2008a). Appropriate conversions of material constants were performed because of differing notation conventions between PAM-CRASH (Material Type 38) and Abaqus (Table 1), and elements from the dura and cord used selective reduced integration for the optimal balance between accuracy and computational complexity.

Table 1.

Material Properties of Spinal Cord and Dura

Tissue	Hyperelastic Ogden constants	Viscoelastic Prony series constants
Spinal cord	μ=40.04 kPa^a	g₁=0.5282^b
	α=4.7^c	τ₁=8 msec^c
	ν=0.45^c	g₂=0.3018 ^b
		τ₂=150 msec^c
Dura	μ=207.41 kPa^a	g₁=0.3182 ^b
	α=16.2^c	τ₁=9 msec^c
	ν=0.45^c	g₂=0.1238 ^b
		τ₂=81 msec^c
		g₃=0.0997 ^b
		τ₃=564 msec^c
		g₄=0.0997 ^b
		τ₄=4.69 sec^c

Adapted from Maikos et al. (2008b) with conversion μ_PC=G_0,AB/α due to differences in notation convention between PAM-CRASH (PC) and Abaqus (AB).

Adapted from Maikos et al. (2008b) with conversion g_i,PC=G_i,AB/G_0,AB.

Taken from Maikos et al. (2008b).

Cerebrospinal fluid

Persson and associates (2011) recently demonstrated the importance of including the incompressible fluid behavior of the CSF in models of SCI, using an ovine FE model with fluid-structure interaction. Previously, FE models of SCI have omitted the CSF (Greaves et al., 2008; Li and Dai, 2009; Scifert et al., 2002), or modeled it as a quasi-fluid using solid elements (Maikos et al., 2008b). This study proposes the smoothed particle hydrodynamics (SPH) method (Monaghan, 1988) as an efficient means to include interactions between the cord, dura, and CSF in impact simulations.

HyperMesh (Altair Engineering) was used to define simple cubic elements distributed regularly in the volume between the dura and cord elements (mesh pitch of 0.075 mm and ∼153,000 elements). These were converted to SPH point elements in PAM-CRASH, and a Murnaghan Equation of State model (Material Type 28) was used to model the fluid, with pressure defined relative to current and initial density (set to 0.001 g/mm³) as: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} P = B \left( \left( \frac { \rho } { \rho_0 } \right) ^7 - 1 \right). \end{align*} \end{document}

This model was proposed previously (Monaghan, 1994) as an efficient means for modeling fluid flow when the fluid velocity is much lower than its speed of sound propagation. A parameter, B, is set to artificially reduce the speed of sound in the fluid in order to decrease the minimum solution time step and thus increase computational efficiency. This strategy is shown to have minimal effects on the density variations and fluid behavior, provided that the reduced speed of sound is maintained at least 10 times the maximum flow velocity (Monaghan, 1994). To be conservative, preliminary simulations were run by reducing the value of B until the fluid simulation was no longer the limiting factor in the minimum simulation time step, yielding an optimal value for B in the current simulations of 200 MPa. This value set the speed of sound in the CSF at roughly 400 times the maximum flow velocity of ∼3 m/sec.

Spinal ligaments

Although spinal ligaments do not play a role in experimental cord contusion injuries, they do play an important role in dislocation injuries, taking on part of the load applied to the vertebrae alongside the intervertebral disc. A nonlinear tension only bar (Material Type 205) was used for all spinal ligaments, with linear elastic properties (Greaves et al., 2008). Cross-sectional areas for the ligaments were scaled from the human values used in the model by Greaves and colleagues (2008), and maximum strains were assigned based on values reported in the literature (Lee et al., 2006; Quinn and Winkelstein, 2007; Yoganandan et al., 2000). All material property data for ligaments are shown in Table 2.

Table 2.

Material Properties of Spinal Ligaments (Material Type 205–Nonlinear Tension Only Bar Element)

Ligament	Cross-sectional area per element A (mm²)	Young's modulus E (MPa)	Stiffness factor k'=E ^ A (N)*	Mass per unit length μ (g/mm)	Failure strain EPSLN_u^e
Anterior longitudinal	0.344^a	35.2^b	12.11	0.000344	0.308^d
Posterior longitudinal	0.396^a	35.7^b	14.14	0.000396	0.182^d
Joint capsule	3.956^a	4.9^b	19.38	0.003956	1.51^c
Ligamentum flavum	0.674^a	3.8^b	2.56	0.000674	0.77^d
Interspinous	1.118^a	5.0^b	5.59	0.001118	0.609^d
Denticulate	0.0215^a	5.8^b	0.1247	0.0000215	0.087
Dural attachments	0.1075^a	35.7^b	3.838	0.0001075	0.182

Scaled from Greaves et al. (2008).

Taken from Greaves et al. (2008).

Taken from Quinn and Winkelstein (2006).

Taken from Yoganandan et al. (2000).

EPSLN_u is the PAM-CRASH parameter name for ultimate Failure Strain, ɛ_u.

Intervertebral discs

The annulus fibrosus of the C4/C5 intervertebral disc was modeled as linear elastic (Material Type 1; E=2.4 MPa), according to the properties used previously by Greaves and associates (2008). The nucleus pulposus was not modeled for this study to reduce complexity, as it is not expected to play a role in the contusion or dislocation mechanisms. Discs were attached to neighboring vertebrae via spot welds to simulate the vertebral endplate connection. Preliminary simulations were performed to calibrate spot-weld rupture criteria to achieve simulated disc endplate failure coinciding with experimental failure as predicted by Choo and colleagues (2007), based on force history measured during injury; the resulting criteria were ultimate tensile and ultimate shear strengths of 0.15 MPa.

Loading and boundary conditions

The loading and boundary conditions for the contusion and dislocation simulations were modeled to recreate the experiments by Choo and associates (2007, 2008) as closely as possible (Fig. 3). All vertebrae (C3–C6) were modeled as rigid, as well as the intervertebral discs, with the exception of the disc directly at injury epicenter (C4/C5). Friction was not included in the model, as Maikos and colleagues (2008b) found this to have negligible influence. The ends of the dura were constrained to prevent axial motion in order to encourage biofidelic membrane behavior and avoid flapping (Maikos, 2008b). No boundary conditions were imposed on the CSF particles at either open end of the model, as preliminary simulations showed minimal fluid leakage in the short time period until peak displacement, and to avoid non-biofidelic reflections at these locations from confusing results. Simulation validation was attempted by comparing applied forces to the experimentally reported values (Choo et al., 2007,2008). The spinal cord contusion force was validated with a more recent study (Choo et al., 2008), in which improvements in instrumentation provided a more accurate measurement of the smaller forces measured during contusion.

FIG. 3.

Sagittal cross-sections of the entire spine model (C3–C6) are shown, demonstrating contusion and dislocation injury mechanism simulations during displacement. Color image is available online at www.liebertonline.com/neu

For the contusion simulation, elements from the C4 and C5 vertebrae were removed to represent the partial laminectomy performed experimentally. This opening made way for a rigid indenter, modeled after the 2-mm spherical-headed steel indenter and aligned normal to and directly at the dural surface. Indenter motion was enforced by applying a velocity ramp from 110 cm/sec (corresponding to the average experimental peak velocity just prior to impact), down to 0 over 3.2725 msec, and continuing to −110 cm/sec to return to the starting position. Note that in the FE simulation, strictly defining the velocity profile over time also results in the corresponding displacement profile being strictly enforced, unlike in experimentally-controlled devices, in which some form of feedback is guiding the control. This trajectory results in a peak indenter displacement of 1.8 mm, and a peak cord compression of 1.08 mm.

Prior to the dislocation simulation, the facet joints and dorsal ligaments between the C4/C5 vertebrae were removed, as was performed experimentally to increase injury repeatability by eliminating residual facet dislocation. The C3 and C4 vertebrae were constrained in all directions, while C5 and C6 were displaced 2.5 mm dorsally from rest by applying velocity ramps up to 95.1 cm/sec over 2.629 msec, down to −95.1 cm/sec, and back to rest.

Regional correlation with tissue damage

Regional zones were outlined in the model and spaced at 1-mm axial slice intervals from the injury epicenter to reflect the dorsal, lateral, ventromedial, and ventrolateral white matter and ventral gray matter regions in which Choo and associates (2007) quantified membrane permeability to dextran. For each injury mechanism, mean membrane permeability values (axons/mm, or percent cells in the gray matter) were plotted against mean peak values of maximum principal strain (mm/mm) for each region by matching data according to slice position.

The data from the white matter regional correlations were pooled to investigate whether an overall relationship existed between strain and tissue damage. Data from the gray matter region were analyzed separately, as these data were quantified differently (percent cells) from white matter regions (Choo et al., 2007). Data corresponding to regions at the injury epicenter for the dislocation mechanism were excluded from all correlations because immediate hemorrhagic necrosis in this region precluded accurate quantification of dextran-positive cells and axons.

Statistical analysis

Linear regression was used to quantify correlations between maximum principal strain in each region, and for the pooled white matter data and the gray matter region, with R² correlation coefficients calculated for each mechanism, along with corresponding p values.

Results

We validated the FE model by comparing predicted loads from the enforced displacement profiles to the corresponding experimental corridors (Fig. 4). The simulated contusion force, applied directly to the dura and cord, followed a very similar pattern to the mean experimental force, but was delayed relative to the experimental traces (Fig. 4A). The peak force for the contusion simulation was 1.4 N, which was close to the experimental mean of 1.5±0.4 N (±SD). Forces applied during dislocation were much higher than direct cord contusion forces, as they were instead applied to the entire vertebral column. The simulated dislocation force demonstrated a time lag that was not observed as strongly experimentally (Fig. 4B), but the force-time curve followed a similar path to the experimental corridor. The peak dislocation force of 17.6 N was within the experimental range, but below the mean of 24.7±5.7 N. In addition, both experimental and simulated dislocation force traces demonstrated multiple local peaks that indicate sequential failure of local soft tissue components such as the intervertebral disc and spinal ligaments.

FIG. 4.

Simulated applied forces for contusion and dislocation injuries are plotted alongside corresponding experimental traces and corridors. (A) The simulated contusion force, applied directly to the dura and cord, follows a very similar pattern to the mean experimental force, but is seen to lag behind. (B) Applied dislocation forces are much higher than direct cord contusion forces, as they are instead applied to the entire vertebral column. While the simulated dislocation force deviates from the experimental corridor, the peak force is within the experimental range, and both experimental and simulated traces demonstrate multiple local peaks that indicate sequential failure of the intervertebral disc and spinal ligaments. Color image is available online at www.liebertonline.com/neu

Mid-sagittal images of the deformed spinal cords with internal strain distributions highlight differences between the mechanisms (Fig. 5). Some of the most striking features were the high strain (>0.16) “tails” seen dorsocaudal and rostroventral to the dislocation injury epicenter (Fig. 5B), showing the local regions subjected to tension due to rotation of the cord at the epicenter. These rostrocaudal asymmetries about the epicenter were evident in the affected dorsal, ventromedial, and ventrolateral white matter region strain plots (Fig. 7F, H, and I). Tissue damage in both ventral white matter regions also showed increases rostral compared to caudal to the epicenter (Fig. 7C and D), though the caudal increase in the dorsal white matter was not evident in the experimental results (Fig. 7A).

FIG. 5.

Distribution of maximum principal strain during contusion (A) and dislocation (B) injury simulations. Sagittal slices are shown at the midline and±1 mm. Contusion plots show a relatively small very-high-strain zone (>0.4 maximum principal strain) localized to the dorsal surface and extending less than 1 mm lateral to the midline, while the very-high-strain zone for the dislocation zone encompasses the full depth of the cord and extends beyond 1 mm laterally, demonstrating differences in mechanical injury severity and extent between the two mechanisms. Other notable features are the high strain (>0.16) “tails” observed dorsocaudal and rostroventral to the dislocation injury epicenter, showing the local regions subjected to tension due to rotation of the cord at the epicenter. Color image is available online at www.liebertonline.com/neu

Strong overall correlations between maximum principal strain and tissue damage indicated the model's biofidelity and corresponding utility. For the pooled white matter regions, maximum principal strain showed significant (p<0.0001) correlations with tissue damage for both contusion (R² =0.86; Fig. 6A) and dislocation (R² =0.52; Fig. 6B) mechanisms, with damage tending to increase with increasing strain (Table 3). For the ventral gray matter, maximum principal strain correlated strongly with tissue damage for both contusion and dislocation mechanisms (both R² =0.93; Fig. 6C and D). Maximum principal strain distributions (Fig. 7) were similar in nature to the distributions of membrane permeability for contusion and dislocation mechanisms, with central peaks flanking the injury epicenter, and decreasing tails toward the caudal and rostral extremes. All non-pooled regions yielded high correlation coefficients for contusion (R² between 0.90 and 0.96), as did three of the five regions for dislocation (R² between 0.61 and 0.93), with notable exceptions being the dorsal and lateral white matter (R² =0.38).

FIG. 6.

Correlations between maximum principal strain and axonal membrane permeability in contusion (left, blue) and dislocation (right, red), using pooled data from all four white matter regions (A and B), and data from the ventral gray matter (C and D). Linear regression lines of best fit are plotted as solid lines, along with 95% confidence intervals as dashed lines. Pooled white matter regions show significant (p<0.0001) and similar correlations between maximum principal strain and tissue damage for contusion (R² =0.86) and dislocation (R² =0.52) mechanisms. Correlations within the gray matter are also strong (R² =0.93), and the slope of the tissue damage versus strain regression is steeper compared to the pooled white matter result, especially for dislocation. Color image is available online at www.liebertonline.com/neu

FIG. 7.

Rostrocaudal distributions of experimentally-measured tissue damage and computed maximum principal strain are plotted at 1-mm intervals for both contusion (blue, solid) and dislocation (red, broken). (A–E) Tissue damage as measured by counts of cells permeable to dextran is generally higher for dislocation than contusion, with peaks near the injury epicenter, trailing down to lower levels rostral and caudal. Local dips directly at the injury epicenter for dislocation injuries are the result of immediate hemorrhagic necrosis (these points were excluded from statistical analyses and correlations). (F–J) The average peak maximum principal strains for each cord region show highest strains near the injury epicenter, as expected, and higher peak strains for dislocation than contusion. The higher average peak strains at the epicenter for dislocation than contusion in the lateral white matter (G) and ventral gray matter (J) regions reflect the larger lateral extent of the very-high-strain zone for dislocation as seen in Figure 5, and are supported by especially low epicenter cell counts due to rampant necrosis in these areas (B and E). In addition, rostrocaudal asymmetries about the dislocation epicenter corresponding to the high-strain tails from Figure 5 are seen in the affected dorsal (F), ventromedial (H), and ventrolateral (I) white matter strain plots, matching the patterns of tissue damage seen in the ventral white matter (C and D), though the caudal increase in the dorsal white matter (F) is not evident in the experimental results (A). Error bars denote standard error for experimental cell permeability counts (A–E), and standard deviation of peak strains for all elements in the simulation regions (F–J). Color image is available online at www.liebertonline.com/neu

Table 3.

Correlation Coefficients for Maximum Principal Strain and Tissue Damage within Cord Regions

	Contusion		Dislocation
Cord region	Correlation coefficient, R²	p Value	Correlation coefficient, R²	p Value
Ventral gray	0.93	1.5e–6	0.93	6.3e–6
Dorsal white	0.94	5.9e–7	0.38	0.056
Lateral white	0.96	1.1e–7	0.38	0.056
Ventromedial white	0.90	8.0e–6	0.82	3.1e–4
Ventrolateral white	0.91	5.6e–6	0.61	0.0075
Pooled white	0.86	2.1e–19	0.52	1.8e–7

Key differences between the contusion and dislocation mechanisms lie in the size and shape of the central zones of very high maximum principal strain (>0.4), with the contusion zone located only near the dorsal surface, and extending less than 1 mm lateral to the midline, while the dislocation zone encompassed the full depth of the cord and extended beyond 1 mm laterally. This corresponded to much higher average peak strains at the epicenter for dislocation than contusion in the lateral white matter (Fig. 7G) and ventral gray matter (Fig. 7J), and was reflected in the experimental dislocation results (Fig. 7B and E) by markedly low epicenter counts due to primary axotomy and widespread necrosis in these areas, respectively.

Discussion

We developed this FE model of contusion and dislocation injuries in the rat cervical spine to address the lack of knowledge regarding the effect of injury mechanism on primary mechanical damage patterns in the spinal cord. It is the first multi-mechanism computational model to be based on experimental injury mechanisms in an animal model, and thus has the advantages of more direct validation and comparison with histological tissue damage. In addition, this study demonstrated the feasibility of using the smoothed particle hydrodynamics method to model the CSF during impact, which may be useful in large animal or human SCI FE models (that involve larger subdural spaces and CSF volumes compared to the rat), in which the CSF has been shown to play an important role to cushion impact to the cord (Jones et al., 2008; Persson et al., 2009).

Overall, the model demonstrated its versatility to simulate both contusion and dislocation injury mechanisms with good biofidelity. The hyperviscoelastic material properties of the spinal cord yielded a realistic contact force during contusion, and the material properties of the intervertebral disc and spinal ligaments (including failure strain limits) resulted in dislocation force profiles that were similar to those measured in animal models. The simulated force history for dislocation, however, deviated somewhat from the experimental corridors, including a prolonged toe region of increasing stiffness at the start of displacement, which indicates that the involved disc and ligaments were not behaving stiffly enough in this initial phase, possibly due to inaccurate pre-tension in ligaments, omission of muscle attachments, or inability of the linear elastic material properties of the disc and ligaments to model behavior accurately in this regime. Because of this, the current ligamentous cervical spine model cannot be considered fully validated and should not be used without further refinement to model more general, external perturbations to the spine, such as rear or head-first impacts. However, this limitation on the biofidelity of gross spinal column forces during dislocation did not affect the time course or amount of spinal cord deformation in our study, as this was determined by contact with displacement-controlled vertebrae. Indeed, maximum principal strain was shown to correlate well to tissue damage for both contusion and dislocation cervical injury mechanisms, extending the previous findings for thoracic weight-drop contusion reported by Maikos and colleagues (2008b).

The correlations between maximum principal strain and tissue damage presented here did not suggest any obvious damage thresholds, such as a minimum strain level below which no damage was observed, that could be used as an injury criterion. However, all regions subjected to at least 0.1 maximum principal strain corresponded to elevated average levels of tissue damage, while averages for regions less than 0.1 strain varied between baseline and moderate levels of damage. The variation below 0.1 strain was especially high within the dislocation results, possibly due to less repeatability for this mechanism. This loose threshold of 10% strain is slightly below the lower bounds of tissue damage thresholds of 13–19% found previously (Bain and Meaney, 2000; Shreiber et al., 1997). Furthermore, steeper slopes for the ventral gray matter correlations compared to the white matter are suggestive of a lower injury tolerance to tensile strain in the gray matter, though it should be noted that the correlations for gray matter are less conclusive due to the small number of samples in that region in the current study; future studies with a focus on better quantifying such differences in injury tolerance are certainly warranted.

Interestingly, although the computational results for contusion injuries show peak strains beneath the tip, histology following contusion shows damage focused in the gray matter. This discrepancy is due to the lower injury threshold of the highly vascularized gray matter in comparison to white matter, as was found previously by Maikos and Shreiber (2007) in weight-drop contusion injuries. At high enough strains, one would expect to see primary damage in the white matter, and indeed this is the case. The models predict greater strain in the lateral white matter during dislocation injury (>0.4 in Fig. 7G), compared to the dorsal white matter during contusion (∼ 0.3 in Fig. 7F). Accordingly, the dorsal white matter is often spared during contusion, whereas primary damage to the lateral white matter is common in our animal models of fracture-dislocation (see loss of axons at the epicenter in Fig. 7B).

In addition to correlating strain and tissue damage and investigating injury tolerances, the model is also a useful tool to compare injury severity between mechanisms. For the mechanisms we studied, which were developed previously by Choo and colleagues (2007), dislocation appears to be much more severe in both peak maximum principal strain intensity and extent. This bears some similarity to the results of Greaves and associates (2008), whose quasi-static human model showed a deeper and wider extent of von Mises strain for dislocation than for contusion. Such comparisons and predictions of injury mechanism severity will be useful for further development of consistent and well characterized injury protocols, and, alongside behavioral survival studies, can increase our understanding of the differences seen in functional deficits and treatment goals between mechanisms.

Several limitations of the current work suggest possible improvements for the future. The hyperviscoelastic cord properties proposed by Maikos and associates (2008b), and further validated by our study, model spinal cord behavior quite well, but were based initially on material testing data, and then adjusted to better match experimental behavior. Further material testing and modeling of the rat cord is necessary, including investigation of white and gray matter inhomogeneity and anisotropy of the tissue, to determine the importance of such factors in modeling.

Furthermore, there is no current consensus on possible differences in white and gray matter material properties (Arbogast and Margulies, 1998; Coats and Margulies, 2006; Ichihara et al., 2001; Maikos et al., 2008b). We omitted such differences due to the ongoing uncertainty in their exact numerical values. Relative differences in the material properties of white and gray matter would cause the two components to deform by different magnitudes, with greater deformation in the softer material than the stiffer material; this difference in deformation would result in an additional shearing strain at the interface between the gray and white matter (Sparrey et al., 2009).

In particular, the viscoelastic properties extracted by Maikos and colleagues (2008b) from dynamic brain tissue tests by Mendis and associates (1995) seem to yield good results in simulating cord behavior in high-velocity impacts (with peak velocities of 0.489–0.690 m/sec for weight-drop contusions in the Maikos study, and 1 m/sec in the current study), but a more detailed characterization of rat cord viscoelasticity may be required to model its behavior accurately over a wider range of impact velocities. In fact, the Ogden and Prony hyperviscoelastic model we employed to simulate the cord properties may ultimately prove inadequate to accurately model cord behavior in all desired scenarios, and development of more complicated material models that include fully nonlinear viscoelasticity and a mechanism for dynamic tissue failure during simulation may be necessary. Because of the difficulty of accurately testing soft tissues like the cord to determine material properties, some combination of material testing and FE simulation may be required to achieve better characterization, as proposed by Morriss and colleagues (2008). Further testing considerations should also include careful attention to the influence of preconditioning on mechanical test results of the cord at various strain rates and magnitudes, an issue recently highlighted by Cheng and colleagues (2009).

A better understanding of the differences in injury patterns between mechanisms, such as that afforded by combined experimental and computational approaches, has the potential to influence future treatment and prevention approaches. For example, therapies might specifically target the central cavitation injury pattern and spared white matter rim associated with contusion injuries, or a full cord width dislocation injury with some rostral and caudal white matter damage, rather than attempting to treat more generic damage patterns. Better stratification of the patient population by factors such as injury mechanism may better identify the strengths and weaknesses of different SCI therapies, and lead to improved clinical trial outcomes overall. Furthermore, with continued improvements to the computational modeling of the spinal cord during traumatic injury, the cord may one day be added to full-body models used to aid design and safety testing of products ranging from automobiles to helmets, raising the profile of SCI in determining safety standards and promoting preventive strategies.

Footnotes

Acknowledgments

The authors thank Andrew Yung for assistance with MRI acquisition and Shivakumara Shetty and Paul Groenenboom for assistance with the initial set-up of smoothed particle hydrodynamics simulations. We are grateful to the Natural Sciences and Engineering Research Council of Canada (NSERC) for financial support.

Author Disclosure Statement

No competing financial interests exist.

References

Arbogast

K.B.

, Margulies

S.S.

1998. Material characterization of the brainstem from oscillatory shear tests. J. Biomech., 31:801–807.

Bain

A.C.

, Meaney

D.F.

2000. Tissue-level thresholds for axonal damage in an experimental model of central nervous system white matter injury. J. Biomech. Eng., 122:615–622.

Bilston

L.E.

, Thibault

L.E.

1996. The mechanical properties of the human cervical spinal cord in vitro. Ann. Biomed. Eng., 24:67–74.

Cheng

, Clarke

E.C.

, Bilston

L.E.

2008. Rheological properties of the tissues of the central nervous system: a review. Med. Eng. Phys., 30:1318–1337.

Cheng

, Clarke

E.C.

, Bilston

L.E.

2009. The effects of preconditioning strain on measured tissue properties. J. Biomech., 42:1360–1362.

Choo

A.M.

, Liu

, Dvorak

, Tetzlaff

, Oxland

T.R.

2008. Secondary pathology following contusion, dislocation, and distraction spinal cord injuries. Exp. Neurol., 212:490–506.

Choo

A.M.

, Liu

, Lam

C.K.

, Dvorak

, Tetzlaff

, Oxland

T.R.

2007. Contusion, dislocation, and distraction: primary hemorrhage and membrane permeability in distinct mechanisms of spinal cord injury. J. Neurosurg. Spine, 6:255–266.

Clarke

E.C.

, Choo

A.M.

, Liu

, Lam

C.K.

, Bilston

L.E.

, Tetzlaff

, Oxland

T.R.

2008. Anterior fracture-dislocation is more severe than lateral: a biomechanical and neuropathological comparison in rat thoracolumbar spine. J. Neurotrauma, 25:371–383.

Coats

, Margulies

S.S.

2006. Material properties of porcine parietal cortex. J. Biomech., 39:2521–2525.

10.

Cullen

D.K.

, LaPlaca

M.C.

2006. Neuronal response to high rate shear deformation depends on heterogeneity of the local strain field. J. Neurotrauma, 23:1304–1319.

11.

Farry

, Baxter

2010. The incidence and prevalence of spinal cord injury in Canada, technical reportRick Hansen Institute and Urban Futures.

12.

Fiford

R.J.

, Bilston

L.E.

2005. The mechanical properties of rat spinal cord in vitro. J. Biomech., 38:1509–1515.

13.

Fiford

R.J.

, Bilston

L.E.

, Waite

, Lu

2004. A vertebral dislocation model of spinal cord injury in rats. J. Neurotrauma, 21:451–458.

14.

Franconi

, Lemaire

, Marescaux

, Jallet

, Le Jeune

J.J.

2000. In vivo quantitative microimaging of rat spinal cord at 7T. Magn. Reson. Med., 44:893–898.

15.

Galbraith

J.A.

, Thibault

L.E.

, Matteson

D.R.

1993. Mechanical and electrical responses of the squid giant axon to simple elongation. J. Biomech. Eng., 115:13–22.

16.

Geddes-Klein

D.M.

, Schiffman

K.B.

, Meaney

D.F.

2006. Mechanisms and consequences of neuronal stretch injury in vitro differ with the model of trauma. J. Neurotrauma, 23:193–204.

17.

Greaves

C.Y.

, Gadala

M.S.

, Oxland

T.R.

2008. A three-dimensional finite element model of the cervical spine with spinal cord: an investigation of three injury mechanisms. Ann. Biomed. Eng., 36:396–405.

18.

Hyun

J.K.

, Kim

H.-W.

2010. Clinical and experimental advances in regeneration of spinal cord injury. J. Tissue Eng., 2010; 650857.

19.

Ichihara

, Taguchi

, Shimada

, Sakuramoto

, Kawano

, Kawai

2001. Gray matter of the bovine cervical spinal cord is mechanically more rigid and fragile than the white matter. J. Neurotrauma, 18:361–367.

20.

Jones

C.F.

, Kroeker

S.G.

, Cripton

P.A.

, Hall

R.M.

2008. The effect of cerebrospinal fluid on the biomechanics of spinal cord: an ex vivo bovine model using bovine and physical surrogate spinal cord. Spine, 33:E580–E588.

21.

Lee

K.E.

, Franklin

A.N.

, Davis

M.B.

, Winkelstein

B.A.

2006. Tensile cervical facet capsule ligament mechanics: failure and subfailure responses in the rat. J. Biomech., 39:1256–1264.

22.

X.-F.

, Dai

L.-Y.

2009. Three-dimensional finite element model of the cervical spinal cord: preliminary results of injury mechanism analysis. Spine, 34:1140–1147.

23.

Maikos

J.T.

, Shreiber

D.I.

2007. Immediate damage to the blood-spinal cord barrier due to mechanical trauma. J. Neurotrauma, 24:492–507.

24.

Maikos

J.T.

, Elias

R.A.I.

, Shreiber

D.I.

2008a. Mechanical properties of dura mater from the rat brain and spinal cord. J. Neurotrauma, 25:38–51.

25.

Maikos

J.T.

, Qian

, Metaxas

, Shreiber

D.I.

2008b. Finite element analysis of spinal cord injury in the rat. J. Neurotrauma, 25:795–816.

26.

Mendis

K.K.

, Stalnaker

R.L.

, Advani

S.H.

1995. A constitutive relationship for large deformation finite element modeling of brain tissue. J. Biomech. Eng., 117:279–285.

27.

Monaghan

1988. An introduction to SPH. Computer Phys. Commun., 48:89–96.

28.

Monaghan

1994. Simulating free surface flows with SPH. J Computational Phys., 110:399–399.

29.

Morriss

, Wittek

, Miller

2008. Compression testing of very soft biological tissues using semi-confined configuration—a word of caution. J. Biomech., 41:235–238.

30.

National Spinal Cord Injury Statistical Center. Spinal Cord Injury Facts and Figures at a Glance. Technical report, 2009.

31.

Persson

, McLure

S.W.D.

, Summers

, Hall

R.M.

2009. The effect of bone fragment size and cerebrospinal fluid on spinal cord deformation during trauma: an ex vivo study. J. Neurosurg. Spine, 10:315–323.

32.

Persson

, Summers

, Hall

R.M.

2011. The importance of fluid-structure interaction in spinal trauma models. J. Neurotrauma, 28:113–125.

33.

Piegl

, Tiller

1997. The NURBS Book. Springer Verlag: New York.

34.

Quinn

K.P.

, Winkelstein

B.A.

2007. Cervical facet capsular ligament yield defines the threshold for injury and persistent joint-mediated neck pain. J. Biomech., 40:2299–2306.

35.

Scifert

, Totoribe

, Goel

, Huntzinger

2002. Spinal cord mechanics during flexion and extension of the cervical spine: a finite element study. Pain Physician, 5:394–400.

36.

Shi

, Blight

A.R.

1996. Compression injury of mammalian spinal cord in vitro and the dynamics of action potential conduction failure. J. Neurophysiol., 76:1572–1580.

37.

Shreiber

, Bain

, Meaney

1997. In vivo thresholds for mechanical injury to the blood–brain barrier. 41st Stapp Car Crash Conference Proceedings, 277.

38.

Sparrey

C.J.

, Keaveny

T.M.

2011. Compression behavior of porcine spinal cord white matter. J. Biomech., 44:1078–1082.

39.

Sparrey

C.J.

, Manley

G.T.

, Keaveny

T.M.

2009. Effects of white, grey, and pia mater properties on tissue level stresses and strains in the compressed spinal cord. J. Neurotrauma, 26:585–595.

40.

Stys

P.K.

2004. White matter injury mechanisms. Curr. Mol. Med., 4:113–130.

41.

Sullivan

P.G.

, Rabchevsky

A.G.

, Waldmeier

P.C.

, Springer

J.E.

2005. Mitochondrial permeability transition in CNS trauma: cause or effect of neuronal cell death? J. Neurosci. Res., 79:231–239.

42.

Tator

C.H.

2006. Review of treatment trials in human spinal cord injury: issues, difficulties, and recommendations. Neurosurgery, 59:957–982discussion 982–987.

43.

Yoganandan

, Kumaresan

, Pintar

F.A.

2000. Geometric and mechanical properties of human cervical spine ligaments. J. Biomech. Eng., 122:623–629.

44.

Yushkevich

P.A.

, Piven

, Hazlett

H.C.

, Smith

R.G.

, Ho

, Gee

J.C.

, Gerig

2006. User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability. Neuroimage, 31:1116–1128.

45.

Zhu

, Prange

, Margulies

2006. Predicting unconsciousness from a pediatric brain injury threshold. Dev. Neurosci., 28:388–395.