Geometrical variations in white and gray matter affect the biomechanics of spinal cord injuries more than the arachnoid space

Abstract

Traumatic spinal cord contusions lead to loss of quality of life, but their pathomechanisms are not fully understood. Previous studies have underlined the contribution of the cerebrospinal fluid in spinal cord protection. However, it remains unclear how important the contribution of the cerebrospinal fluid is relative to other factors such as the white/gray matter ratio. A finite element model of the spinal cord and surrounding morphologic features was used to investigate the spinal cord contusion mechanisms, considering subarachnoid space and white/gray matter ratio. Two vertebral segments (T6 and L1) were impacted transversely at 4.5 m s⁻¹, which demonstrated three major results:

While the presence of cerebrospinal fluid plays a significant contributory role in spinal cord protection (compression percentage decreased by up to 19%), the arachnoid space variation along the spine appears to have a limited (3% compression decrease) impact.

Differences in the white and gray matter geometries from lumbar to thoracic spine levels decrease spinal cord compression by up to 14% at the thoracic level.

Stress distribution in the sagittal spinal cord section was consistent with central cord syndrome, and local stress concentration on the anterior part of the spinal cord being highly reduced by the presence of cerebrospinal fluid.

The use of a refined spinal cord finite element method showed that all the geometrical parameters are involved in the spinal cord contusion mechanisms. Hence, spinal cord injury criteria must be considered at each vertebral level.

Introduction

Thoracic and lumbar spinal cord injuries (SCIs), which account for 20% of the cases ¹ seen in traumatic injuries, carry a very high societal cost and have a great impact on the quality of life. ² Clinical trials, which aim to improve SCI patient rehabilitation, investigate new drugs and treatments which enhance neural tissue recovery.^3,4 However, no current treatment has thus far resulted in a positive and predictable outcome. ³ This underlines the need for further understanding of SCI, particularly contusion injuries following burst fracture, which show a high prevalence. ⁵ Both this prevalence and the severity of the injury highly vary, depending on the vertebral level of the injury and the type of vertebral fracture involved.^1,2,5 While this latter may be relatively well known today, the contribution of the geometrical variations in the human spinal cord and canal along the spine in SCI mechanisms is yet to be studied extensively. The geometrical parameters of particular interest are the spinal canal cross-sectional area, the amount of cerebrospinal fluid (CSF), and the white matter (WM)/gray matter (GM) proportion, since they are recognized to vary significantly along the spine.^6–9

Experimental studies have suggested that CSF plays an important role in spinal cord protection by reducing its deformation during impact, through its damping effects. ¹⁰ These findings have been confirmed by finite element (FE) studies, which have, in particular, allowed a description of internal stress and stress fields in the spinal cord during SCI.^11,12 Recent FE models have shown that a relationship between subarachnoid space diminution and spinal cord impairment increases following contusion-type loading.^12,13 However, Persson et al. ¹³ experimentally revealed low variations in spinal cord compression under transverse loading when considering subarachnoid spaces at the C6 or T12 human vertebral levels.

The spinal nervous system geometrical parameters considered in this study include subarachnoid space and WM and GM morphology in the transverse plane. The latter has been shown to vary along the spine,^6,7 and its influence on spinal cord contusion mechanisms has not been studied extensively. Yan et al. ¹⁴ have developed an FE model with realistic three-dimensional (3D) variations in the spinal cord to study SCI following burst fracture. However, they have so far only analyzed a single vertebral level.

The objective of this work was thus to develop a refined finite element method (FEM) allowing a description of SCI mechanisms, considering two representative vertebral levels. The T6 and L1 levels were chosen as they are levels with high burst fracture occurrence, ¹ and they offer geometrical variability, in terms of spinal cord occupation ratio within the spinal canal, as well as WM and GM ratio and geometry. Furthermore, the FE model used in this study was developed respecting the required characteristics identified through previous studies. These include a representation of several morphological features (GM and WM, pia and dura maters, denticulate ligaments, and CSF) and their variations along the cord, the use of nonlinear viscoelastic mechanical properties, as well as fluid–structure interactions.

Materials and methods

The geometry of WM and GM used in this study was based on histological cadaveric spinal cord cross sections obtained from the literature. ⁷ In order to obtain the sagittal profile of the spinal cord, each cross section was positioned in the spinal canal of a 50th percentile thoracic and lumbar spine (T1–L5) FE model (Spine Model for Safety and Surgery—SM2S)^15–17 using the correspondence between spinal and vertebral levels defined in the literature. ¹⁸ Pia mater was modeled as the external surface of the WM. The dura mater geometry was built by offsetting the pia mater’s elements. In this version, nerve roots were not included in the model, and holes in the dura mater corresponding to root sleeves were filled. Denticulate ligaments were attached along the lateral sides of the pia mater, and to the dura mater between each vertebral foramen. The WM and GM were meshed using solid tetrahedral elements and three-node shell elements were used for pia mater, dura mater, and denticulate ligaments. The elements’ characteristic lengths were chosen similar to the SM2S model for computational purposes, and decreased where needed, in order to allow a proper representation of all the geometric features (Table 1). Finally, a frictionless contact interface was defined between pia mater and dura mater to avoid interpenetration of the membranes.

Table 1.

Mesh characteristics and material properties of the spinal cord finite element model.

	Type of elements	Characteristic length (mm)	Number of elements		Material properties	Material parameters (units: g, mm, ms)
			Thoracic	Lumbar
Gray matter	Tetrahedral	0.8	38,293	34,531	Stress–strain tabulated	See Figure 1
White matter	Tetrahedral	1.2	90,763	56,389	Stress–strain tabulated	See Figure 1
Pia mater	Tria	1.2	17,290	11,212	Linear elastic	ρ = 0.001; ν = 0.45; E = 2.3
Dura mater	Tria	2	5246	3675	Linear elastic	ρ = 0.001; ν = 0.45; E = 5
CSF	Cubic	2	626,400	486,680	ALE	ρ = 0.001; Cl = 2089; νl = 8.9e − 4
Dentate ligaments	Tria	1.5	1042	866	Linear elastic	ρ = 0.001; ν = 0.45; E = 10

E: Young’s modulus; ρ: density; ν: Poisson’s modulus; Cl: liquid bulk modulus; νl: shear kinematic viscosity; CSF: cerebrospinal fluid; ALE: arbitrary Lagrangian–Eulerian.

Since the strain rate is known to influence the biomechanical behavior of the spinal cord, a strain rate–dependent tabulated law was used to represent WM and GM behaviors, as suggested by Sparrey et al. ¹⁹ As no dynamic mechanical properties of the spinal cord were reported in the literature, the tabulated law used in this work was built from a linear extrapolation of force–deformation curves reported in the literature;^20–25 hence, the 1 s⁻¹ stress versus strain curve was obtained by applying a 2.09 scale factor on the 5e⁻⁴ s⁻¹ available stress–strain curve (Figure 1).

Figure 1.

Stress–strain curves for white matter (circle marker) and gray matter (square marker) at 5e⁻⁴ s⁻¹ strain rate (dotted line) and 1 s⁻¹ (continuous line). Properties for strain rates between 5e⁻⁴ s⁻¹ and 1 s⁻¹ are interpolated by the solver.

The linear elastic properties were used for the pia mater, dura mater, and dentate ligaments, based on Young’s modulus and Poisson’s ratio obtained in the literature (Table 1). With CSF being composed of about 99% water, it was assumed to follow water properties. The viscosity of 0.89 MPa s⁻¹ used was included in the CSF viscosity range presented by Bloomfield et al. ²⁶

An arbitrary Lagrangian–Eulerian (ALE) formulation was used to reproduce the fluid–structure interaction between the CSF, the pia mater, and the dura mater. The ALE formulation, which has already been used in biomechanical studies, including fluid–structure interactions, ²⁷ allows a representation of the movement of the material (fluid) through a static mesh. All structures of the thoracic and lumbar nervous system were contained in a block of cubic ALE elements, and contact interfaces were defined between the fluid and solid membranes (pia mater and dura mater).

Contusions were simulated on a thoracic vertebral level (T5–T7 segment, impact at T6 vertebra) and a lumbar vertebral level (T12–L1 segment, impact at L1 vertebra). To ensure consistency with previous experimental works, and to support model verification, four different model definitions were used (Figure 2):

• Model 1: spinal cord (GM, WM, pia mater) only;

• Model 2: spinal cord and dura mater;

• Model 3: spinal cord, dura mater, and CSF (including dentate ligaments). In this configuration, the thickness of the subarachnoid space containing the CSF was defined to match the human morphology after positioning of dura mater relative to the inner surface of the SM2S spinal canal, with an epidural space of approximately 1.5 mm to account for the presence of surrounding soft tissues. The spinal cord/CSF antero-posterior diameter (∅AP) ratio measured in the resulting model was 0.42 at level T6 and 0.47 at level L1;

• Model 4: same as model 3, but with a subarachnoid space thickness corresponding to bovine morphology in order to allow a comparison with the reference experimental study from Persson et al. ¹⁰ The spinal cord/CSF ∅AP ratio was manually adjusted here to 0.78 at both the thoracic and lumbar vertebral levels.

Figure 2.

Finite element model of the spinal cord: (a) sagittal section of the complete model (spinal cord in red, dura matter in yellow), with transverse cross section of the spinal cord for both vertebral levels (T6 and L1) tested in this study and (b) representation of the different models used for the sensitivity study.

A fixed rigid wall representing the posterior elements of the vertebrae was modeled posterior to the spinal cord models. Horizontal cylindrical impactors of three different cross-sectional areas (impactor 1: 314 mm², impactor 2: 167 mm², and impactor 3: 83.5 mm²), with the same weight (7 g), were used to simulate an anterior contusion on the center of the spinal cord, at an initial velocity of 4.5 m s⁻¹, as described in experimental investigations. ¹⁰

All the simulations were performed with an explicit solver (Radioss v.10.0; Altair Engineering Inc., Troy, MI, USA), on a 6-ms time range. The spinal cord AP diameter was measured at maximum compression and divided by the original diameter to calculate the compression percentage (CP). The compressive force at maximum compression was also measured. The Spearman rank test was performed to evaluate the correlation between the CP and Von Mises stresses in the WM and GM. Wilcoxon matched paired tests were used to compare the CP results to reference data from the literature and to compare the CP and Von Mises stresses between the T6 and L1 vertebral levels and between models 3 and 4. Mann–Whitney tests were used to compare the CP and Von Mises stresses between the simulations, with and without fluid. Finally, the biomechanics of SCI was described through stress distribution and impactor speed and displacement in time observed in the sagittal medial plane.

Results

Spinal cord contusion mechanisms

Spinal cord CP and maximum Von Mises stress results are presented in Tables 2 and 3, respectively. The similarities of the CPs with experimental data ¹⁰ were a function of the vertebral level. CP values were included in the experimental range for 89% of the simulations at the T6 level (p = 0.28) and for only 33% at the L1 level (p = 0.007). Impactor 3 led to the highest CPs, followed by impactor 2, and then impactor 1. Considering all levels, maximum Von Mises stresses ranged between 0.33 and 0.96 MPa in the WM and between 0.45 and 1.81 MPa in the GM (Table 3) and were significantly correlated with CP values (p < 0.001, r = 0.68 for GM; p < 0.001, r = 0.74 for WM). In all the simulations, maximum stresses were higher in GM compared to WM and were located in anterior areas of WM and GM, as can be seen in Figure 3.

Table 2.

Compression of the spinal cord in antero-posterior direction grouped by models, impactor type, and vertebral level.

Compression percentage		Level T6			Level L1			Mean
		Impactor			Impactor
		1	2	3	1	2	3
Model	1	0.52	0.59	0.67	0.38	0.53	0.59	0.55
	2	0.52	0.59	0.66	0.43	0.52	0.58	0.55
	3	0.34	0.42	0.49	0.28	0.35	0.55	0.41
	4	0.39	0.47	0.54	0.30	0.38	0.57	0.44
Mean		0.44	0.52	0.59	0.35	0.44	0.57
		0.52			0.45

Table 3.

Maximum Von Mises stresses in WM and GM grouped by models, impactor type, and vertebral level.

Maximum Von Mises stress (MPa) (GM; WM)		Level T6			Level L1			Mean
		Impactor			Impactor
		1	2	3	1	2	3
Model	1	1.47; 0.84	1.58; 0.84	1.54; 0.84	0.91; 0.70	1.29; 0.87	1.25; 0.84	1.34; 0.82
	2	1.30; 0.89	1.81; 0.79	1.60; 0.89	1.06; 0.81	1.33; 0.91	1.34; 0.96	1.41; 0.87
	3	0.57; 0.34	0.77; 0.46	0.91; 0.59	0.45; 0.33	0.71; 0.54	0.85; 0.64	0.71; 0.48
	4	0.77; 0.47	0.90; 0.56	1.10; 0.71	0.46; 0.45	0.67; 0.52	0.84; 0.64	0.79; 0.56
Mean		1.03; 0.63	1.26; 0.66	1.29; 0.76	0.72; 0.57	1.00; 0.71	1.07; 0.77
		1.19; 0.68			0.93; 0.68

GM: gray matter; WM: white matter.

Figure 3.

Typical Von Mises stress (MPa) field in (a) gray matter and (b) white matter. This figure was obtained at maximum cord compression for configuration 3, impactor 3, at L1 level.

For models 1 and 2, injury kinematics included a single compression, with a decrease in the impactor speed starting at the first contact with the spinal cord. As shown in Figure 4, there are two phases of spinal cord compression when including the CSF (models 3 and 4). The first phase starts at the first contact between the impactor and the dura mater (Figure 4b) and ends with the contact between the impactor, dura mater, cord, and posterior wall (no CSF is left at the impactor level (Figure 4c)). The second phase ends at the maximum cord compression (Figure 4d), after which the impactor is pushed back by the spinal cord’s elasticity.

Figure 4.

Typical change in spinal cord antero-posterior (AP) diameter through time during compression (lumbar spinal cord segment, impactor 3). Von Mises stresses are seen on transverse cross sections of the spinal cord (a) at initial state, (b) at dura matter contact, (c) at spinal cord contact, (d) at maximum compression, and (e) after impactor’s rebound.

Relative influence of morphological features: WM versus GM ratio and CSF influence

The CP was significantly different (paired T-test on CP for each level, grouped by impactor, p < 0.01) between the thoracic and lumbar levels and was 7% overall, and up to 14% (model 1, impactor 1) higher at T6 compared to the L1 vertebral level. The maximum compressive force during impact stresses was not significantly different between the two vertebral levels for WM, but were higher at T6 compared to L1 for GM (p = 0.002).

The presence of CSF significantly decreased the compressive force, from 55.2 ± 11.8 N (models 1 and 2, both levels) to 34.9 ± 8.6 N (models 3 and 4, both levels). For both vertebral levels, the compressive force and CP were significantly higher without CSF, surrounding the spinal cord (p < 0.001). The compressive force was 55.2 ± 11.8 N for all models without fluid and 34.9 ± 8.6 N with fluid. Models 1 and 2 obtained the highest CP, followed by model 4 and model 3. The maximum stress was significantly decreased by the presence of CSF (comparison of models 2 and 3) for both WM (49% decrease) and GM (45% decrease; p < 0.001). However, even where it was significantly different, CP was only 3% higher in model 4 compared to model 3. However, this slight overall increase in CP resulted in a stress increase of up to 28% (WM, impactor 1, level T6).

Discussion

This article presents an original study of the biomechanics of spinal cord contusion, using a detailed FEM of the spinal cord, with a realistic representation of the spinal cord morphology and behavior under traumatic loading. SCIs following thoracic and lumbar burst fractures were reproduced, including biofidelic viscoelastic material behavior for the WM and GM and fluid–structure interaction between the CSF and the spinal cord and canal. Levels T6 and L1 were chosen since they present a high burst fracture occurrence, ¹ in addition to a distinct WM/GM geometry allowing an investigation of the effect of cross-sectional geometry at the thoracic and lumbar vertebral levels, as described in the literature^6,7 (about 13% of the total cord occupied by the GM at T6 vertebra, 41% at L1). For verification purposes, the different models that were used were defined similar to a previous experimental study. ²⁸

This study has confirmed previous findings^12,13 concerning the role played by CSF presence in reducing the spinal cord CP and Von Mises stresses in WM and GM. However, the results of this study show that slight variations in the subarachnoid space thickness, such as between bovine and human configurations or between human vertebral levels, have a slight impact on spinal cord compression. These findings suggest that the spinal cord is indeed protected by fluid, and the protection is subject to low variations when the subarachnoid space exceeds a minimum threshold. When the subarachnoid space is situated below this threshold, such as in murine spine, ²⁹ fluid protection is insignificant. The minimum subarachnoid space threshold for spinal cord protection should be defined in future studies. This result highlights the importance of developing a human or an animal FEM of the spine relative to species-specific characteristics.

An analysis of injury kinematics showed that when the subarachnoid space thickness is sufficient, CSF leakage in the first phase of compression initiates spinal cord movement within the spinal canal, significantly reducing the maximum Von Mises stress in both WM and GM. Clinically, these results support the importance of spinal cord decompression in the case of spinal stenosis for preventing secondary neurologic trauma.

The results also suggest that spinal cord CP is significantly different between the thoracic and lumbar vertebral levels. This difference did not seem to originate from the variations in CSF thickness, as the thoracic spinal cord/CSF ∅AP ratio was only 5% lower than at the lumbar level. The dependency was most likely due to a higher stiffness in GM compared to WM. With the cross-sectional surface ratio of GM over WM being higher at the L1 level than at the T6 level, the lumbar segment was globally stiffer than the thoracic segment. The behavior at the T6 level was close to the cervical spinal cord behavior reported by Persson et al. ¹⁰ which is consistent with morphometric studies ⁷ reporting that the cross-sectional surface ratio of GM over WM at the thoracic levels (T6: 13.4%) is closer to what is seen at the cervical levels (C4: 14.5%) than at the lumbar levels (L1: 41.1%). This result shows the importance of accurate modeling of WM and GM, in terms of geometry and mechanical properties’ assignment. In addition, previously, the severity of SCI has mostly been linked with spinal level-specific external loading characteristics. The results also indicate that the geometric properties of the spinal cord must also be taken into account when analyzing burst fracture resulting in neural injury. From a clinical perspective, no stipulations can be made regarding post-traumatic radiographs, as they do not provide information on bone fragment displacements during the fracture occurrence. However, for a given spinal cord deformation, the results of this study suggest that at the lumbar spinal levels, the GM cell bodies could undergo more stress, compared to what occurs at the thoracic and cervical levels.

The maximum Von Mises stresses in WM and GM were higher than in other FE studies^30,31 due to the higher impact forces and deformations simulated in this study. Previous studies have shown that GM is stiffer and more fragile than WM. ²³ Consequently, the stress distribution in the spinal cord was consistent with central cord syndrome. Also, the results showed that the Von Mises stresses in the GM varied with the vertebral level, indicating an influence of either GM geometry variations or the higher CPs at the thoracic levels. On the other hand, it appeared that the maximum Von Mises stress in the WM was mainly caused by contact with the impactor and was not affected by the vertebral level (0.71 MPa in T6 vs 0.74 in L1 on average). This result shows that the presence of CSF decreases stress concentrations resulting from the bone fragment impacting the cord by limiting direct contact with the anterior part of the spinal cord. This particular point has an impact on the protection of motor pathways, which are predominant in the anterior funiculi. This also shows that CSF should not be modeled with solid elements, which prevent direct contact between the impactor, dura mater, and pia mater; rather, specific computational fluid dynamics (CFD) formulations, such as ALE or smoothed particle hydrodynamics, should be preferred.

Although the results are in agreement with the experimental results, FE simulation results need to be analyzed with care. Spinal cord geometry was taken from post-mortem histology data ⁷ and has been shown to be a good representation of in vivo geometry, even though post-mortem spinal cord dimensions may be slightly smaller.^6,7 Moreover, geometrical characteristics, such as spinal cord eccentricity in the spinal canal, nerve roots anchoring, arachnoid matter presence, and spinal cord anterior median fissure mechanical characteristics, should be considered when developing a spinal cord FEM. Also, the spinal cord anisotropic properties were not represented in this study. Considering them would allow a more phenomenological analysis of SCI. Future works should focus on a comprehensive representation of the spinal cord tissue anisotropy and hence provide the corresponding anisotropic stress and strain analysis. Despite these limitations, the model was found to be suited to relative comparisons of different impact conditions. The impactors were defined to obtain different contact surfaces, although they are not representative of the variations in bone fragment found after burst fractures, as they do not vary in shape and weight. However, the results confirm that a higher contact surface between a given bone fragment and the spinal cord will decrease the contact pressure and thus the cord’s deformation. An explicit solver was used to take into account inertia and damping effects, as well as to facilitate the management of contact and material non-linearities. In order to ensure simulation accuracy, an energy-based error criterion was observed and did not exceed 10%. Tetrahedral elements were selected, as they allowed a representation of small and irregular geometries, as well as the variations in the geometry of WM and GM over all the vertebral levels. The use of tetrahedral elements was combined with the calibration of dedicated material properties, allowing a proper fit with reference experimental results. However, tetrahedral elements are known to increase the stiffness of a structure under great strains. This phenomenon probably prevented the model from achieving spinal cord compressions higher than 70%. Also, a mesh convergence study (see supplementary data file) revealed a stable behavior in terms of compression percent and stress levels when modifying the elements’ characteristic length (from 1.5 to 0.25 mm). The material properties of the spinal cord were taken from in vitro experiments reported in the literature,^20–25 which could have been altered by post-mortem degradation. ³² Therefore, further experimental studies should be undertaken, ideally on fresh specimens, in order to define viscoelastic properties for the spinal cord under dynamic loading. In future developments of the model, the use of poroelastic material should be investigated in order to take into account the spinal cord’s vascularization, and the implementation of a damage criterion for the dura matter would allow the consideration of dural tears.

Conclusion

The behavior of the FE model of the central nervous system developed and used in this study was consistent with reference data and allowed spinal cord contusion analysis at two representative vertebral levels. The results showed that CSF presence is important for spinal cord protection, but that the variations in the subarachnoid space have little impact on SCI biomechanics. In addition, the kinematics of spinal cord contusion was described as vertebral level dependent due to the geometrical variations in WM and GM.

This work demonstrates the importance of geometrical characteristics such as subarachnoid space and WM and GM morphologies when developing a representative model for SCI simulation. The results support the importance of rapidly reducing spinal canal stenosis in order to prevent secondary cord injury. From a computational perspective, this study has shown the importance of modeling the CSF as a fluid, with non-solid elements. Here, the ALE formulation allows the fluid to completely flow out of the impact area and reveals high stresses in WM caused by contact between the impactor, dura mater, and pia mater. Further analyses should be carried out in order to describe other SCI mechanisms, such as rotational injuries (type C in Magerl’s vertebral fractures classification ³³ ), which are an important cause of SCI, but have rarely been studied.