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.

Introduction

Based on tangent modulus measurements of uniaxial tensile tests on both human and bovine spinal samples [14], 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 [57]. 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 [79], uniaxial tension in the longitudinal and circumferential directions [4,9,10], and quasi-static biaxial tension [11]. 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) [9].

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,1216]. 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 [1416], including mechanisms to simulate dynamic tissue damage [16]. 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 [1720]. Specifically, a variety of models have been used to describe the failure of vasculature [2123]; tendon and ligament [24,25]; vaginal [26]; and rectus sheath [27] 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 [9], 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

Specimen Preparation.

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 [2831]. 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,79], 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.

Mechanical Testing.

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 [27]. 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 $λ=(L0+ΔL)/L0$, where $L0$ is the preload length and $ΔL$ is the recorded actuator displacement. Cauchy stress was obtained using $σ=λ*(F/A0)$, where $F$ is the force (N) and $A0$ 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.

Constitutive Modeling.

The nonlinear directional damage model for fibrous biological soft tissues first proposed by Calvo et al. [19] and further developed in 2008 [32] and 2012 [27] was used to fit the data in this study. This is an uncoupled directional damage model which defines different damage characteristics for the collagen and the ground substance material (matrix). Assuming that the damage process only affects the isochoric elastic part of the deformation [17], the free-energy density ($ψ$) of the continuum damage model can be expressed as follows:
$ψ(C,M,Dm,Df)=ψvol(J)+(1−Dm)ψ¯0m(C¯)+(1−Df)ψ¯0f(C,¯M)$
(1)
where $M$ is defined as the tensor product of the fiber direction vector in the undeformed configuration, $C¯$ is the modified right Cauchy–Green tensor (defined by $C¯=F¯TF¯$, where $F¯$ is the product of the deformation gradient tensor, $F$, and the cube root of its determinant, $J−13$), and $ψ¯0m$ and $ψ¯0f$ denote the elastic response of the undamaged matrix and collagen fibers, respectively [18,27]. The terms $(1−Dm)$ and $(1−Df)$ represent the reduction factors for damage as a function of the deformation (stretch) [17]; the damage parameters, $Dm$ and $Df$, are normalized [0,1] values related to how the matrix and fibers, respectively, accumulate damage. As proposed by Simo [17], the evolution of damage parameters are given by a series of piecewise and irreversible functions $Dk(Ξtk)$, with $k$ representing either the matrix ($m$) or the fibers ($f$). The interior functions, $Ξtk$, are defined as time ($s$) functions
$Ξsk=2ψ¯0k[C¯(s)]$
(2)
The maximum values of $Ξsk$ over the past time history $s∈(−∞,t)$ are defined as $Ξtk$
$Ξtk=maxs∈(−∞,t)2ψ¯0k[C¯(s)]$
(3)
The criterion for damage to occur is given by [17]
$Φk(C(t),Ξtk)=2ψ¯0k[C¯(t)]−Ξtk≤0$
(4)
meaning, at any time ($t$) of the loading procedure, that if the above equality is fulfilled for the matrix or fibers, then damage increases in that component of the tissue. The function describing the evolution of damage [$Dk(Ξtk)$] is given by
$Dk(Ξtk)={0,ifΞtk<Ξmin0kξ2[1−βk(ξ2−1)],ifΞmin0k≤Ξtk≤Ξmax0k1,ifΞtk>Ξmax0k$
(5)

where $ξ=[(Ξtk−Ξmin0k)/(Ξmax0k−Ξmin0k)]$ and is a dimensionless variable, and $Ξmin0k$ and $Ξmax0k$ 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 $Ξmin0k$, while the damage accumulation properties are characterized by the behavior after the $Ξmin0k$ threshold is passed. Given the irreversible nature of the damage, a constraint is imposed that the damage evolution function $Dk(Ξtk)$ must be monotonically increasing with $Ξtk$, $βk∈[−1.0,1.0]$ [26,32]. The $βk$ terms are affected by the other fitting parameters and the size of the damage region of the curves [27].

Fitting Procedure.

The spinal dura was assumed to be incompressible (i.e., $I3=J2=1$), which simplifies the free-energy density function given in Eq. (1) [27]
$ψ=(1−Dm)ψ̂0m(C)+(1−Df)ψ̂0f(C,M)$
(6)
For the case of uniaxial tension in the x3-direction, $λ3=λ$, $λ1=λ2=λ−12$ and $I4=λ2$. Therefore, the Cauchy stress tensor is diagonal with $σ33=σ$ and $σ11=σ22=0$. The total strain energy ($ψ$) is the superposition of the strain energies of the isotropic matrix ($ψ̂isom$) and the anisotropic collagen fibers ($ψ̂anif$)
$ψ=ψ̂isom+ψ̂anif$
(7)
The isotropic response of the dural matrix was modeled with the exponential strain energy function proposed by Demiray et al. [33], while the anisotropic response of the dural collagen fibers was modeled with the piecewise function proposed by Calvo et al. [26]
$ψ̂isom=c1c2{ exp[c22(I1−3)]−1}$
(8)

$ψ̂anif={0,ifI4I4ref$
(9)

In this formulation, it is assumed that the anisotropic fiber term only contributes to the total strain energy when the fibers are stretched ($I4>I40$). The terms $I40$ and $I4ref$ characterize the location and length of the toe region of the response. In the previous equations, $c1,c3>0$ represent parameters that are analogous to the stiffnesses of the matrix and fibers, respectively, while $c2,c4>0$ are dimensionless parameters that characterize the matrix and fiber nonlinearity, respectively. It should be noted that $c5,c6,c7$ 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 $c1$, $c2$, $c3$, $c4$, $I40$, $I4ref$, $Ξminm$, $Ξmaxm$, $Ξminf$, $Ξmaxf$, $βm$, and $βf$. 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 [34]. 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, $c1$, to that of the fibers, $c3$), Shapiro–Wilk tests were performed for normality, followed by paired t-tests. A p-value of 0.05 was selected to define statistical significance.

Results

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 $c3$ (stiffness of the fibers), $c4$ (nonlinearity of the fibers), $Ξminm$ (associated with the strain energy at the initiation of matrix damage), $Ξmaxf$ (associated with the strain energy at complete fiber damage), and $βf$ (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: $c1$ and $c3$ (stiffness of the matrix and fibers, respectively); $c2$ and $c4$ (nonlinearity of matrix and fibers, respectively); $Ξminm$ and $Ξminf$ (associated with the strain energy at the initiation of damage in the matrix and fibers, respectively); and $Ξmaxm$ and $Ξmaxf$ (associated with the strain energy at complete matrix damage and fiber damage, respectively).

The average $c3$ 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 $c3$ between the 1 mm/s and the quasi-static groups was not statistically significant (p = 0.898). The $c4$ parameter exhibited the opposite pattern; the average $c4$ 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 $c1$ and $c3$ for within strain-rate groups, $c3$ was significantly greater than $c1$ for the 6 mm/s loading rate (p = 0.018). When comparing the nonlinearity parameters $c2$ and $c4$ within strain-rate groups, $c4$ was significantly less than $c2$ 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 $Ξminm$ (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 $Ξminm$ 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 $Ξminm$ for the quasi-static and 1 mm/s loading-rate groups was not statically significant (p = 0.278). The value of $Ξmaxf$ (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. $Ξminf$ was significantly greater than $Ξminm$ at the quasi-static rate (p < 0.001), while $Ξmaxf$ was significantly greater than $Ξmaxm$ 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 $βf$ 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 $βm$ parameter. No significant differences in the $I40$ and $I4ref$ parameters, which relate to the nonlinear toe region of the response curve, were found between strain rate groups (p = 0.644 for $I40$ and p = 0.911 for $I4ref$). No other significant differences were found between strain rates groups or related parameters within strain rate groups.

Discussion

Despite multiple studies in the literature on the tension-to-failure properties of spinal dura mater from both cadaveric and animal sources [4,810,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, $ψ̂anif$ (Eq. (9)), the $c3$ term relates to fiber stiffness, while the $c4$ term relates to the nonlinearity of the elastic fiber response. The significant differences found between strain-rate groups for $c3$ and $c4$, 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,3741]. 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 [9]. 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 ($c1$) and nonlinearity ($c2$) of the isotropic matrix, $ψ̂isom$ (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 $c1$ (matrix stiffness) and $c3$ (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 $c2$ (matrix nonlinearity) and $c4$ (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 [4245].

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, $I40$ and $I4ref$, 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 ($c3$ and $c4$), 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: $I40=1.00−1.44$ [27]; vaginal tissue: $I40=1.05−1.31,I4ref=1.46−2.10$ [26]). However, the toe region identified by $I40$ and $I4ref$ 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, $Ξminm$, $Ξmaxm$, $Ξminf$, $Ξmaxf$, $βm$, and $βf$) also showed significant differences between strain-rate groups and between related parameters within the same strain-rate group. Although the change in $Ξminm$ from the quasi-static to the 1 mm/s loading rate was not significant, the increase in $Ξminm$ across all three strain rates suggests that the initiation of damage to the matrix may be delayed at higher strain rates. $Ξminf$ (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 $Ξminf$ 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 $Ξmin$ terms is assumed (Fig. 3), it is not surprising that the only significant difference between $Ξminm$ and $Ξminf$ 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 $Ξmaxm$ and $Ξmaxf$ increased with increasing strain rate, but the only differences that reached significance were between the $Ξmaxf$ of the 6 mm/s loading-rate group and that of the other two groups. In comparing the fiber and matrix terms, $Ξmaxf$ was significantly greater than $Ξmaxm$ 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 [27] and vaginal tissue [26]. Figure 4 shows these results in graphical form by plotting the damage parameters $Dm$ and $Df$ (defined in Eq. (5) as functions of $Ξminm$, $Ξmaxm$, $βm$ and $Ξminf$, $Ξmaxf$, $βf$, 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 $Ξminm$, 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 $Ξmaxf$, 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 ($Ξmin$ versus $Ξmax$) also provides insight into the failure process of each component. While there is a relatively large gap in $Ξminm$ and $Ξmaxm$ 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 $Ξminf$ and $Ξmaxf$, 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 $Df$ curves are nearly vertical lines (representing a brittle-like behavior), whereas almost all of the $Df$ 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 [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 [8], 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 [4648]. 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 [4951], 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 [10]. 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).

Conclusion

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.

Funding Data

• National Institute of Biomedical Imaging and Bioengineering (Grant No. EB012048).

References

References
1.
Bilston
,
L. E.
, and
Thibault
,
L. E.
,
1995
, “
The Mechanical Properties of the Human Cervical Spinal Cord In Vitro
,”
Ann. Biomed. Eng.
,
24
, pp.
67
74
.
2.
Mazuchowski
,
E. L.
, and
Thibault
,
L. E.
,
2003
, “Biomechanical Properties of the Human Spinal Cord and Pia Mater,” Summer Bioengineering Conference, Key Biscayne, FL, June 25–29, Paper No.
1205
.http://www.tulane.edu/~sbc2003/pdfdocs/1205.PDF
3.
Oakland
,
R. J.
,
Hall
,
R. M.
,
Wilcox
,
R. K.
, and
Barton
,
D. C.
,
2006
, “
,”
Proc. Inst. Mech. Eng. H.
,
220
(4), pp.
489
492
.
4.
Runza
,
M.
,
Pietrabissa
,
R.
,
Mantero
,
S.
,
Albani
,
A.
,
Quaglini
,
V.
, and
Contro
,
R.
,
1999
, “
Lumbar Dura Mater Biomechanics: Experimental Characterization and Scanning Electron Microscopy Observations
,”
Anesth. Analg.
,
88
(6), pp.
1317
1321
.
5.
Jones
,
C. F.
,
Kroeker
,
S. G.
,
Cripton
,
P. A.
, and
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
(
17
), pp.
E580
E588
.
6.
Hall
,
R. M.
,
Oakland
,
R. J.
,
Wilcox
,
R. K.
, and
Barton
,
D. C.
,
2006
, “
Spinal Cord-Fragment Interactions Following Burst Fracture: An In Vitro Model
,”
J. Neurosurg. Spine
,
5
(3), pp.
243
250
.
7.
Wilcox
,
R. K.
,
Allen
,
D. J.
,
Hall
,
R. M.
,
Limb
,
D.
,
Barton
,
D. C.
, and
Dickson
,
R. A.
,
2004
, “
A Dynamic Investigation of the Burst Fracture Process Using a Combined Experimental and Finite Element Approach
,”
Eur. Spine J.
,
13
(6), pp.
481
488
.
8.
Maikos
,
J. T.
,
Elias
,
R.
, and
Shreiber
,
D. I.
,
2008
, “
Mechanical Properties of Dura Mater From the Rat Brain and Spinal Cord
,”
J. Neurotrauma
,
25
(1), pp.
38
51
.
9.
,
C.
,
Evans
,
S.
,
Marsh
,
R.
,
Summers
,
J. L.
, and
Hall
,
R. M.
,
2010
, “
Poisson's Ratio and Strain Rate Dependency of the Constitutive Behavior of Spinal Dura Mater
,”
Ann. Biomed. Eng.
,
38
(
3
), pp.
975
983
.
10.
Mazgajczyk
,
E.
,
Ścigała
,
K.
,
Czyż
,
M.
,
Jarmundowicz
,
W.
, and
Będziński
,
R.
,
2012
, “
Mechanical Properties of Cervical Dura Mater
,”
Acta Bioeng. Biomech.
,
14
(
1
), pp.
51
58
.https://www.ncbi.nlm.nih.gov/pubmed/22742530
11.
Shetye
,
S. S.
,
Deault
,
M. M.
, and
Puttlitz
,
C. M.
,
2014
, “
Biaxial Response of Ovine Spinal Cord Dura Mater
,”
J. Mech. Behav. Biomed. Mater.
,
34
, pp.
146
153
.
12.
Li
,
X. F.
, and
Dai
,
L. Y.
,
2009
, “
Three-Dimensional Finite Element Model of the Cervical Spinal Cord: Preliminary Results of Injury Mechanism Analysis
,”
Spine
,
34
(
11
), pp.
1140
1147
.
13.
Sparrey
,
C. J.
,
Manley
,
G. T.
, and
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
(4), pp.
585
595
.
14.
Greaves
,
C. Y.
,
,
M. S.
, and
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
(
3
), pp.
396
405
.
15.
Scifert
,
J.
,
Totoribe
,
K.
,
Goel
,
V.
, and
Huntzinger
,
J.
,
2002
, “
Spinal Cord Mechanics During Flexion and Extension of the Cervical Spine: A Finite Element Study
,”
Pain Phys.
,
5
(
4
), pp.
394
400
16.
Russell
,
C. M.
,
Choo
,
A. M.
,
Tetzlaff
,
W.
,
Chung
,
T. E.
, and
Oxland
,
T. R.
,
2012
, “
Maximum Principal Strain Correlates With Spinal Cord Tissue Damage in Contusion and Dislocation Injuries in the Rat Cervical Spine
,”
J. Neurotrauma
,
29
(8), pp.
1574
1585
.
17.
Simo
,
J. C.
,
1987
, “
On a Fully Three-Dimensional Finite-Strain Viscoelastic Damage Model: Formulation and Computational Aspects
,”
Comput. Methods Appl. Mech. Eng.
,
60
(2), pp.
153
173
.
18.
Alastrué
,
V.
,
Rodríguez
,
J. F.
,
Calvo
,
B.
, and
Doblaré
,
M.
,
2007
, “
Structural Damage Models for Fibrous Biological Soft Tissues
,”
Int. J. Solids Struct.
,
44
(18–19), pp.
5894
5911
.
19.
Calvo
,
B.
,
Pena
,
E.
,
Martinez
,
M. A.
, and
Doblaré
,
M.
,
2007
, “
An Uncoupled Directional Damage Model for Fibred Biological Soft Tissues: Formulation and Computational Aspects
,”
Int. J. Numer. Methods Eng.
,
69
(10), pp.
2036
2057
.
20.
Rodríguez
,
J. F.
,
Cacho
,
F.
,
Bea
,
J. A.
, and
Doblaré
,
M.
,
2006
, “
A Stochastic-Structurally Based Three Dimensional Finite-Strain Damage Model for Fibrous Soft Tissue
,”
J. Mech. Phys. Solids
,
54
(4), pp.
864
886
.
21.
Li
,
D.
, and
Robertson
,
A. M.
,
2009
, “
A Structural Multi-Mechanism Constitutive Equation for Cerebral Arterial Tissue
,”
Int. J. Solids Struct.
,
46
(14–15), pp.
2920
2928
.
22.
Gasser
,
T. C.
,
2011
, “
An Irreversible Constitutive Model for Fibrous Soft Biological Tissue: A 3-D Microfiber Approach With Demonstrative Application to Abdominal Aortic Aneurysms
,”
Acta Biomater.
,
7
(6), pp.
2457
2466
.
23.
Balzani
,
D.
,
Schröder
,
J.
, and
Gross
,
D.
,
2006
, “
Simulation of Discontinuous Damage Incorporating Residual Stresses in Circumferentially Overstretched Atherosclerotic Arteries
,”
Acta Biomater.
,
2
(6), pp.
609
618
.
24.
Liao
,
H.
, and
Belkoff
,
S. M.
,
1999
, “
A Failure Model for Ligaments
,”
J. Biomech.
,
32
(2), pp.
183
188
.
25.
De Vita
,
R.
, and
Slaughter
,
W. S.
,
2007
, “
A Constitutive Law for the Failure Behavior of Medial Collateral Ligaments
,”
Biomech. Model. Mechanobiol.
,
6
(3), pp.
189
197
.
26.
Calvo
,
B.
,
Peña
,
E.
,
Martins
,
P.
,
Mascarenhas
,
T.
,
Doblaré
,
M.
,
Natal Jorge
,
R. M.
, and
Ferreira
,
A.
,
2009
, “
On Modelling Damage Process in Vaginal Tissue
,”
J. Biomech.
,
42
(5), pp.
642
651
.
27.
Martins
,
P.
,
Peña
,
E.
,
Natal Jorge
,
R. M.
,
Santos
,
A.
,
Santos
,
L.
,
Mascarenhas
,
T.
, and
Calvo
,
B.
,
2012
, “
Mechanical Characterization and Constitutive Modelling of the Damage Process in Rectus Sheath
,”
J. Mech. Behav. Biomed. Mater.
,
8
, pp.
111
122
.
28.
Woo
,
S. L. Y.
,
Orlando
,
C. A.
,
Camp
,
J. F.
, and
Akeson
,
W. H.
,
1986
, “
Effects of Postmortem Storage by Freezing on Ligament Tensile Behavior
,”
J. Biomech.
,
19
(
5
), pp.
399
404
.
29.
Huang
,
H.
,
Zhang
,
J.
,
Sun
,
K.
,
Zhang
,
X.
, and
Tian
,
S.
,
2011
, “
Effects of Repetitive Multiple Freeze-Thaw Cycles on the Biomechanical Properties of Human Flexor Digitorum Superficialis and Flexor Pollicis Longus Tendons
,”
Clin. Biomech.
,
26
(
4
), pp.
419
423
.
30.
Suto
,
K.
,
Urabe
,
K.
,
Naruse
,
K.
,
Uchida
,
K.
,
Matsuura
,
T.
,
Mikuni-Takagaki
,
Y.
,
Suto
,
M.
,
Nemoto
,
N.
,
Kamiya
,
K.
, and
Itoman
,
M.
,
2012
, “
Repeated Freeze-Thaw Cycles Reduce the Survival Rate of Osteocytes in Bone-Tendon Constructs Without Affecting the Mechanical Properties of Tendons
,”
Cell Tissue Bank.
,
13
(
1
), pp.
71
80
.
31.
Szarko
,
M.
,
Muldrew
,
K.
, and
Bertram
,
J. E.
,
2010
, “
Freeze-Thaw Treatment Effects on the Dynamic Mechanical Properties of Articular Cartilage
,”
BMC Musculoskelet. Disord.
,
11
(
1
), p.
231
.
32.
Peña
,
E.
,
Calvo
,
B.
,
Martinez
,
M. A.
, and
Doblaré
,
M.
,
2008
, “
On Finite-Strain Damage of Viscoelastic‐Fibred Materials: Application to Soft Biological Tissues
,”
Int. J. Numer. Methods Eng.
,
74
(7), pp.
1198
1218
.
33.
Demiray
,
H.
,
Weizsäcker
,
H. W.
,
Pascale
,
K.
, and
Erbay
,
H. A.
,
1988
, “
A Stress-Strain Relation for a Rat Abdominal Aorta
,”
J. Biomech.
,
21
(
5
), pp.
369
374
.
34.
Moore
,
D.
, and
McCabe
,
G. P.
,
1998
,
Introduction to the Practice of Statistics
,
W. H. Freeman
, New York.
35.
Patin
,
D. J.
,
Eckstein
,
E. C.
,
Harum
,
K.
, and
Pallares
,
V. S.
,
1993
, “
Anatomic and Biomechanical Properties of Human Lumbar Dura Mater
,”
Anesth. Analg.
,
76
(3), pp.
535–540
.https://www.ncbi.nlm.nih.gov/pubmed/8452262
36.
Zarzur
,
E.
,
1996
, “
Mechanical Properties of the Human Lumbar Dura Mater
,”
Arq. Neuropsiquiatr.
,
54
(
3
), pp.
455
460
.
37.
Rashid
,
B.
,
,
M.
, and
Gilchrist
,
M. D.
,
2014
, “
Mechanical Characterization of Brain Tissue in Tension at Dynamic Strain Rates
,”
J. Mech. Behav. Biomed. Mater.
,
33
, pp.
43
54
.
38.
Fiford
,
R. J.
, and
Bilston
,
L. E.
,
2005
, “
The Mechanical Properties of Rat Spinal Cord In Vitro
,”
J. Biomech.
,
38
(7), pp.
1509
1515
.
39.
Danto
,
M. I.
, and
Woo
,
S. L.
,
1993
, “
The Mechanical Properties of Skeletally Mature Rabbit Anterior Cruciate Ligament and Patellar Tendon Over a Range of Strain Rates
,”
J. Orthop. Res.
,
11
(1), pp.
58
67
.
40.
Pioletti
,
D. P.
,
Rakotomanana
,
L. R.
, and
Leyvraz
,
P. F.
,
1999
, “
Strain Rate Effect on the Mechanical Behavior of the Anterior Cruciate Ligament-Bone Complex
,”
Med. Eng. Phys.
,
21
(2), pp.
95
100
.
41.
Yamamoto
,
S.
,
Saito
,
A.
,
Kabayama
,
M.
,
Sugimoto
,
S.
,
Nagasaka
,
K.
,
Mizuno
,
K.
, and
Tanaka
,
E.
,
2003
, “Strain-Rate Dependence of Mechanical Failure Properties of Rabbit MCL and ACL,” Summer Bioengineering Conference, Key Biscayne, FL, June 25–29, Paper No.
0109
.http://www.tulane.edu/~sbc2003/pdfdocs/0109.PDF
42.
Hirsekorn
,
M.
, and
Grail
,
G.
,
2011
, “
On the Role of Matrix Nonlinearity in Mechanical Modeling of Long-Fiber Reinforced Composites
,”
International Conference on Composite Materials
, Jeju, South Korea, Aug. 21–26, Paper No.
F10-2-IF0552
.http://www.iccm-central.org/Proceedings/ICCM18proceedings/data/2.%20Oral%20Presentation/Aug26%28Friday%29/F10%20Multi-scale%20Modeling%20II/F10-2-IF0552.pdf
43.
Bocchieri
,
R. T.
, and
Schapery
,
R. A.
,
2000
, “
Nonlinear Viscoelastic Behavior of Rubber-Toughened Carbon and Glass/Epoxy Composites
,”
Time Dependent and Nonlinear Effects in Polymers and Composites
,
American Society for Testing and Materials
,
West Conshohocken, PA
, pp.
238
265
.
44.
Holmes
,
G. A.
,
Peterson
,
R. C.
,
Hunston
,
D. L.
,
McDonough
,
W. G.
, and
Schutte
,
C. L.
,
2000
, “
The Effect of Nonlinear Viscoelasticity on Interfacial Shear Strength Measurements
,”
Time Dependent and Nonlinear Effects in Polymers and Composites
,
American Society for Testing and Materials
,
West Conshohocken, PA
, pp.
98
117
.
45.
Modniks
,
J.
,
Joffe
,
R.
, and
Andersons
,
J.
,
2011
, “
Model of the Mechanical Response of Short Flax Fiber Reinforced Polymer Matrix Composites
,”
Procedia Eng.
,
10
, pp.
2016
2021
.
46.
Bilston
,
L. E.
,
1994
, “The Biomechanics of the Spinal Cord During Traumatic Spinal Cord Injury,”
Ph.D. dissertation
, University of Pennsylvania, Philadelphia, PA.http://repository.upenn.edu/dissertations/AAI9521002/
47.
McKenzie
,
J. A.
, and
Williams
,
J. F.
,
1971
, “
The Dynamic Behaviour of the Head and Cervical Spine During ‘Whiplash’
,”
J. Biomech.
,
4
(6), pp.
477
90
.
48.
Panjabi
,
M. M.
,
Cholewicki
,
J.
,
Nibu
,
K.
,
Babat
,
L. B.
, and
Dvorak
,
J.
,
1998
, “
Simulation of Whiplash Trauma Using Whole Cervical Spine Specimens
,”
Spine
,
23
(
1
), pp.
17
24
.
49.
Kural
,
M. H.
,
Cai
,
M.
,
Tang
,
D.
,
Gwyther
,
T.
,
Zheng
,
J.
, and
Billiar
,
K. L.
,
2012
, “
Planar Biaxial Characterization of Diseased Human Coronary and Carotid Arteries for Computational Modeling
,”
J. Biomech.
,
45
(5), pp.
790
798
.
50.
Wagner
,
H. P.
, and
Humphrey
,
J. D.
,
2011
, “
Differential Passive and Active Biaxial Mechanical Behaviors of Muscular and Elastic Arteries: Basilar Versus Common Carotid
,”
ASME J. Biomech. Eng.
,
133
(
5
), p.
051009
.
51.
Palevski
,
A.
,
Glaich
,
I.
,
Portnoy
,
S.
,
Linder-Ganz
,
E.
, and
Gefen
,
A.
,
2006
, “
Stress Relaxation of Porcine Gluteus Muscle Subjected to Sudden Transverse Deformation as Related to Pressure Sore Modeling
,”
ASME J. Biomech. Eng.
,
128
(5), pp.
782–787
.