Poynting eﬀect of brain matter in torsion †

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 m = 900 (cid:2) 312 Pa and for the second Mooney–Rivlin parameter of c 2 = 297 (cid:2) 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.


Introduction
The brain is an extremely soft and fragile tissue, making it hard to test experimentally using the standard protocols in place for soft matter such as elastomers or rubbers. 1 The inflation test is not appropriate for obvious reasons. The uni-axial tensile test requires a dogbone geometry and clamping, which cannot be realised in practice for brain tissue; instead, the end faces of a cylinder have to be glued to the tension plates, which leads rapidly to an inhomogeneous deformation. 2,3 The compression test can achieve homogeneous deformation with lubrication of the plates, but only up to about 10% strain, after which it starts to bulge out. 4 By contrast, the simple shear test works well 5,6 up to a 451 tilting angle leading to a maximal stretch of more than 60%. Similarly, as we show in this paper, the torsion test can be implemented readily for brain matter.
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 . 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.

Materials and methods
In this section, we give a brief description of the procedure for preparation and testing of the brain samples.

Tissue preparation
Two fresh porcine heads from freshly killed 22-week-old mixed sex pigs were obtained from a BRC-regulated (British Retail Consortium) food processing facility (Dawn Pork & Bacon, Waterford, Ireland). The animals did not die for this study; therefore, we did not need to require ethical approval from the Ethics Committee at University College Dublin.
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 AE 0.3 cm, 6.4 AE 0.3 cm and 3 AE 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.

Mechanical testing
A Discovery HR2 Hybrid rheometer with parallel plates was used for all mechanical testing. This device has a torque resolution of 0.1 nN m and a normal force resolution of 0.5 mN. To enable easy removal and to protect the platens, masking tape was applied to both faces of the parallel plates, 3,4,6 and then the samples were glued rigidly to the tape using a high viscosity glue. The tape thickness (r1 mm) is negligible compared to the height of the samples. All testing was performed at room temperature (231-251) and the samples were kept hydrated until the beginning of the test. A cylindrical Peltier plate with radius r p = 10 mm was used. The rheometer was controlled through the TRIOS Software (v4.3.1). Each of the 9 samples was tested once. The test consisted of a single stress growth step: a ramp test, where the upper plate is rotated at a constant twist rate for a duration of 10 s. The distance H between the top plate and the bottom of the instrument was adjusted until the normal force was zero at the beginning of each test. The twist rate (angular velocity of the upper plate per unit height) was _ f = 300 rad m À1 s À1 . The testing protocol was validated on silicon gel samples (see the ESI, † for details).

Experimental results
In this section, we describe the filtering procedure required to get a clean set of data, ready for model fitting and parameter estimation. Fig. 3 shows a set of typical output data from the rheometer. In Fig. 3(a), the torque t and the normal force N z are plotted against the twist f and the twist squared, respectively. The output twist f in the plots is the angle of rotation a per unit length f = a/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.

View Article Online
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 _ f = 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.

Modelling
To fit the experimental data and get a quantitative estimation of the behaviour of the brain in torsion, we first analyse the data,  and then reproduce the mechanical tests theoretically. Finally, we perform finite element simulations in Abaqus.

Theory
Here, we calculate the torque t and the normal force N z required to maintain a cylindrical sample of initial radius R 0 and initial height L 0 in a state of torsion. 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: where l is the (tensile or compressive) pre-stretch, f = a/(lL 0 ) is the twist per unit height and a is the angle of rotation in radians.
Here, (R,Y,Z) and (r,y,z) identify the coordinates of vector points X = RE r + YE Y + ZE Z and x = re r + ye y + ze z in the undeformed (initial) and deformed (twisted) configurations of the sample, respectively. The deformation gradient F = qx/qX associated with the deformation in eqn (1) is then: 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 where c 1 and c 2 are constants, I 1 = tr[B], I 2 = tr[B À1 ] and B = FF T . For this model, the shear modulus is m = 2(c 1 + c 2 ). The corresponding constitutive equation for the Cauchy stress r reads: where p is the Lagrange multiplier introduced to enforce incompressibility and I is the identity matrix. 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 l 2 and l 3 are obtained by solving the following equations: The elastic equilibrium of the deformation is translated as the following problem: d dr s rr ðrÞ þ s rr ðrÞ À s yy ðrÞ r ¼ 0; s rr ðr 0 Þ ¼ 0; with solution: These formulas were first established by Rivlin, 7 see Appendix A for details. Now, the torque t ¼ 2p Ð R 0= ffiffi l p 0 r 2 s yz ðrÞdr and the normal force rs zz ðrÞdr that have to be applied to the cylinder to maintain the deformation in (1) are: where the constants A, B, and C are introduced, the parameters c 1 and c 2 are linked by the following relations: and the pre-stretch l is the unique real and positive root of the following cubic: 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 m = 2(c 1 + c 2 ) is positive, it follows that A and ÀC are positive when l o 1 (when the cylinder is contracted prior to the twist). The coefficient B is associated with the Poynting effect displayed by the sample (B o 0 corresponds to the normal Poynting effect) and is due almost entirely to the twist, whereas the coefficient C accounts for the pre-stretch only. When l = 1, i.e. in pure torsion, we have A = pR 0 4 (c 1 + c 2 ), B = ÀpR 0 4 (c 1 /2 + c 2 ), and C = 0. The values of l in Table 2 show that the samples S 2 and S 3 experience less than 1% pre-stretch; hence, for those two samples, B provides an effective measure of the exact Poynting effect, i.e. in the absence of a normal compressive force, the samples would expand axially.

Parameter estimation
To fit the data in Fig. 4 (9). Moreover, the fit on {f 2 ,N z } uses a weighted (with respect to f 2 ) least squares method. The weights are given by the values of f 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 A, B and C 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 m = 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 c 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 l max 2 , attained at breaking point close to the top face.

Computational validation
We performed brain torsion simulations using ABAQUS Standard 6.14-1 to validate our analytical modelling of the deformation.
The initial cylindrical geometry of the sample was obtained by setting the radius R 0 = 10 mm and the height L 0 = c 0 /l, calculated according to Table 2. We used a mesh of 78 750 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  Table 2 Geometry of the samples after pre-compression, prior to twisting: the estimated axial stretch l, the length c 0 = lL 0 (measured by the instrument), the radius r 0 = l À1/2 R 0 of the nine samples and the maximum value of the greatest principal stretch l 2 before sample breaking  View Article Online 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.

Discussion and conclusions
Here, we note that the results presented in Table 1 show that cylindrical samples of brain matter with initial radius R 0 = 10 mm twisted at a twist rate _ f = 300 rad m À1 s À1 behave as Mooney-Rivlin materials with a shear modulus m = 900 AE 312 Pa. This value is in the same range of values found by Rashid et al. 16 when a block of porcine brain matter was sheared at a shear rate of _ k ¼ 30 s À1 , estimated to be conducive to diffuse axonal injury. 17,18 Moreover, from Table 2, we note that the values of the principal stretch (greatest stretch) l 2 at the breaking point are in the range of 1.67-2.14, which corresponds to extensions of 67% to 114%. These values of strain are well above the estimated axonal strain thresholds associated with diffuse axonal injury (40.05) 19 and to white matter damage in the optical nerve (40.34). 20 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 m 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 s 33 (where r is the Cauchy stress) throughout the brain is shown for a rotational angle of a = 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 m and are five to ten times larger than the shear stress component s 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.

Conflicts of interest
There are no conflicts to declare.

Appendix A
By integrating the first of (6) together with the initial condition, we obtain: Then, by substituting (2) into (13) and then into (12) we obtain (7). The results are plotted in Figure 1: the torque and the normal force data (points) are fitted with a Mooney Rivlin model (solid lines) and the coefficient of determination R 2 is shown for each set of data. The estimated elastic parameters µ and c 2 are reported in Table 1b. The normal force data clearly show that a Positive Poynting effect is measured for silicon cylinders as well. (The sign has been inverted from the original output to be consistent with the modelling convention and with the main paper).  Figure 2: Original collected data on Sample S4: twist rate against time, the ramp data are shown in black; the ramp time is approximately 0.25 sec.
3 Rotational head impact simulations 23 [Pa] [Pa] Here the orange and blue areas correspond to peak stress magnitude of approximatively 4 kPa, same as for the vertical stress component σ 33 , see Figure 6 in the main paper.

A remark on differences and similarities between simple shear and torsion
We briefly compare the results obtained here for torsion tests with those obtained elsewhere for simple shear tests and for torsion modelled as simple shear. We begin by recalling that the deformation gradient for uni-axial compression in the Z direction, followed by simple shear of amount κ in the Y Z plane, has the form ?
Hence, we see from comparison with Equation (2) in the main article that there is a formal connection between torsion and simple shear. However, the equivalence is local only, as simple shear is homogeneous but torsion is not: the amount of "shear" experienced by an element in torsion (κ = rφ = rα/H) depends on the dimensions of the sample and the position of the element. Thus, it does not make sense to compare the amount of shear and the shear rate experienced by all elements in a simple shear experiment with the amount of "shear" and the "shear" rate experienced by a given element at a given location for a given sample dimension in a torsion experiment. Despite this disconnect, finite shear and torsion are often confused in the literature, and torsion experiments in rheometers are routinely modelled as simple shear, see for example the papers cited in the extensive review by Chatelin et al. ? .