As the strongest of the meningeal tissues, the spinal dura mater plays an important role in the overall behavior of the spinal cord-meningeal complex (SCM). It follows that the accumulation of damage affects the dura mater's ability to protect the cord from excessive mechanical loads. Unfortunately, current computational investigations of spinal cord injury (SCI) etiology typically do not include postyield behavior. Therefore, a more detailed description of the material behavior of the spinal dura mater, including characterization of damage accumulation, is required to comprehensively study SCI. Continuum mechanics-based viscoelastic damage theories have been previously applied to other biological tissues; however, the current work is the first to report damage accumulation modeling in a tissue of the SCM complex. Longitudinal (i.e., cranial-to-caudal long-axis) samples of ovine cervical dura mater were tensioned-to-failure at one of three strain rates (quasi-static, 0.05/s, and 0.3/s). The resulting stress–strain data were fit to a hyperelastic continuum damage model to characterize the strain-rate-dependent subfailure and failure behavior. The results show that the damage behavior of the fibrous and matrix components of the dura mater are strain-rate dependent, with distinct behaviors when exposed to strain rates above that experienced during normal voluntary neck motion suggesting the possible existence of a protective mechanism.
Based on tangent modulus measurements of uniaxial tensile tests on both human and bovine spinal samples [1–4], the dura mater is up to 100 times stiffer than the spinal cord and pia mater tissues that it encases. Accordingly, the dura mater plays an important stability role in the overall behavior of the spinal cord-meningeal complex (SCM). Together with the cerebrospinal fluid it contains, the dura mater functions to shield the weaker spinal cord from excessive mechanical loads and reduce cord deformation during traumatic loading scenarios such as vertebral burst fracture events [5–7]. Given its functional importance, multiple research groups have characterized the dura's mechanical properties under various loading conditions, including uniaxial tension at various strain rates [7–9], uniaxial tension in the longitudinal and circumferential directions [4,9,10], and quasi-static biaxial tension . Taken together, these studies have demonstrated that the dura mater exhibits a nonlinear viscoelastic response. However, unlike other viscoelastic soft tissues, it has been reported that the dura's mechanical response is not rate dependent at relatively low strain rates (i.e., 0.01–1 s−1) .
The material models developed from the aforementioned studies have provided properties that can be implemented into computational models of the SCM in order to improve our understanding and simulation of spinal cord injuries (SCI). Specifically, due to the complex loading that often leads to SCI, finite element (FE) models have traditionally been used as a robust means to determine the internal stresses and strains necessary to link the local mechanical environment to resulting tissue damage and/or long-term injury severity. These relationships between internal mechanical parameters and resultant damage are critical to the study, treatment, and prevention of SCI. Numerous research groups have developed FE models of the SCM for investigating the mechanical underpinnings of contusion (during vertebral burst fracture events), distraction, and dislocation events that lead to SCI [7,12–16]. Despite the known functional importance of the dura mater in protecting the spinal cord, most material models found in these FE models use simplifying assumptions of isotropic or anisotropic elasticity [7,14,15]. While such FE models have provided valuable insight into spinal mechanics and internal stress/strain patterns, many researchers in this area have concluded that additional characterization of the components of the SCM is necessary to accurately model the etiology of SCI [14–16], including mechanisms to simulate dynamic tissue damage . Therefore, a thorough description of the material behavior of spinal dura mater, including damage characterization during the subfailure and failure regimes, is necessary to comprehensively investigate SCI.
One method for incorporating post yield softening behavior of soft tissues into FE models is to implement continuum viscoelastic damage formulae [17–20]. Specifically, a variety of models have been used to describe the failure of vasculature [21–23]; tendon and ligament [24,25]; vaginal ; and rectus sheath  tissues. However, this technique has not yet been extended to describe the failure behavior of any tissue of the SCM. Therefore, the aim of this study is to characterize the damage accumulation behavior of the spinal dura mater under uniaxial loading. Additionally, while it has been reported that the hyperelastic parameters fitted to dura tension-to-failure tests are not strain-rate dependent , it is unknown whether the parameters used to describe damage accumulation are rate dependent. Therefore, a secondary objective of this study is to compare the damage accumulation process at three different strain rates.
Materials and Methods
Four ovine cervical spines (C2–C7) were collected from animals euthanized at Colorado State University's Surgical Research Laboratory for unrelated studies. Until mechanical testing was performed, spines were wrapped in saline-soaked gauze and stored at −20 °C. A single freeze-thaw cycle preservation technique has been shown to not significantly affect the mechanical properties of the dura mater [4,9] and other fibrous soft tissues [28–31]. On the day of testing, the spine was thawed at room temperature and the SCM was carefully removed via gross dissection, transection of the pedicles, and resection of the nerve roots. A longitudinal cranial-to-caudal incision was performed on the dura mater and the denticulate ligaments were severed to release the dura from the underlying meninges and cord.
As controversy still exists over the ultrastructural and mechanical differences between longitudinal and circumferential samples of spinal dura [4,7–9], a simple uniaxial tension-to-failure test in the dura's longitudinal direction (i.e., aligned with the long-axis of the SCM) was chosen as the procedure for modeling the damage accumulation in spinal dura. Longitudinal strips measuring 35 mm by 5 mm were cut from the sheet of dura mater, avoiding areas where exiting nerve roots would produce localized effects. This procedure yielded approximately 15 test specimens per spine that were subsequently randomized to one of the three strain rate groups. In order to keep the tissue adequately hydrated, samples waiting to be tested were placed in a saline bath at room temperature.
Uniaxial testing was performed using a custom-built test stand consisting of a 44.5-Newton load cell (model 31, Sensotec, Honeywell, Columbus, OH), a linear actuator with a 0.16 μm resolution (T-LLS105, Zaber Technologies, Vancouver, BC, Canada), two thin film grips (FC-20, IMADA, Northbrook, IL), and a 4.2-megapixel camera (Grasshopper3, Point Grey, Richmond, BC, Canada) positioned directly above the grips (Fig. 1(a)). To prevent slippage at the grip interface, 7.5 mm squares of water proof 180 grit sandpaper were attached with a small amount of cyanoacrylate to either end of each dural sample before it was placed in the grips. The sample was then loaded to a 0.5 N preload and digital images were acquired with the grips at 0 deg (the test orientation) and turned 90 deg (as shown in Fig. 1(b)). ImageJ (NIH, Bethesda, MD) was used to obtain five thickness measurements from the 90 deg image (Fig. 1(b)), five width measurements from the 0 deg image, and one length measurement (also from the 0 deg image). These measurements were used to determine the mean sample cross-sectional area and original length for post-hoc stress and strain rate calculations, respectively. Using displacement control of the linear actuator, each sample was tensioned to failure at 0.01 mm/s, 1 mm/s, or 6 mm/s (representing strain rates of 0.0005±7.27 × 10−5, 0.051±0.008, and 0.284±0.044 s−1, respectively). Figures 1(c)–1(e) show a sample at the preload, just prior to midsubstance failure and just after failure.
In order to follow the damage testing procedure outlined by Martins et al., the dural samples were not preconditioned prior to the pull-to-failure test . Reaction forces were recorded at 1000 Hz for the 1 mm/s and 6 mm/s speeds and at 100 Hz for the quasi-static speed (0.01 mm/s). Images were collected at 45 Hz for all tests (to monitor for grip slip and verify site of failure). Saline irrigation (1 drop per minute) was used to keep samples hydrated during the 0.01 mm/s test; the short duration of the 6 mm/s and 1 mm/s tests (<2 min) precluded the need for intratest hydration.
Only samples that failed midsubstance with no evidence of slippage were retained for analysis. The final sample sizes for the quasi-static, 1 mm/s, and 6 mm/s loading-rate groups were 14, 21, and 10, respectively. As global stretch measurements are needed for the modeling of failure behavior, actuator displacement was converted to global stretch using , where is the preload length and is the recorded actuator displacement. Cauchy stress was obtained using , where is the force (N) and is the cross-sectional area (mm2). Statistical analysis of dimensional measurements revealed no significant differences between groups, with an overall average initial length of 20.81±3.17 mm, width of 4.29±0.90 mm, and thickness of 0.18±0.04 mm.
where and is a dimensionless variable, and and are variables associated with the strain energy at which initial and total damage, respectively, occur during the loading procedure [17,32]. Therefore, the elastic response of the tissue is given by the behavior up to , while the damage accumulation properties are characterized by the behavior after the threshold is passed. Given the irreversible nature of the damage, a constraint is imposed that the damage evolution function must be monotonically increasing with , [26,32]. The terms are affected by the other fitting parameters and the size of the damage region of the curves .
In this formulation, it is assumed that the anisotropic fiber term only contributes to the total strain energy when the fibers are stretched (). The terms and characterize the location and length of the toe region of the response. In the previous equations, represent parameters that are analogous to the stiffnesses of the matrix and fibers, respectively, while are dimensionless parameters that characterize the matrix and fiber nonlinearity, respectively. It should be noted that are not independent as they ensure continuity of the strain field, the stress field, and the derivative of stress [26,27].
Fitting of the experimental data was performed with the fmincon function in matlab (version R2012b, The MathWorks, Natick, MA) via constrained nonlinear optimization of , , , , , , , , , , , and . In order to facilitate adequate fitting of the failure phase of the curves, which occurred much more rapidly than the ramp phase, the data were resampled to create a balance between the ramp and failure phases. Quality of the fits was assessed by root-mean-square error (RMSE) between the experimental Cauchy stress and stress given by the fitted damage model strain energy function. The results presented are the parameters that resulted in the lowest RMSE for each test.
Large standard deviations were noted for some of the fitted parameters. Accordingly, an outlier analysis was performed on each parameter group according to the method described by Moore and McCabe . Specifically, any datum point more than 1.5 times the interquartile range above the third quartile value or below the first quartile value was identified as an outlier. Tests containing any parameter identified as an outlier were excluded, resulting in final sample sizes of 9, 14, and 7 for the quasi-static, 1 mm/s, and 6 mm/s speed groups, respectively. These reduced data were used to determine statistical differences.
sigmaplot software (version 13.0, Systat Software, San Jose, CA) was used for all statistical analyses. To determine statistical differences between strain rate groups, Shapiro–Wilk and Brown–Foresythe tests were performed to test for normality and equal variance, respectively. Data that passed both tests were analyzed with a one-way ANOVA and post-hoc Student-Newman–Keuls tests. Data that failed either the normality or equal variance test were analyzed with a Kruskal–Wallis one-way ANOVA on ranks and post-hoc Dunn's tests. To determine statistical differences between related parameters within a strain rate group (e.g., comparing the stiffness of the matrix, , to that of the fibers, ), Shapiro–Wilk tests were performed for normality, followed by paired t-tests. A p-value of 0.05 was selected to define statistical significance.
All stretch–stress curves demonstrated strong nonlinearity, which is characteristic of hydrated collagen fiber-reinforced soft tissues. The majority of test data displayed a smooth increase in stress up to a maximum value, followed by a rapid decrease to zero. A small subset displayed minor decreases in stress during the ramp phase and/or plateaus in stress during the failure phase (Fig. 2). No significant differences were found between the three strain rate groups with respect to the maximum failure stress (p = 0.313) or the stretch at the maximum stress (p = 0.598). The data from the two highest speed groups fit the continuum damage model well (Fig. 2), with maximum RMSE values of 0.27 MPa (<1.5% of average maximum stress) and 0.47 MPa (<3% of average maximum stress) for the 6 mm/s and 1 mm/s loading-rate groups, respectively. The RMSE values for the quasi-static strain-rate fits were slightly higher, with a range of RMSE values from 0.36 to 1.48 MPa. The average RMSE of 0.95 MPa represents approximately 6% of the quasi-static strain-rate group's average maximum stress.
The average fitted parameters for the reduced data set (that excluding identified outliers) are shown in Table 1. Statistically significant differences were found between the strain-rate groups for (stiffness of the fibers), (nonlinearity of the fibers), (associated with the strain energy at the initiation of matrix damage), (associated with the strain energy at complete fiber damage), and (relates to size of the damage region of the fibers). Statistically significant differences within strain-rate groups were also found between the following parameter pairs: and (stiffness of the matrix and fibers, respectively); and (nonlinearity of matrix and fibers, respectively); and (associated with the strain energy at the initiation of damage in the matrix and fibers, respectively); and and (associated with the strain energy at complete matrix damage and fiber damage, respectively).
The average value of the 6 mm/s loading-rate group was significantly greater than that of quasi-static (p < 0.001) and the 1 mm/s (p = 0.003) groups, whereas the difference in between the 1 mm/s and the quasi-static groups was not statistically significant (p = 0.898). The parameter exhibited the opposite pattern; the average value of the 6 mm/s loading-rate group was significantly less than both the quasi-static (p < 0.001) and the 1 mm/s (p = 0.002) groups, while the difference between the 1 mm/s and quasi-static groups was not statistically significant (p = 0.99). When comparing the stiffness parameters and for within strain-rate groups, was significantly greater than for the 6 mm/s loading rate (p = 0.018). When comparing the nonlinearity parameters and within strain-rate groups, was significantly less than at all speeds (p = 0.001, p < 0.001, p = 0.003 for the quasi-static group, 1 mm/s loading-rate group, and 6 mm/s loading-rate group, respectively).
The value of (associated with the strain energy at the initiation of damage to the matrix) showed significant differences between the strain rate groups; the 6 mm/s loading-rate was significantly greater than that of both the quasi-static (p = 0.003) and 1 mm/s loading-rate (p = 0.006) groups. As with other parameters, the difference between the for the quasi-static and 1 mm/s loading-rate groups was not statically significant (p = 0.278). The value of (associated with the strain energy at complete fiber damage) followed the same pattern with respect to strain rate; the value for the 6 mm/s loading-rate group was significantly greater than that of both the quasi-static group (p = 0.007) and the 1 mm/s group (p = 0.020), but the difference between the quasi-static and 1 mm/s groups was not significant (p = 0.214). The minimum and maximum terms within a strain-rate group also showed significant differences. was significantly greater than at the quasi-static rate (p < 0.001), while was significantly greater than at all three rates (p = 0.005, p < 0.001, p = 0.019 for the quasi-static, 1 mm/s loading-rate, and 6 mm/s loading-rate groups, respectively).
The parameter, which describes the size of the damage region for the fibers, was significantly greater for the 1 mm/s loading-rate group than the quasi-static group (p = 0.008), but no other significant differences were found for this or the parameter. No significant differences in the and parameters, which relate to the nonlinear toe region of the response curve, were found between strain rate groups (p = 0.644 for and p = 0.911 for ). No other significant differences were found between strain rates groups or related parameters within strain rate groups.
Despite multiple studies in the literature on the tension-to-failure properties of spinal dura mater from both cadaveric and animal sources [4,8–10,35,36], this is the first study to apply a damage constitutive model to the dura and relate the model's parameters to strain-rate effects. These data and damage formulation can be implemented into finite element computational models of the SCM to improve the accuracy of simulations of spinal cord injury scenarios and dynamics.
As the continuum damage model used in this investigation is phenomenological in nature, the results obtained by application of the model (i.e., the damage parameters) are not easily interpreted. However, given the form used to model the anisotropic response of the fibers, (Eq. (9)), the term relates to fiber stiffness, while the term relates to the nonlinearity of the elastic fiber response. The significant differences found between strain-rate groups for and , therefore, follow the expectant strain-rate-dependent behavior of viscoelastic materials wherein a stiffer and more linearized response is characteristically obtained at higher strain rates [9,37–41]. Persson et al. fit uniaxial tensile tests of bovine dura matter at three strain rates to an Ogden model and reported an increase in G (i.e., stiffness) and a decrease in α (i.e., nonlinearity) with increasing strain rate for longitudinal samples, but the differences failed to reach statistical significance . Our data extended these results by explicitly finding that parameters associated with fiber stiffness and nonlinearity were significantly affected by strain rate.
The lack of significant strain-rate dependence with respect to the stiffness () and nonlinearity () of the isotropic matrix, (Eq. (8)), suggest that the matrix material is less rate-sensitive than the fibers, and that the differences seen in global properties are mainly due to the viscoelastic fiber response. The differences between (matrix stiffness) and (fiber stiffness) also support the conclusion that the fibers are primarily responsible for globally observed strain-rate dependency. At 6 mm/s, the stiffness of the fibers was, as expected, significantly higher than that of the matrix. However, at 1 mm/s, this difference was not significant and at the quasi-static loading-rate the stiffness of the fibers approximated that of the matrix. This pattern suggests that at very slow strain rates the matrix and fiber stiffness is nearly equivalent, but any increase in strain rate creates a nonproportionally greater increase in stiffness in the fibers. The significant difference between (matrix nonlinearity) and (fiber nonlinearity) was also as expected for a fiber-reinforced composite; while collagen fiber straightening (i.e., uncrimping) does contribute to composite nonlinearity, the majority of composite nonlinearity can usually be attributed to intrinsic properties of the matrix [42–45].
As the other eight fitted parameters are not directly relatable to the tissue's physical properties, it is slightly more difficult to draw conclusions about their significance. The two parameters related to the elastic response of the tissue not discussed earlier, and , were not significantly different between strain-rate groups. These parameters are associated with the stretch values that define the beginning and end of the curve's toe region. The lack of significance between strain-rate groups show that, despite differences in behavior within the toe region ( and ), the location and length of the toe region seems to be unaffected by strain rate. The values obtained for both parameters were similar to those reported in studies applying the same damage model to other soft tissues (rectus sheath: ; vaginal tissue: ). However, the toe region identified by and in the current study comprised a larger portion of the total stretch-stress curve than those found in the aforementioned studies. It is expected that the characteristics of the global toe region would depend not only on the nonlinearity of the matrix and fibers, but also on the distribution, density, and alignment of the fibers. Differences in these tissue-specific properties could account for the differences seen between studies.
The parameters related to the damage response of the tissue (namely, , , , , , and ) also showed significant differences between strain-rate groups and between related parameters within the same strain-rate group. Although the change in from the quasi-static to the 1 mm/s loading rate was not significant, the increase in across all three strain rates suggests that the initiation of damage to the matrix may be delayed at higher strain rates. (related to the initiation of damage to the fibers) showed the opposite trend with respect to strain rate, although none of the differences reached significance; the decrease in with increasing strain rate suggest that the fibers incur damage earlier in the loading process at higher strain rates. If a linear relationship between strain rate and the terms is assumed (Fig. 3), it is not surprising that the only significant difference between and was observed at the quasi-static loading rate. Interestingly, the slopes of the proposed linear relationships are equal in magnitude (to within two decimal places). According to these data, at relatively low strain rates the matrix incurs damage before the fibers, but at strain rates above 15%/s damage is initiated in the fibers before the matrix.
Examining the terms related to complete damage, both and increased with increasing strain rate, but the only differences that reached significance were between the of the 6 mm/s loading-rate group and that of the other two groups. In comparing the fiber and matrix terms, was significantly greater than at all three strain rates; this suggests that the matrix fails before the fibers such that the stress just prior to tissue failure is fully supported by the fibrous component. This pattern was also reported in damage model fits for rectus sheath  and vaginal tissue . Figure 4 shows these results in graphical form by plotting the damage parameters and (defined in Eq. (5) as functions of , , and , , , respectively) at all three strain rates. The initiation of damage relates to the deviation of the damage parameter from zero, while complete damage relates to the maximum value of the damage parameter (typically one). For a minority of the tests, especially those at the quasi-static strain rate, the load did not return completely to zero following midsubstance failure; for these cases, the damage parameter did not reach one and complete damage was taken as the maximum value obtained. In line with the strain-rate-dependent differences found for , the initiation of damage to the matrix occurs at higher stretch levels for the 6 mm/s loading-rate compared to the quasi-static and 1 mm/s loading rates. Also, in line with the strain-rate-dependent differences found for , complete fiber damage occurs at higher stretch levels for the 6 mm/s loading rate compared to the quasi-static and 1 mm/s loading rate.
Examining the difference between the initiation of damage and failure ( versus ) also provides insight into the failure process of each component. While there is a relatively large gap in and at the two lower strain rates, this gap is greatly reduced at the 6 mm/s loading rate (from an average of 0.48 MPa−1 to 0.08 MPa−1). This can also be seen in Fig. 4. The majority of matrix damage parameter curves at the quasi-static and 1 mm/s loading rate display an exponential shape, while all curves at the 6 mm/s loading rate display very steep slopes. This suggests that when increasing the loading rate from 1 mm/s to 6 mm/s, the matrix transitions from accumulating damage gradually to experiencing sudden failure with rapid damage accumulation. As for the gap between and , it appears to increase almost linearly as a function of strain rate from 0.01 MPa−1 at the quasi-static rate to 1.11 MPa−1 at the 6 mm/s (or 0.284 s−1) rate. In Fig. 4, this is shown as a gradual shift from very steep slopes to a more gradual accumulation of damage. All of the quasi-static curves are nearly vertical lines (representing a brittle-like behavior), whereas almost all of the curves for the 6 mm/s loading-rate group show an exponential shape. This suggests that at higher strain rates the fibers display a more ductile behavior, taking on more energy before failing.
With significant strain-rate effects between the damage behavior of both components at 6 mm/s compared to the other two speeds, these findings may represent a protective mechanism of the dura mater. During normal voluntary neck motion, the spinal tissues are exposed to strain rates between 0.04 and 0.24 s−1 . The strain rates in this study (0.0005, 0.051, and 0.284 s−1) were chosen to represent the quasi-static rate, a rate typically experienced during voluntary motion, and a rate slightly above the typical range. At the 0.284 s−1 rate, the initiation of damage to the matrix and the complete damage of the fibers were significantly delayed (occurring later in the loading cycle). This suggests that the tissue is able to support a higher strain energy without failing at higher loading rates. Therefore, in determining if tissue damage has occurred, the rate of loading (not only the maximum load or displacement reached) is significant. This distinction is important for any study of SCM injuries and suggests that computational models that use only static loading and/or static tissue damage properties are not capable of fully capturing the underlying tissue injury mechanisms.
There are limitations to the above work that should be noted. For this initial study, only the longitudinal loading direction was analyzed as it has greater implications for the injury mechanism of interest, i.e., whiplash (neck hyperflexion). In order to accurately model damage accumulation from circumferential loading, such as from excessive cerebrospinal fluid pressure, the aforementioned study should be repeated using the orthogonal loading direction. Similarly, based on structural and mechanical differences between spinal and cranial dura mater , it is unclear if the models developed above would be appropriate for use in cerebral investigations (e.g., models of traumatic brain injury or subdural hematoma). The strain rates utilized, while slightly above those seen in the tissue during normal voluntary motion, are still well below those reported for injurious levels [46–48]. Therefore, it is unknown if the differences seen between 1 mm/s and 6 mm/s can be extrapolated to higher speeds. Also, each group contained a relatively small sample size. It is possible that more significant differences would be found if additional samples were included in the analysis. Finally, as with most fitting procedures, there can be uncertainties regarding the identification of local versus global minimums and the uniqueness of the fitted results. This is especially true when the number of fitted parameters is relatively high. The application of an outlier analysis allowed for identification of fits that fell outside the typical range, and while this is an accepted practice when working with biological tissues [49–51], the differences identified between parameters may vary if the full data set had been analyzed instead of the reduced data set. However, of the nine statistically significant differences between strain-rate groups (as identified from the reduced data set), six remained significantly different (p < 0.05) or trended that way (p < 0.1), when the full data set was examined. This suggests that there are indeed differences between the groups that are obscured by one or two relatively large or small fitted coefficients. Again, additional samples could affect the outlier analysis or reduce the number of outliers.
Despite the above limitations, the results presented in the current study encourage additional work with the presented damage constitutive model. Future work will include testing of dura mater at higher strain rates, including those indicative of injury (up to 20 s−1), to explicitly investigate the role that protective mechanisms may play in dural damage accumulation. The small sample size precluded an analysis of anterior versus posterior differences, but there is evidence to suggest that the anterior dura mater may behave differently than the posterior dura mater . Specifically, since major spinal cord injury scenarios load the anterior and posterior aspects of the dura differently, future work will include a regional comparison of the damage accumulation process. Finally, the damage model will be applied to other tissues of the SCM, including the pia mater and the spinal cord (treating longitudinal axons as fibers).
In conclusion, this work is the first to report the application of a constitutive damage model to the dura mater. The results show distinct damage behaviors for the matrix and fiber constituents and that the damage effects vary with applied strain rate. These differences suggest a possible protective mechanism occurring at strain rates above what the tissue experiences during normal voluntary neck motion. Given these findings, it is imperative that the formulation presented herein be implemented into finite element computational models of the SCM in order to improve the accuracy of simulations of spinal cord dynamics and injury/damage scenarios.
National Institute of Biomedical Imaging and Bioengineering (Grant No. EB012048).