Valentina
Balbi
*^{a},
Antonia
Trotta
^{b},
Michel
Destrade
^{ab} and
Aisling
Ní Annaidh
^{b}
^{a}School of Mathematics, Statistics and Applied Mathematics, NUI Galway, University Road, Galway, Ireland. E-mail: vbalbi@nuigalway.ie; v.balbiv@gmail.com
^{b}School of Mechanical & Materials Engineering, University College Dublin, Belfield, Dublin 4, Ireland
First published on 13th June 2019
We investigate experimentally and model theoretically the mechanical behaviour of brain matter in torsion. Using a strain-controlled rheometer, we perform torsion tests on fresh porcine brain samples. We quantify the torque and the normal force required to twist a cylindrical sample at constant twist rate. Data fitting gives a mean value for the shear modulus of μ = 900 ± 312 Pa and for the second Mooney–Rivlin parameter of c_{2} = 297 ± 189 Pa, indicative of extreme softness. Our results show that brain always displays a positive Poynting effect; in other words, it expands in the direction perpendicular to the plane of twisting. We validate the experiments with finite element simulations and show that when a human head experiences a twisting motion in the horizontal plane, the brain can experience large forces in the axial direction.
Simple shear and torsion tests are particularly useful to study the Poynting effect, a typical nonlinear phenomenon displayed by soft solids. When sheared or twisted, those materials tend to elongate (positive Poynting effect^{7}) or contract (negative Poynting effect^{8}) in the direction perpendicular to the shearing or twisting plane. This phenomenon has been observed for brain matter in simple shear tests.^{6} However, a practical limitation of simple shear tests is that to date there are no shearing devices able to measure and quantify the normal force, which limits the determination of material parameters.
An alternative test is torsion. It can be performed by glueing a cylindrical sample between two parallel plates and then applying a twist to the sample by rotating one plate with respect to the other. In the past, torsion tests on brain matter have been performed using a rheometer, at constant strain rates (from 0.05 to 1 s^{−1}),^{9} to measure the elastic properties of the tissue and to investigate its viscoelastic behaviour dynamically over a range of frequencies (20–200 Hz).^{10} A comprehensive summary of the results of mechanical tests on brain tissue can be found in two recent reviews.^{11,12} However, in these studies, only the torques were recorded. Moreover, torsion was modelled as simple shear, an equivalence that, as we show in the ESI,† is only valid locally. To our knowledge, the role played by the normal forces arising during torsion of brain matter has not been investigated yet.
In this work, we perform torsion tests on cylindrical porcine brain samples and measure the torque and the axial force required to twist the samples at a constant twist rate of 300 rad m^{−1} s^{−1}. In Section 2, we describe the experimental protocols for preparing and testing the brain specimens. In Section 3, we present the collected and filtered data and describe the filtering strategy adopted to keep meaningful experimental measurements. We then accurately model the data with the Mooney–Rivlin model in Section 4 and obtain material parameters of brain matter that compare well with those found from other tests. To further validate the analytical modelling, we implement finite element (FE) simulations in Abaqus to mimic the experiments and finally, we use the estimated mechanical parameters to simulate a rotational head impact.
Our main finding is that brain matter exhibits the normal Poynting effect, i.e. it tends to expand along its axis when twisted.^{13} As noted by Rivlin,^{7} the Poynting effect is a nonlinear elastic effect par excellence and cannot be explained by the linearised theory. It was present in all the cylindrical samples we tested, and when we simulated the rapid twisting of a head in a finite element model, we found that large vertical stresses developed in the whole brain also.
The scalp was removed using a scalpel and the cranial bone was removed using an oscillating saw. Following removal of the skull, the meninges tissue was removed using surgical scissors. Finally, following resection of connective and vascular tissue and separation from the spinal cord, the undamaged brain was placed in phosphate buffered saline (PBS) solution.
The brains had an average (maximum) length, width and height of 7.5 ± 0.3 cm, 6.4 ± 0.3 cm and 3 ± 0.3 cm, respectively. From each brain, three slices of approximately 15 mm were obtained from the frontal, parietal and occipital portions, respectively. Then, a stainless steel cylindrical punch of 20 mm diameter was used to remove up to four cylinders from each slice (depending on the available surface of the slice) as shown in Fig. 2. Each long cylindrical sample was then cut to samples of approximately 10 mm using a scalpel and template. The exact height of the specimen was measured again prior to testing. A total of 9 mixed grey and white matter samples were then used for testing and modelling in this study. After extraction, the cylindrical samples were placed in PBS solution in multi-well plates of 20 mm diameter and placed in a fridge for less than 2 hours while all samples were prepared.
Fig. 3 shows a set of typical output data from the rheometer. In Fig. 3(a), the torque τ and the normal force N_{z} are plotted against the twist ϕ and the twist squared, respectively. The output twist ϕ in the plots is the angle of rotation α per unit length ϕ = α/H. From the data shown in Fig. 3(a), we identify three regions, for both the torque and the force data: (i) a noisy region (in black), at the very beginning of the test; (ii) a linear region (in purple) bounded by a maximum (minimum), and (iii) a decaying region (in orange) towards the end of the experiment.
The initial noisy region is due to the delay time of the instrument in reaching a constant twist rate. The plot in Fig. 3(b) shows the twist rate against the twist and clearly highlights the initial region where the upper plate is accelerating to reach the constant twist rate = 300 rad m^{−1} s^{−1}. The plots in Fig. 3(c) further show that some data were actually generated at twist rates lower than 300 rad m^{−1} s^{−1}. Only the data generated at 300 rad m^{−1} s^{−1} were thus considered in the following analysis.
The other filtering criterion is the breaking point of the sample, which is identified clearly by a steep drop and rise in the plots of Fig. 3(a), indicating that an irreversible change in the mechanical response of the tissue occurred. Therefore, all data points after the breaking point were discarded.
The remaining “good” data obtained from nine samples S_{1},…,S_{9} are shown in Fig. 4. Note that the rheometer reported positive values for the normal force, showing that the samples tend to expand axially during the twist (positive Poynting effect). Here, we changed the sign of the data to be consistent with the modelling approach presented in the next section.
As mentioned in Section 2.2, the normal force was set to zero before commencing each test. However, the force transducer of the HR2 rheometer has a sensitivity of 0.01 N, so that variations of the force within that range are not detected by the instrument. We therefore expect that the sample undergoes a small contraction prior to the transducer picking up a meaningful value for the force. Mathematically, we superpose an axial contraction to the actual rotation so that the total deformation is written in cylindrical coordinates as follows:
(1) |
(2) |
We are interested in the elastic behaviour of brain matter, which we assume to be isotropic and incompressible. In view of the linear dependence of the torque with respect to the twist and of the normal force with respect to the twist squared highlighted by the results in Fig. 4, we conclude that the constitutive behaviour of the brain must^{14} be modelled with a Mooney–Rivlin strain energy function:^{6,15}
W = c_{1}(I_{1} − 3) + c_{2}(I_{2} − 3), | (3) |
σ = 2c_{1}B − 2c_{2}B^{−1} − pI, | (4) |
The principal stretches are the square roots of the eigenvalues of B. The intermediate stretch is associated with the radial direction, and the maximum and minimum stretches λ_{2} and λ_{3} are obtained by solving the following equations:
(5) |
The elastic equilibrium of the deformation is translated as the following problem:
(6) |
(7) |
These formulas were first established by Rivlin,^{7} see Appendix A for details.
Now, the torque and the normal force that have to be applied to the cylinder to maintain the deformation in (1) are:
(8) |
(9) |
(10) |
2(λ^{3} − 1) − R_{0}^{2}λ^{2} = 0. | (11) |
Note, as expected for the Mooney–Rivlin model, the linear dependence of the torque on the twist and of the normal force on the twist squared, see (8) and (9).
Because the shear modulus μ = 2(c_{1} + c_{2}) is positive, it follows that and − are positive when λ < 1 (when the cylinder is contracted prior to the twist). The coefficient is associated with the Poynting effect displayed by the sample ( < 0 corresponds to the normal Poynting effect) and is due almost entirely to the twist, whereas the coefficient accounts for the pre-stretch only. When λ = 1, i.e. in pure torsion, we have = πR_{0}^{4}(c_{1} + c_{2}), = −πR_{0}^{4}(c_{1}/2 + c_{2}), and = 0. The values of λ in Table 2 show that the samples S_{2} and S_{3} experience less than 1% pre-stretch; hence, for those two samples, provides an effective measure of the exact Poynting effect, i.e. in the absence of a normal compressive force, the samples would expand axially.
Moreover, the fit on {ϕ^{2},N_{z}} uses a weighted (with respect to ϕ^{2}) least squares method. The weights are given by the values of ϕ^{2}, which increase over time as the test progresses. Therefore, data collected at the beginning of the test (which are more uncertain) are weighted less than those collected at the end of the test. Finally, we input the coefficients , and into eqn (10) to get the elastic parameters c_{1} and c_{2}.
The results of the linear regression are shown in Table 1. The mean values for the elastic parameters are μ = 900 Pa and c_{2} = 297 Pa, respectively. The dimensions of the nine samples after the compression are summarised in Table 2, where the length _{0} and radius r_{0} are given for each sample, as well as the pre-stretch computed from (11) and the maximum values of the highest principal stretch λ^{max}_{2}, attained at breaking point close to the top face.
Sample | μ [Pa] | c _{2} [Pa] | R _{ τ } ^{2} | R _{ N z } ^{2} |
---|---|---|---|---|
S_{1} | 1232.50 | 294.45 | 0.999 | 0.946 |
S_{2} | 1092.31 | 235.84 | 0.998 | 0.939 |
S_{3} | 766.95 | 310.60 | 0.988 | 0.966 |
S_{4} | 491.14 | 201.22 | 0.996 | 0.958 |
S_{5} | 656.96 | 59.68 | 0.988 | 0.87 |
S_{6} | 644.59 | 113.75 | 0.994 | 0.92 |
S_{7} | 952.36 | 347.26 | 0.993 | 0.947 |
S_{8} | 803.12 | 401.41 | 0.997 | 0.925 |
S_{9} | 1460.17 | 710.29 | 0.995 | 0.962 |
Mean ± SD | 900 ± 312 | 297 ± 189 |
Sample | λ | _{0} [mm] | r _{0} [mm] | λ ^{max}_{2} |
---|---|---|---|---|
S_{1} | 0.93 | 12.62 | 10.36 | 1.82 |
S_{2} | 0.99 | 16.02 | 10.05 | 1.73 |
S_{3} | 0.99 | 15.92 | 10.05 | 1.81 |
S_{4} | 0.89 | 12.95 | 10.59 | 2.14 |
S_{5} | 0.94 | 13.26 | 10.31 | 1.85 |
S_{6} | 0.85 | 14.19 | 10.84 | 2.04 |
S_{7} | 0.89 | 9.51 | 10.59 | 1.87 |
S_{8} | 0.89 | 12.22 | 10.59 | 1.68 |
S_{9} | 0.95 | 10.14 | 10.26 | 1.69 |
The initial cylindrical geometry of the sample was obtained by setting the radius R_{0} = 10 mm and the height L_{0} = _{0}/λ, calculated according to Table 2. We used a mesh of 78750 hexahedral elements (C3D8) with hybrid formulation to reproduce exact incompressibility and we assigned the Mooney–Rivlin model with material parameters in Table 1 to account for hyperelasticity.
To simulate the twist, we first defined a reference point at the centre of the top surface of the cylinder, which we then coupled with all points on the surface, and finally, we assigned a rotational displacement around the longitudinal axis (ramp form of amplitude 3 rad) whilst setting the other degrees of freedom to zero. The bottom surface of the cylinder was encastred. The output variables were the resultant axial force (RF3) and torque (RM3).
To perform the simulations, we chose two specimens: S_{1} and S_{2}. An additional step, prior to the torsion, was added to simulate the 7% pre-compression (see Table 2) undergone by S_{1}.
The results are shown in Fig. 5. We see that the numerical simulations validate the predictions of the analytical model described in Section 4.1. The torque and the normal force calculated in Abaqus are consistent with the analytical predictions and the measured data for both cases with and without pre-compression (S_{1} and S_{2}, respectively). We note that there is a small mismatch between the analytical and the numerical normal force for S_{1}. This is due to the longitudinal bulging of the sample occurring during the 7% pre-compression phase, which results in a non-homogeneous deformation along the axis of the cylinder.
Fig. 5 Comparison of the resultant torque τ and normal force N_{z} for S_{1} and S_{2}. Results of the numerical simulations in Abaqus (triangles), analytical predictions with the models (8) and (9) (solid lines) and experimental data (red circles for S_{1} and orange circles for S_{2}). |
In addition, here, we are able to directly estimate the second Mooney–Rivlin coefficient c_{2} from the normal force data and thus provide a direct measure of the Poynting effect. This quantity cannot be measured from shear stress data alone, although its sign (positive for porcine brain matter) can be deduced by piercing a hole in one of the platens.^{6} We note also that in contrast to simple shear, where a neo-Hookean material (c_{2} = 0) does not display the Poynting effect,^{6} the same material does have a non-zero normal force in torsion.^{21}
Finally, it should be noted that the influences of tissue orientation and anisotropy have not been considered in this study, as the 9 tested samples were all excised from slices in the coronal plane. The orientation of axons within the white matter might play a role in determining the constitutive behaviour of the brain at different length-scales.^{22,23} Testing of brain samples excised from different planes could therefore provide useful information on the effect of anisotropy both at the micro and macro scales.
Another limitation of the study is that all tests were performed at room temperature. Although a general consensus on the matter has not been reached yet, discussions on the effect of temperature on brain stiffness can be found in ref. 24–26.
Our conclusion is that brain matter exhibits a real and large Poynting effect in torsion, which is bound to lead to the development of large normal forces in an impacted brain.
In order to investigate the existence and magnitude of axial forces during twisting head impacts, we used Abaqus/Explicit to simulate a rotational impact with the University College Dublin Brain Trauma Model (UCDBTM) developed by Horgan and Gilchrist.^{27}
We applied a rotational acceleration in the axial plane to the centre of gravity of the head, peaking at 2170 rad s^{−2}. This value of rotational acceleration is in the range of accelerations experienced in boxing.^{28} We used the mean values of the estimated parameters μ and c_{2} in Table 1 to model the mechanical behaviour of the brain matter. In Fig. 6, the distribution of the axial stress component σ_{33} (where σ is the Cauchy stress) throughout the brain is shown for a rotational angle of α = 0.52 rad. The three sections in the sagittal, coronal and axial planes, respectively, highlight areas where the stress reaches peaks of magnitude in the thousands of Pa. These values indicate that vertical stresses developing during a rotational head impact are of the same order of magnitude as the stiffness μ and are five to ten times larger than the shear stress component σ_{23}, even when the hydrostatic pressure is removed, see ESI,† for more details.
Therefore, the high normal stresses developing during rotational impacts could potentially contribute to traumatic brain injury (TBI) during this type of impact. The results presented in this work thus open the path towards further studies to quantify the role played by normal forces in TBI, in particular with a view to defining more accurate threshold criteria for TBI.
(12) |
(13) |
Then, by substituting (2) into (13) and then into (12) we obtain (7).
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9sm00131j |
This journal is © The Royal Society of Chemistry 2019 |