Rate- and Region-Dependent Mechanical Properties of Göttingen Minipig Brain Tissue in Simple Shear and Unconfined Compression

Traumatic brain injury (TBI) is a major cause of injury and death. In the United States (U.S.), TBI is responsible for over two million emergency department visits per year and contributes to the deaths of over 50,000 people, with the incidence of TBI-related emergency department visits in the U.S. increasing in recent years [1]. TBI is also an increasing problem in modern military conflicts, with an estimated 300,000 of the 1.6 million U.S. service members deployed to Iraq and Afghanistan sustaining at least one mild TBI between 2001 and 2007 [2].

Finite element (FE) simulations are an important tool to understand the complex mechanics of TBI [3,4], but the selection and definition of their material models significantly influence their predictions [5]. Clearly, the use of constitutive models not validated for strains and strain rates relevant to TBI may introduce error [6]. Experimental characterization of brain tissue began in the late 1960s. Experiments since have examined a variety of loading modes, directions, and rates; brain regions; and animal species. Properties for some combinations of these variables still need characterization, but it is notable that comparable experimental studies show wide variation in material properties [7]. This variation is evidence of the difficulty of characterizing brain tissue, likely because of its high water content [8,9], relatively soft structure (stiffness of the order of 100–1000 Pa [10]), and strong dependence on strain rate [11–13].

Explosive blast produces complex loading modes in the brain, causing both increased intracranial pressure and shearing [14]. As a result, constitutive models used to study blast TBI should consider both compressive and shear properties of brain tissue. Further, while many investigations of brain tissue include evaluation at high strain rates, most do not include rates relevant to blast exposure. Strain rates in finite element models of human blast events have been reported of the order of ten [15,16] to several hundred per second [17] depending on the model used. Strain rate data from blast simulations on porcine brains is lacking in the literature, but due to the smaller size of the porcine brain it is reasonable to expect rates to be larger than those seen in human brain simulations. Modeling of impact injuries in porcine and human brains has shown that peak strain rates are indeed higher in pigs than humans [18], with peak strain rates being an order of magnitude higher in pigs (up to 235 s⁻¹) compared to humans (up to 65 s⁻¹). Given this, there is a need for the development of validated constitutive models at high strain rates based on experimental data from both shear and compression to aid in modeling both impact and blast injuries.

Few authors have examined the response of brain tissue ramps in simple shear, and most of these studied the tissue at strain rates below 50 s⁻¹ ([19–23]). A more substantial body of literature reports the properties of brain tissue in unconfined compression, but like shear, many of these papers focus on lower rates (below 50 s⁻¹) or quasi-static (QS) compression [21–31]. To date, no author has performed the experiments necessary to fully define and model the rate dependent response of brain tissue in both simple shear and unconfined compression over the range of strain rates relevant to blast injury.

Most studies on brain tissue focus on tissue harvested from the cerebrum. Few authors have examined tissue from the brainstem, despite both its physiological importance and its propensity for injury in TBI [32]. The majority of the published brainstem studies tested the tissue in simple shear [22,23,33–35]. No compression data exist for the brainstem above a rate of 50 s⁻¹. In addition, few authors have investigated the material properties of the cerebellum. While the cerebellum is broadly a mix of white and gray matter, much like the cerebrum, material properties measured by magnetic resonance elastography show that the cerebellum is softer than the cerebrum [36], suggesting the need to model it as a separate material. Other authors have examined murine cerebellar tissue exposed to high-rate shear [35] and compression [23,25] as well as porcine and murine tissue using low-rate micro-indentation [37] (producing boundary conditions similar to confined compression [10]). No high-rate data exists for porcine cerebellar tissue in either shear or unconfined compression at strain rates above 50 s⁻¹.

In this study, we deformed brain tissue from the cerebrum, cerebellum, and brainstem of adolescent male Göttingen minipigs in shear and compression at a quasi-static strain rate of 0.02 s⁻¹, an intermediate strain rate of 150 s⁻¹, and a high strain rate of 300 s⁻¹. While strain rates in porcine blast may be higher than the 300 s⁻¹ rate, we chose this rate as it was on the upper end of the velocities achievable by our test equipment while allowing for larger tissue samples that minimize the influence of boundary conditions. Tests were conducted with the objective of defining a constitutive model that would be appropriate for blast simulations in these animals. Mooney–Rivlin and both first and second-order Ogden hyperelastic models, each coupled with a one-term Prony series to capture rate dependence, were evaluated as potential constitutive representations of the measured material response. Taken together, this research provides a set of experimentally fit constitutive models appropriate over a wide range of strain rates that can be used in computational simulations of TBI [3].

Thirty-five juvenile (aged 4–5 months) male Göttingen minipigs were euthanized via intravenous phenytoin and pentobarbital. Both the Animal Care and Use Review Office of the U.S. Army Medical Research and Development Command, Fort Detrick, MD, and the Institutional Animal Care and Use Committee at the University of Utah approved all experimental protocols. The braincase was carefully opened with a hammer and chisel. The brain was then freed from the dura and cranial nerves, removed from the skull, and placed in a PEG buffer (7.5% by weight polyethylene glycol in phosphate buffered saline) to minimize swelling relative to the native tissue state [38]. To minimize the potential influence of tissue temperature [31,39], all tissue was stored and tested at room temperature (21 °C).

Tissue samples were harvested from the cerebrum, cerebellum, and brainstem immediately prior to mechanical testing (completed within 8 h of death). The brain was first sectioned into thick (roughly 10 mm) coronal slices and then cut using a custom square die having inner dimensions of 8 mm × 8 mm, with a height of 7 mm for cerebrum and cerebellum samples and 5 mm for brainstem samples. A scalpel was run over the surface of the die to ensure that the surface of the specimen was uniformly cut to the correct height.

Samples from the brainstem (Fig. 1(a)) were harvested from the medulla and cut so that the direction of shear was applied in the rostral-caudal direction and compression was applied in the inferior-superior direction. Samples from the cerebellum (Fig. 1(b)) were cut so that the direction of shear was applied in the inferior-superior direction and compression applied in the rostral-caudal direction. Samples from the cerebrum (Fig. 1(c)) were cut from the corona radiata so that the direction of shear was applied in the inferior-superior direction and compression applied in the rostral-caudal direction. For each brain, 2 to 3 samples were harvested from both the cerebellum and brainstem, and 2 to 4 samples were harvested from the cerebrum, for a total of n = 8 samples per group. While the relatively large size of the human brain allows investigators to examine isolated white and gray matter samples [19], studies using smaller animal brains generally test mixed white and gray matter samples [11,23–26,28–30,35,40–49] due to size constraints.

Tissue testing was performed on a custom soft tissue tester previously used by our group [50] to test cerebral blood vessels at high strain rates (over 1000 s⁻¹). Tissue was mounted (Fig. 2, top row) to produce either shear or compression between the sled and either a 250-gram load cell for shear tests (model 31 Low, Honeywell, Golden Valley, MN) or a 1000-gram load cell for compression tests (model 31 Mid, Honeywell, Golden Valley, MN). Displacement was measured either via a hall effect sensor built into a voice coil (MGV52-25-1.0, Akribis Systems, Singapore) for QS tests or a linear variable differential transformer (LVDT, DC-SE 2000, TE Connectivity, Schaffhausen, Switzerland) and a high-speed camera (Phantom Micro-EX4; Vision Research, Wayne, NJ) for high-rate tests. Sensor data were acquired using a DAQ card capable of simultaneous sampling (PXI-6133, National Instruments, Austin, TX) controlled by a custom program in LABVIEW (National Instruments, Austin, TX).

Using a voice coil actuator directly connected to a low-friction sled, the tester was capable of QS motions, which was directly controllable via a labview VI. High-rate deformations were achieved by replacing the voice coil with a drop tube (Fig. 2, bottom row) connected to the sled via steel cabling attached to a foam-padded can. Velocities corresponding to strain rates of 150 and 300 s⁻¹ were consistently achieved by dropping a mass into the can from a predefined height.

Samples were mounted on the tester after securing the sled in place with a set screw. Prior to mounting, all samples were measured with calipers and photographed. For shear tests, the samples were glued between two parallel shear plates with a thin layer of cyanoacrylate adhesive. For compression tests, the sample was placed on a Teflon tape covered compression platen; no adhesive was used.

Quasi-static shear tests were conducted to deformations of K = 1 at a rate of approximately 0.02 s⁻¹. For high-rate shear tests, the sled was secured in place with masking tape before loosening the set screw. A steel ball was then dropped from a height of 25 or 81 cm within the drop tube to shear the tissue to failure at a target strain rate of 150 or 300 s⁻¹, respectively.

For compression tests, the sled was secured in place with a set screw so that the compression platen was either 1, 10, or 15 mm above the tissue for QS, 150 s⁻¹, or 300 s⁻¹ tests, respectively. In QS tests, samples were compressed to deformations of 40% at a strain rate of approximately 0.2 s⁻¹. For the 150 and 300 s⁻¹ tests, the steel ball was dropped from a height of 20 cm, with the difference in velocity controlled by securing the moving platen 10 or 15 mm, respectively, above the fixed platen. An aluminum stop was positioned to arrest high-rate compression 3 mm from the bottom of the tissue, producing maximum compressive deformations of 60% for the brainstem and 43% for the cerebellum and cerebrum.

Noise in the load and displacement signals was smoothed using the Butterworth, 4-pole, phaseless filter specified in the SAE J211 standard [50] with a cutoff frequency of 81 Hz. Load and displacement data were then converted into First Piola–Kirchoff (first PK) stress. The first PK stress was calculated by dividing the load signal by the cross-sectional area of the tissue, defined as the product of the measured length and width. Deformation was quantified using either the F₁₁ (⁠$λ$⁠) or F₁₂ (K) component of the deformation gradient for compression or shear tests, respectively. F₁₁ was computed by dividing the measured compressive displacement by the specimen height, while F₁₂ was computed by dividing the measured shear displacement by the specimen height.

For each group, an average stress versus deformation curve was created. This was done by interpolating stress values from each test at deformation gradient increments of 0.01 between 0 and 0.35 for shear and 1.0 and 0.7 for compression (0 to 30%), using the interpl function in MATLAB (MathWorks, Natick, MA). The upper value for shear was chosen due to tissue failure occurring at about K = 0.35 in many high-rate shear tests. For compression, the upper threshold of 0.30 was chosen based on the compression stop height. Average stress values and the corresponding standard deviations at each interpolated point were then determined for each group by averaging the interpolated stress data from individual tests within each group.

Constitutive models were fit to our experimental data in multiple steps. First, we fit three different hyperelastic models, that have previously been shown to provide good quality fits in brain tissue, to our quasi-static data. For each brain region, we performed unidirectional fits, where shear or compression model parameters were optimized using only shear or compression data (two sets of parameters per region), as well as bidirectional fits, where we required the models to fit both the compression and shear data together (one set of parameters per region). We then used the best-performing hyperelastic model as the basis of a hyper-viscoelastic model, which was optimized across all strain rates, again for both unidirectional and bidirectional fits. All constitutive models chosen here are available in commonly used finite element packages.

Hyperelastic models were fit to the QS experimental data using the one-term and two-term Ogden models [51] and the Mooney–Rivlin model [52]. The one- and two-term Ogden models were chosen due to their previously demonstrated ability to fit brain tissue well [11,19,45], while the Mooney–Rivlin model was chosen to evaluate the response with a simpler constitutive model.

in compression [45], where $λ$ is the stretch value from the $F11$ component of the deformation gradient.

in compression.

Here, $Pij$ is the experimental stress and $Pijψ$ is the predicted stress.

and the fit parameter $τ$ is a characteristic time scale.

Three versions of the hyper-viscoelastic model were used to fit the two unidirectional and one bidirectional datasets by minimizing the objective functions in Eqs. (7) and (8), respectively, in matlab. Hyperelastic parameters were fixed to the values determined during the previous quasi-static fitting, so only viscoelastic parameters were allowed to vary during the optimization. Initial values for the viscoelastic parameters were evaluated over several orders of magnitude via a grid search.

For both compression and shear data, statistical analysis was performed on the one-term Ogden parameters that were fit to each sample, such that n = 8 sets of parameters were present per group. manova tests were performed for each mode of loading using the dependent variables of stiffness (⁠$μ$⁠) and nonlinearity (⁠$α$⁠), with independent variables of strain rate and brain region. Where a statistical difference of p < 0.05 was detected, a Tukey Test was performed to examine pairwise differences in either stiffness or nonlinearity between rates and regions.

The average shear stress–strain data for all brain regions at all rates are presented in Fig. 3 (results for individual samples are shown in Supplemental Figure 1 available in the Supplemental Materials on the ASME Digital Collection.). QS data generally demonstrate a concave-down shape, regardless of tissue region, while the shapes of the high-rate curves were variable. Brainstem behavior was concave-up at both 150 and 300 s⁻¹, where the 300 s⁻¹ average demonstrated a larger degree of curvature. The cerebellum showed a concave-down curve for the 150 s⁻¹ group and a sigmoidal-like curve for the 300 s⁻¹ group, while the cerebrum was concave-down for all groups. Stress levels were comparable between all regions. Cerebrum data showed the greatest amount of scatter, followed by cerebellum; variability was highest in high-rate tests. Additionally, high-rate curves showed a stiffer response than QS curves. Note that for the brainstem the two high-rate curves were similar.

Statistical characterization (Table 1) showed that shear strain rate had a significant effect on $μ$ (p < 0.001), with established differences between each of the high rates and the QS rate, but no difference between the two high rates. None of the other variables had a significant effect on $μ$⁠. The only variable that significantly affected $α$ was brain region (p = 0.004), with the comparison between brainstem and cerebrum being significant.

Average compressive stress–strain curves for all regions and rates are presented in Fig. 4 (results for individual samples are included in Supplemental Figure 2 available in the Supplemental Materials). All results demonstrated a concave-up response. The brainstem showed immediate exponential stiffening as soon as deformations were applied, while the cerebrum and cerebellum did not stiffen until about 10% compression. Stress values (and data scatter) were significantly higher for brainstem than for cerebellum and cerebrum, which were similar to one another. It is notable that compressive stress values were an order of magnitude higher than shear stresses. As with the shear data, the high-rate curves were substantially stiffer than the QS curves. For the brainstem, the 150 and 300 s⁻¹ curves were not distinct, producing almost identical responses. For the cerebellum, the 150 s⁻¹ curve was unexpectedly stiffer than the 300 s⁻¹ curve after 12% compression, before which the curves overlap. In contrast, the cerebrum showed a stiffer response in the 300 s⁻¹ curve after 12% compression than the 150 s⁻¹ curve, prior to which the curves overlap, though none of the observed differences between the high-rate curves were statistically significant.

Statistical evaluation revealed that compressive strain rate, as with shear, had a significant effect on $μ$ (p < 0.001), again with established differences between each of the high rates and the QS rate but no difference between the two high rates (Table 2). No other variables had a significant effect on $μ$⁠, and there were no differences in $α$ between compression groups.

The one- and two-term Ogden models resulted in good quality fits for all regions and loading modes (Supplemental Figures 3 and 4 available in the Supplemental Materials on the ASME Digital Collection). The Mooney–Rivlin model had lower quality fits, especially in compression. For all groups, the one- and two-term Ogden models provided equivalent fits in both shear and compression. The Mooney–Rivlin model performed poorly for all groups, resulting in a linear fit in all cases. The one-term Ogden and Mooney–Rivlin model parameters for each fit are shown in Supplemental Table 1 available in the Supplemental Materials, while the two-term Ogden parameters are shown in Supplemental Table 2 available in the Supplemental Materials.

None of the models fit the combined data (bidirectional fits) as well as they did unidirectional data. For all regions except the cerebrum, the two-term Ogden model performed best (Figs. 5 and 6 and Supplemental Tables 3 and 4 available in the Supplemental Materials). In the cerebrum, the two-term model performed best in shear, but the one-term model performed best in compression. The one-term Ogden model performed almost identically to the two-term model for all regions but the cerebrum. The Mooney–Rivlin model once again performed the poorest, resulting in linear predictions of both shear and compression.

As the one-term Ogden model provided almost identical fits to the slightly better performing two-term Ogden model and has 2 less parameters, this model was coupled with a Prony series to simulate the rate dependent (viscoelastic) response of the brain tissue. Unidirectional viscoelastic model fits to data across all strain rates for shear and compression, for each brain region are shown in Figs. 7 and 8, respectively. Model parameters for each of these fits are shown in Table 3. Note that six separate sets of model parameters were derived, one for each combination of brain region and loading mode. Generally, the model was able to capture the response of the experimental data. In particular, it was able to simulate the experimentally observed rate-stiffening behavior for both shear and compression, but the shapes of the predicted curves were not always consistent with the experimental data, at least in shear, likely because these shapes were variable in the experimental data. While model fits showed deviation from experimental values in several cases (most notably in QS shear in the brainstem and cerebellum), the curves remained within one standard deviation of the mean for all cases except for 300 s⁻¹ shear in the brainstem. R² values for the viscoelastic fits in shear ranged from 0.74 to 0.99. Compression models generally fit better than shear models, with R² values between 0.85 to 0.99. For compression, the lowest R² values for each region were the QS rate for the brainstem, 150 s⁻¹ for the cerebellum, and 300 s⁻¹ for the cerebrum, but all compression curves followed the trends of the experimental data.

Bidirectional viscoelastic model fits are shown in Figs. 9 (shear) and 10 (compression). Corresponding model parameters for the combined shear and compression fits are shown in Table 4. In this case, three sets of model parameters were determined, one for each brain region. The bidirectional fits produced notably worse models than the unidirectional fits. This was especially apparent for brainstem and cerebrum, where the QS shear stress was substantially underpredicted, and for cerebellum, where the QS compressive stress was substantially overpredicted. For the high-rate groups, the model did a better job predicting the stress response in both shear and compression. Still, overprediction can be seen in the cerebellum and cerebrum in shear at 150 s⁻¹, and notable underprediction occurred after compressive stretches of 20% in the cerebellum and cerebrum at high rates. R² values for the bidirectional fits in shear ranged from −1.99 to 0.99. The lowest R² values in shear for each region were the QS rates for the brainstem and cerebrum and the 150 s⁻¹ rate for the cerebellum. In contrast to the unidirectional fits, compression models performed similarly to shear models with R² values ranging from −2.75 to 0.98. The lowest R² in each compression group corresponded to the QS rate for the brainstem and the cerebellum, and the 300 s⁻¹ rate for the cerebrum.

The objective of this work was to quantify and model the material behavior of Göttingen minipig brain tissue from the cerebrum, cerebellum, and brainstem in low- and high-rate shear and compression. Results reveal rate stiffening in all regions for both loading modes. The presented viscoelastic model, a combination of a one-term Ogden model and a one-term Prony series, provides a reasonably good fit of the stiffening behavior for shear and compression separately, but it is less accurate when modeling both loading modes together. This is the first study to quantify the response of minipig brain tissue from multiple regions at strain rates relevant to blast injury using a unified set of experimental methods.

It is notable that the observed rate stiffening (i.e., dependence of μ on strain rate) across all brain regions and both loading modes was significant in the comparison between both high rates and the QS rate but not between the high rate (150 and 300 s⁻¹) groups. Generally, the 300 s⁻¹ group had a slightly larger peak stress than the 150 s⁻¹ specimens, but the comparison was not statistically significant. This observation suggests that rate dependent effects may be lost in observational error between the two high-rate groups, as the 150 s⁻¹ rate is several orders of magnitude higher than the 0.02 s⁻¹ rate used for QS tests, while being only half the 300 s⁻¹ rate. It is also possible that a threshold exists above which the tissue no longer stiffens at increasing strain rates, as suggested by Fung's theory of quasi-linear viscoelasticity [55]. Multiple studies have shown that brain tissue subject to oscillatory shear at low strain levels exhibits a nonlinear increase in dynamic modulus, for which parameters fit at only two strain rates likely would not be able to fully capture [56]. Future work should aim to quantify the response of brain tissue at lower intermediate strain rates to more comprehensively characterize the evolution of rate dependence [11,57].

Shear data displayed notable variability in concavity from group to group, both among regions and rates. This phenomenon was especially apparent in the QS shear experiments. While most tests resulted in a concave-down response, a few in each group demonstrated a concave-up behavior. This concavity variation is also present in experimental shear data in brain tissue from multiple species published by Rashid et al. [11] (porcine), Mihai et al. [58] (human), and Haslach et al. ([35,59]) (murine). Notably, Haslach examined not just simple shear but also shear tests where the tissue was subjected to varying degrees of precompression. These experiments showed that tissue subjected to no transverse compression exhibited an initial concave-down response followed by a concave-up response, while samples subjected to transverse compression exhibited a purely concave-down response throughout the test. These results suggest that precompression may have been present for some of our samples. While we attempted to minimize this, a small amount of initial compression was required to get the tissue to adhere to both plates. After lowering the upper plate to generate this force, the plate was again raised to return the tissue to its original height, but our measurement of this was approximate, so it is possible that there was some variability in the degree of tissue compression. Preliminary finite element modeling of quasi-static shear tests subjected to varying degrees of precompression showed increasing softening and a more concave-down stress–strain response as the degree of precompression increased. Considering this, future work should examine this effect, both through additional experiments examining varying degrees of precompression as well as finite element modeling. This is likely to be especially valuable in understanding the response of brain tissue as the tissue seems to expand slightly after it is cut, suggesting that it is under some degree of precompression in vivo.

For both the unidirectional and bidirectional QS datasets, the best-performing hyperelastic model was the two-term Ogden model, though it was only slightly better at fitting data than the one-term Ogden model. Notably the two-term model performed better or equivalent to the one-term model for all regions but the bidirectional fits for the cerebrum, where the one-term performed best in compression. This appears to largely be a result of the better fit of the two-term model in shear, though fit quality between the one- and two-term models remains close. While all models could not always capture the concavity of the shear tests, the Mooney–Rivlin model was especially unable to do so, as in simple shear it reduced to a linear function. The one-term Ogden models presented by Rashid et al. [11] also showed an inability to fit the initial concavity of shear tests, producing a mostly linear fit. In contrast, the one-term Ogden models presented by Budday et al. [19] showed a greater ability to predict the concavity in shear for unidirectional fitting and in both shear and compression for bidirectional fitting. However, constitutive model fits in that study were only fit to shear strains of K = 0.2 and up to 10% compression, compared to the shear strains of K = 0.0–0.35 and compressive strains of up to 30% reported here. As the experimental shear data becomes roughly linear at strains above 0.2, the linear region may be influencing fit quality more so than in other studies with lower strains.

Few studies have attempted to model both shear and compressive loading of brain tissue with a single set of parameters. However, both Budday et al. [19,60] and Moran et al. [61] fit tri-directional experimental data from isolated white and gray matter samples of human brain tissue using shear, tension, and compression data. Interestingly, both groups presented models with a greater fit quality in both shear and compression than shown here. It is possible that this is due to increased variability resulting from the use of mixed white and gray matter samples in our study, as opposed to the isolated gray or white matter samples used by these researchers. Additionally, the inclusion of tensile data in both studies may have contributed to greater fit quality.

In this work, the optimized one-term Ogden model, based on QS shear of cerebrum, resulted in a shear modulus (μ) of 540 Pa and nonlinearity parameter (α) of 1.43 × 10⁻⁶. Several other authors have reported Ogden model fits for porcine cerebral tissue loaded in simple shear at a variety of strain rates. Rashid et al. [11] reported moduli of 1038 Pa at a rate of 30 s⁻¹ and 2073 Pa at a rate of 120 s⁻¹ (corresponding α values of 2.77 and 3.23). Prange and Margulies [22] reported stiffness values of 296 Pa in thalamic gray matter and 254.2 Pa in corona radiata white matter at a strain rate of 4 s⁻¹ (corresponding α values of 0.04 and 0.06), while Coats and Margulies [20] reported stiffnesses of 208 Pa in thalamic gray matter and 190 Pa in white matter from the corona radiata at a rate of 3 s⁻¹ (α values of 0.04 and 0.06). As expected, the modulus of our tissue is substantially lower than that of the much higher rate tests performed by Rashid but relatively close to the lower rate tests performed by Coats and Prange. While these tests have lower stiffness values, it is worth noting that they have higher values for α, which may lead to a lower optimized value of μ.

In compression, our reported QS moduli were 5680 Pa in the brainstem, 2890 Pa in the cerebellum, and 2150 Pa in the cerebrum, with corresponding nonlinearity parameters of 14.3, 12.5, and 11.6. Several authors report comparable values for porcine brain tissue loaded at multiple strain rates. Rashid et al. [45] reported moduli of 6000 Pa at 30 s⁻¹ and 12642 at 90 s⁻¹ (corresponding α values of 0.59 and 5.05), while Singh et al. [62] reported a value of 36 Pa (corresponding α of 9.97) for quasi-static compression. Additionally, Li et al. [63] report shear modulus values, for quasi-static compression of the brainstem, cerebellum, and cerebrum of 8-week old piglets, of 300 Pa, 390 Pa, and 360 Pa, respectively (α values of 7.30, 6.95, and 7.70). Interestingly, our quasi-static results seem most like the 30 s⁻¹ results reported by Rashid, with cerebral samples being roughly 1/3rd the stiffness compared to at least an order of magnitude stiffer than the more similar tests performed by Li and Singh. However, it is worth noting that as Li's tests were in piglets, the brain tissue may be substantially softer than in older pigs [56]. Reported material properties for human brains in quasi-static compression [19] show a stiffness of 1630 Pa and a nonlinearity parameter of 16.5, which is similar to the porcine compression properties we report here.

The unidirectional viscoelastic fits were generally able to capture the response of brain tissue at all rates for all regions, especially at higher rates relevant to head trauma. However, these fits still could not capture the response of the tissue in shear, with the poorest performance at QS rates. The shear fits were more linear than the experimental data, especially for the concave-up experimental curves of the brainstem at high rates. Published viscoelastic models for high-rate simple shear of brain tissue over large deformations are limited ([11,64]). These previous studies used stress relaxation data to fit viscoelastic parameters and reported higher levels of fit quality (R² > 0.89 compared to the R² > 0.74 seen here). These studies also used higher-order Prony series having more than the two viscous parameters used in this work. Many studies model brain tissue using at least a two-term Prony series [65] in order to capture both the short- and long-time viscous response [34]. However, due to the large difference in strain rates examined here, we opted to use a single term and to allow the QS response to be governed by the $G∞$ term in the model, thus reducing the number of parameters needed to be fit and reduce the risk of overfitting.

The viscoelastic models were unable to adequately predict the response of tissue in both compression and shear (bidirectional data), suggesting that additional modeling considerations are needed. Previous work suggests that a biphasic model with a hyper-viscoelastic solid component may be useful [26]. Unfortunately, our experiments were not designed to determine the tissue permeability critical to the biphasic framework [66]. To our knowledge, no models capturing brain response for multiple loading modes at high strain rates have been published to date. However, Wu et al. [67] modeled brain tissue data from shear, compression, and tension tests at strain rates of 0.5 to 30 s⁻¹ using a viscoelastic model. They showed relatively good quality fitting in all loading modes but 30 s⁻¹ tension, though results from their compression model showed a notable difference in shape compared to experimental data, and their model in shear tended to deviate from experimental results at higher strains. In contrast, our model is generally able to capture the response of brain tissue in both shear and compression over a larger range of strain rates. However, a direct comparison between these models is difficult due to the difference in tested strain rates as well as the inclusion of tensile test data in the Wu model, which is not present in ours. Budday et al. ([60,68]) modeled both shear and compression using a viscoelastic framework based on cyclic shear, compression, and tension experiments, but only at low strain rates. Their models showed substantial underestimation of initial cycles of both shear and compression, similar to what was seen with the models for the QS rates in the cerebrum and brainstem in our work. Interestingly their models showed better predictive accuracy for subsequent conditioned cycles. This suggests that the modeling of viscoelasticity alone may be unable to capture the unconditioned response of brain tissue. It may also be that biaxial data, where the tissue is loaded at various levels of compression and shear simultaneously, would be more useful for developing models that capture both responses.

Other considerations may also improve model performance. For example, approaches which capture anisotropy, possibly by explicitly modeling fiber contributions, may improve the fit quality for shear samples, especially in the brainstem and cerebrum where QS stresses were substantially under-predicted by our isotropic viscoelastic model. Future work should aim to model brain tissue as a biphasic material with viscoelastic reinforcing fibers to better model the physical response of brain tissue in addition to relying on the purely phenomenological viscoelastic and hyperelastic models used here.

Our work had several limitations that warrant discussion. First, the experimental stress–strain response includes a significant amount of scatter. This is especially notable in the shear tests of tissue from the cerebellum and cerebrum, though the QS brainstem results also show a large degree of scatter. Li et al. [23] reported standard deviations similar to ours for the brainstem in shear, but our standard deviations for the cerebrum and cerebellum in shear are more than two times higher than theirs for the same tissues. Interestingly, tensile tests from their study show higher spread for the cerebellum than the brainstem or the cerebrum. It is possible that differences in experimental configuration and the exact location of sample harvest from each of the three regions may contribute to this difference.

It is unlikely that the spread seen in our stress and strain data can be explained by outliers that do not necessarily capture the average tissue response. Prior to performing any constitutive model fitting, we performed outlier analysis on the μ and α values of the individual one-term Ogden model fits as well as on the slope of a linear fit to the data from each sample. We plotted all values for each of these parameters in a group together in a boxplot and any sample was excluded if either both the Ogden parameters or the slope of the linear fit was 1.5 times the interquartile range below the 25th percentile value or above the 75th percentile value was excluded if either the stiffness. Test video from all excluded samples was inspected and notable issues with tests such as poor glue adherence or presence of meninges on the outside of the sample was noted in all excluded tests.

To some degree, this level of variation in the cerebellum and cerebrum is expected because the tissue samples were made up of differing proportions of white and gray matter which were not consistent between samples. Multiple studies have demonstrated a difference in the mechanical properties of white and gray matter [19,37,44,69,70], so different proportions would be expected to introduce variability into the results. The effect of this variation may be especially pronounced in shear compared to compression, because compressive response is expected to be strongly influenced by fluid flux, while shear behavior is likely a direct result of the tissue matrix itself. Additionally, variability in cerebellar and cerebral shear could result in slip planes created by the presence of sulci in the samples. This effect would be especially pronounced in the cerebral samples due to the larger number of sulci in the cerebrum. Future work could aim to quantify and model the response of isolated white and gray matter specimens at high rates either using smaller samples to measure each region independently (though these may be subject to boundary conditions), or attempting to capture the contributions of each through inverse modeling.

As the brainstem is comprised almost entirely of white matter, differing contributions of white and gray matter do not explain the high degree of scatter seen in brainstem tests. Instead, potential sources of error may be due to the influence of anisotropy or precompression. While we attempted to keep the orientation of samples the same, we did not apply any markings to the surface of the tissue sample prior to cutting, so it is possible that some samples may have been tested in a slightly different orientation than others. As the brainstem has been shown to be transversely isotropic [33], this could be partially responsible for the spread seen in the brainstem data, especially in shear. Additionally, as previously mentioned, unintentional precompression of the tissue may lead to a softer response and a change in the shape of the stress–strain curve, further compounding error. High-rate tests appear to have a larger degree of variability than QS tests. In shear, all samples exhibit a higher standard deviation, as a percentage of the mean, at 300 s⁻¹ than QS tests, though this is not consistently the case between the QS and 150 s⁻¹ groups. This suggests that some source of experimental error may be amplified at 300 s⁻¹, potentially due to an increased level of vibration in the tester.

An additional source of error in this work may be due to variations in sample geometry. Due to the ultrasoft nature of brain tissue and the complex geometry of the brain, cutting uniformly sized and shaped samples proved difficult, and brain tissue samples would often slightly deform in an unpredictable way once removed from the tissue punch, though we attempted to minimize swelling by using the PEG buffer. While sample measurements were taken after the tissue was removed from the punch, and either glued to one of the two shear platens or placed on the compression platen, some degree of creep may have also occurred prior to testing. Tissue dimension measurements and stress calculations also relied on the assumption that samples were a perfect cuboid, which is likely not the case. To deal with both inhomogeneous deformations and inconsistent sample geometry, future experiments should be designed to capture fiducial markers on multiple planes of the specimen and used to generate a deformation map using digital image correlation (DIC) before fitting a constitutive model using inverse finite element methods. DIC coupled with inverse finite element modeling has been demonstrated as an effective way to fit anisotropic hyper-viscoelastic models of the brainstem in low-rate compression [71] as well as hyperelastic models of soft tissue phantoms [72] and lung tissue [73].

Test videos suggested that induced deformations were not perfectly homogeneous, due to both boundary conditions and inertial effects at high rates. Evaluation of shear test video showed some degree of deformation inhomogeneity at lower strain levels (⁠$K≤0.05$⁠), with areas closer to the moving plate appearing to begin moving several frames before the rest of the tissue. This phenomenon was evident in both QS and high-rate tests for all tested regions. As the level of applied strain increased, the tissue generally began to respond more uniformly. In compression, sample deformation was uniform in QS tests, consistent with a lack of boundary effects on the Teflon-coated platens. At high rates, some samples expanded laterally more at the top (moving) platen than bottom (fixed) platen. No obvious wave propagation was observed in the high-rate shear test videos, but preliminary finite element analysis showed clear stress wave propagation in high-rate simulations in both shear and compression, suggesting that high-rate deformations were not homogenous as assumed. The DIC and inverse modeling approach suggested in the previous paragraph may help clarify this in future research.

Finally, this work did not address anisotropy of the brain tissue. Brain stiffness has been shown to correlate with myelin content [48], which may lead to axons acting as reinforcing fibers and associated anisotropy in aligned regions of the brain. Other studies have previously shown that the brainstem exhibits transversely isotropic behavior in shear [33], though anisotropy in other regions of the brain remains poorly defined, with some authors finding cerebral white matter to exhibit direction dependence [21,39,64,74] and others reporting no direction dependence [19,75,76]. Due to the physiological importance of the brainstem as well as the susceptibility of these regions to diffuse axonal injury in blast [77], future work should aim to better quantify the response of brain tissue in multiple loading directions.

The work presented here comprises the first set of experimentally derived constitutive models for minipig cerebrum, cerebellum, and brainstem in response to simple shear and unconfined compression at rates consistent with blast exposure. Our results show that the response of brain tissue can be well modeled in either compression or shear over a range of strain rates using a hyper-viscoelastic model with a one-term Ogden hyperelastic component. Bidirectional models, which fit both shear and compression data to a single set of parameters, generally perform worse, suggesting additional parameters may be needed to capture the compression and shear response together.

The opinions and assertions contained herein are the private views of the authors and are not to be construed as official or as reflecting the views of the U.S. Army, the U.S. Department of Defense (DoD), or The Henry M. Jackson Foundation for the Advancement of Military Medicine, Inc. (HJF). This paper has been approved for public release with unlimited distribution.

• National Science Foundation (Award No. 2027367; Funder IDs: 10.13039/100000001 and 10.13039/100000147).

• U.S. Army Medical Research and Materiel Command (Award No. W81XWH-17-2-0008; Funder IDs: 10.13039/100000182 and 10.13039/100016156).

• U.S. DoD (Funder ID: 10.13039/100000005).

• U.S. Army Military Operational Medicine Research Program Area Directorate (Funder ID: 10.13039/100000090).

$C$ =

right Cauchy-Green tensor

$C1$ and $C2$ =

coefficients of the Mooney–Rivlin hyperelastic model (Pa)

$F$ =

deformation gradient

$G∞$ =

long term time response coefficient of the viscoelastic model

$G1$ =

viscoelastic coefficient scaling response at relaxation time $τ$

$K$ =

amount of shear defined as the component of $F$ in the 1,2 (shear) direction

$I1$ and $I2$ =

first and second invariants of $C$

$P$ =

first Piola–Kirchoff (first PK) stress (Pa)

$W$ =

strain energy

$μ$ =

Ogden hyperelastic model stiffness coefficient (Pa)

$λ$ =

component of $F$ in the 1,1 (compressive) direction

$α$ =

Ogden hyperelastic model non-linearity coefficient

$λi$ =

ith principle strain

$χ$ =

objective function value

$t$ =

time (s)

$τ$ =

relaxation time

Subscripts and Superscripts

bi =

bidirectional

e =

elastic

max =

maximum value

MR =

Mooney–Rivlin hyperelastic model

Ogd =

n-term Ogden hyperelastic model

uni =

unidirectional

$ψ$ =

predicted value

Acronyms and Abbreviations

BS =

brainstem

CB =

cerebellum

C =

cerebrum

First PK =

first Piola-Kirchoff stress

MR =

Mooney–Rivlin hyperelastic model

PEG =

polyethylene glycol

QS =

quasi-static

TBI =

traumatic brain injury

1T Ogden =

one-term Ogden hyperelastic model

2T Ogden =

two-term Ogden hyperelastic model