Cognition based bTBI mechanistic criteria; a tool for preventive and therapeutic innovations

Blast-induced traumatic brain injury has been associated with neurodegenerative and neuropsychiatric disorders. To date, although damage due to oxidative stress appears to be important, the specific mechanistic causes of such disorders remain elusive. Here, to determine the mechanical variables governing the tissue damage eventually cascading into cognitive deficits, we performed a study on the mechanics of rat brain under blast conditions. To this end, experiments were carried out to analyse and correlate post-injury oxidative stress distribution with cognitive deficits on a live rat exposed to blast. A computational model of the rat head was developed from imaging data and validated against in vivo brain displacement measurements. The blast event was reconstructed in silico to provide mechanistic thresholds that best correlate with cognitive damage at the regional neuronal tissue level, irrespectively of the shape or size of the brain tissue types. This approach was leveraged on a human head model where the prediction of cognitive deficits was shown to correlate with literature findings. The mechanistic insights from this work were finally used to propose a novel protective device design roadmap and potential avenues for therapeutic innovations against blast traumatic brain injury.

duration, and 5 kPa minimum underpressure. Rats were positioned outside the shock tube within one tube diameter to maintain shock overpressure conditions (4), as confirmed previously by shadowgraphy high-speed-camera imaging (unpublished data). The rats' heads were restrained to limit gross head motion via a stereotaxic restraint device. These conditions have been validated to result in a mild severity post-injury phenotype (1,5).

Brain deformation recordings
Traces of point deformations, published previously (2) in the living rat brain during blast exposure were used to validate the computational model. Briefly, an implantable soft magnet transducer was placed on the surface of the rat brain prior to blast exposure. A three-member giant magnetoresistance (GMR) reference sensor array was then affixed to the surface of the rat skull. The GMR array was empirically calibrated with a micromanipulator stage to track the soft magnet transducer with 10 µm spatial resolution and 25 kHz temporal resolution. The calibration results were confirmed and replicated with a 2-equation, 3-variable Gaussian computational model (2). The GMR array wirelessly tracked the relative position of the implanted soft magnet in real-time by detecting dynamic changes in the magnetic field. Relative displacements between the surface of the brain and the skull ranging from hundreds of microns up to a millimeter were detected on the surface of the brain using this method under the mild bTBI conditions described above (2).

Acquisition of anatomical and diffusion-weighted magnetic resonance images of the rat brain
T1-and T2-weighted anatomical magnetic resonance imaging (MRI) data were acquired on a 7 Tesla (T) small animal scanner (Bruker Biosystems) using the Paravision 6.0.1 software platform. The same probe (RF RES 300 1H 112/086 QSN TO AD, Bruker Biosystems), gradient coil (BA-GA12SHP BC 70/30), and surface coil (RF SUC 300 1H R BR QSN RO AD, Bruker Biosystems) were used for all image acquisitions. Rat anaesthesia was initiated with 4% isoflurane in air via an inhalation box and anaesthesia delivery system (SomnoSuite Low-Flow Anesthesia System, Kent Scientific). After adequate anaesthetic depth was achieved, the rat was moved to the scanning platform and secured in plastic, MRI-compatible stereotaxic head fixation apparatus (custom made via 3D printing) with bite bar and ear bars. Isoflurane was continuously administered via nosecone between 1-2% to maintain anaesthesia throughout the duration of the scanning procedure. Respiration rate (30-40 breaths/minute) and body temperature (37 ⁰C) were monitored and maintained via minor anaesthesia and warming surface temperature adjustments throughout the scanning procedure. After positioning the animal in the scanner, wobbling was performed to tune and match the scanner, after which a simple localizer scan was run with B0 adjustment. Upon completion of the localizer, a rapid low-resolution T2 scan was run to visualise the brain in its entirety and incorporate an ellipsoid mapshim procedure. The ellipsoid was aligned with the brain and manually adjusted to encompass only brain tissue while excluding skin, muscle, and skull tissues. T1-and T2-weighted images were then acquired. For both scan types, a field of view (FOV) of 40 mm × 40 mm × 30 mm (x, y, z where x = lateral, y = dorsoventral, and z = rostrocaudal) was used. The FOV was composed of 512 × 512 × 60 voxels, each of size 0.078 mm × 0.078 mm × 0.5 mm, and was acquired in interleaved coronal slices acquired dorsal to ventral with a FOV saturation pulse used on the ventral aspect of the brain to prevent edge artefacts. T2 scans were acquired with an echo time of 12.615 ms, repetition time of 7308.8 ms, RARE factor of 8, 6 averages, and 1 repetition. T1 scans were acquired with an echo time of 10.625 ms, repetition time of 3215.9 ms, 6 averages, and 1 repetition. While the same animal was still anaesthetised, diffusion-weighted images of the brain were captured using a multi-segment gradient echoplanar imaging (EPI) sequence. The acquisition settings were adapted from published datasets (6)(7)(8). For diffusion-weighted image acquisition, the centre of mass remained aligned with the T1-and T1-weighted images, the FOV was reduced to 36 mm × 18 mm × 27 mm (voxel dimensions 0.2 mm × 0.2mm × 0.5 mm) in order to minimise array size and reduce image acquisition time. Scans were acquired with 5 A0 images (b = 0 s/mm 2 ; referred heretofore as "B0 volumes") followed by diffusion-weighted images for 30 gradient directions (b = 650 s/mm 2 ). The entire scan was completed in 6 segments with 2 repetitions. Echo time was 29 ms. Repetition time was 6,500 ms. The gradient duration and dwell time were 3 ms and 10 ms, respectively. Images were visually compared with published datasets and found to be qualitatively similar (6)(7)(8). Major fibre tracts with known anatomical orientation (i.e., genu of corpus callosum, cingulate gyrus) were inspected and confirmed to possess anatomically accurate orientations.

Alignment and processing of diffusion-weighted images
All image processing was performed in FSL 5.0.9 (9). Eddy current correction was applied to the raw diffusion image set using the FDT Diffuson Toolbox (10). Brain masking (for both anatomical and diffusion data) was manually performed in FSLView to isolate the brain from surrounding tissues. The masked B0 volumes were extracted from the rest of the diffusion dataset and averaged. This averaged, masked B0 volume was aligned to the masked T1-weighted anatomical reference image using the FLIRT Toolbox (11). Spline-interpolated rigid body (6 parameters) registration was used. The transformation matrix generated by the prior registration was then saved and applied to the entire diffusion dataset (B0 and gradient direction volumes). The diffusion tensor image (DTI) fibre orientation and fractional anisotropy (FA) data were then obtained using the FDT Diffusion Toolbox and "dtifit" function.

Grey and white matter of brain
The grey matter generally exhibits an isotropic mechanical response (12)(13). However, the white matter bundles of axons induce a transversely isotropic mechanical behaviour. In this work, the mechanical behaviour of the grey matter was modelled as isotropic and, the mechanical behaviour of the white matter was modelled by adding the contribution of an anisotropic response associated to the axonal response to the underlying glial matrix (12,14,15). In addition, many authors have observed a strong dependence of both grey and white matter mechanical behaviours on strain rate (13,16). Following these observations, a viscous contribution that incorporates strain rate dependency was added to the constitutive model.
The strain-energy functions from which the stress expressions can be derived read as: where Ψ iso m and Ψ iso a are the strain-energy functions associated to the glial matrix and to the axons respectively, for the white matter, while Ψ iso a can simply be discarded for the grey matter. Subsequently, only the white matter model is presented with the assumption that the grey matter model can readily be obtained by considering that Ψ iso a = 0 and that the grey matter behaves similarly to the white matter glial matrix as a first approximation.
The isotropic glial matrix contribution is defined by using a Gent strain-energy function commonly used in biological tissues (17) and dissociating the strain rate-independent Ψ iso mri ( * ) from the strain rate-dependent Ψ iso mrd ( * ) contributions: where μ m and μ v are material parameters associated with the shear modulus, j m and j v are dimensionless parameters controlling the limited chain extensibility, I 1 * = tr( * ) is the general deviatoric first strain invariant and I 1 e * = tr( * ) is the viscous deviatoric first strain invariant. The distortional right Cauchy-Green deformation tensors are defined by * = ( * ) T * and * = ( * ) T * , where the distortional part of the deformation gradient is split into elastic and viscous parts as * = * * .
To complete the formulation of the glial matrix response, the definition of the viscous flow is needed and its evolution is determined by: where is the viscous component of the velocity gradient, γ̇v is the viscoelastic multiplier and is the direction of the viscoelastic flow defined as: and where γ̇o v is a dimensional scaling constant, 0 < α ≪ 1 is a constant incorporated to eliminate singularity and σ VT , C and n are material properties.
The free energy associated to the axonal contribution is defined by the free energy function proposed by Nolan et al. (19): where k 1 and k 2 are material parameters quantifying the increase of stiffness in the axonal direction. The deviatoric fourth strain invariant I 4 * = tr� � * � is defined depending on a structure tensor � , in turn function of a single dispersion parameter ξ and the preferred axon orientation � . Gasser et al. proposed the following compact form (20): The definition of the structure tensor presented in Equation (S.8) has been previously used by Wright et al. to incorporate neural tract alignment through DTI (21). Following the same approach, a functional dependence on FA is proposed for the dispersion parameter ξ: When ξ adopts a value of 1/3, an isotropic distribution of the axons is considered and the structure tensor is spherical. When ξ adopts a value of 0, an ideal coalignment of the axons is considered and the structure tensor reduces then to � = � ⨂� . Thus, the anisotropic behaviour of the white matter arising from axonal orientation and axonal dispersion can be taken into account in the constitutive modelling through � by connecting experiments with modelling.
Using the methodology followed in Garcia-Gonzalez et al. (22), and noting that = , the following expression for the isochoric Cauchy stress tensor is obtained: where dev( * ) and dev( * ) are the deviatoric parts of the distortional left Cauchy-Green deformation tensors * = * ( * ) T and * = * ( * ) T .
Finally, the volumetric response of the tissue is defined by using the Tait equation of state (EOS), commonly used to model fluids and brain tissue (23): where B and Λ o are material constants and ρ o is the reference density. The parameter B can be computed through its relation with the reference bulk modulus K o : This constitutive model was calibrated against experimental data for grey and white matter with a good agreement between experiments from literature and model predictions for both quasi-static and very high rate loadings, see Fig. S5 (13). Note that the strain rate sensitivity is accounted for by the isotropic response of the glial matrix. Also, as observed elsewhere (13,16), anisotropy due to axonal orientation is not relevant at high rates since the stiffening effect due to strain rate effects plays the predominant role in the mechanical behaviour of the tissue. The anisotropic parameters were defined following the work of Wright et al. (21) by assuming the same ratio of the fibre reinforcement parameter to the quasistatic shear modulus (k 1 /μ m ) as in the fit proposed by Velardi et al. (12). Finally, the volumetric parameters were taken from the work of Moore et al. (23). The model parameters for both grey and white matters are provided in Table S1. Note that, due to the lack of available data for human or rat tissue at the high strain rates reached during blast events, the model was calibrated with porcine tissue. Therefore, although the tendencies and overall response of the tissue can be assumed the same, some quantitative differences are expected.

Skin/fat and scalp
The mechanical behavior of skin/fat and scalp was defined herein following Equations (1) and (2) (main text) for the volumetric contribution and using a Neo-Hookean strain-energy function the isochoric stress response: where μ is the shear modulus. The expression for the isochoric Cauchy stress tensor can be obtained by The material parameters for skin/fat and scalp tissues are provided in Table S2. These constants were taken from the literature (23)(24)(25).

Skull
The mechanical response of the skull was described by several authors as isotropic and linear elastic as a first approximation (24,26,27). However, its shock response has also been previously described with nonlinear model, e.g., Moore et al. (23) defined the volumetric response of the skull by using a Mie-Grüneisen EOS. Here, as the deformations observed in the simulations within the skull are very small, the volumetric stress response of this tissue was simply defined by a linear approximation as: The material parameters for skull were taken from previous work and are provided in Table S3 (24).

CSF and ventricles
The isochoric mechanical behavior of CSF and ventricles was considered as driven by its water content and thus characterized by a Newtonian viscosity: = 2η̇′ (S. 16) where η is the dynamic viscosity and ̇′ is the deviatoric strain rate.
The volumetric response is defined here through the Mie-Grüneisen EOS: where ϛ = 1 − ρ o /ρ is the nominal volumetric compressive strain with ρ o as the initial density. c o is the speed of sound in water, s is the slope of the u s -u p curve in the Hugoniot formulation, where u s and u p are the shock and particle velocities, respectively. Γ o is the Grüneisen coefficient and E m is the internal energy per unit mass.
The material parameters for CSF and ventricles were identified from literature and are provided in Table  S4 (24).

Air
The isochoric mechanical behavior of air was defined with a Newtonian viscosity following Equation (S.16), while the volumetric response was defined by the ideal gas EOS: where P a is the atmospheric pressure, R the gas constant and T the absolute temperature. To complete the definition of the volumetric response, the specific energy E is defined as: where E o and T o are the specific energy and temperature in the initial state, and C v is the specific heat at constant volume.
The material parameters for air were taken from the literature and are provided in Table S5 (28).

Injury criteria
The maximum of the following quantities for the whole duration of the simulations was calculated. The simulation time was taken to be long enough so as to ensure that the maximum was reached.

Pressure
The pressure criterion is based on the maximum value of P, Equation (1) of the main text, reached in each region during the deformation process.

Equivalent strain
The equivalent strain is defined as: is the Green-Lagrangian strain tensor.

Shear energy rate
The shear energy rate per unit volume Ψ̇i so is the time derivative of Ψ iso . In this work, the deviatoric energy was computed following Equation (S.1).

Volumetric energy rate
The volumetric energy rate per unit volume Ψ̇v ol is the time derivative of Ψ vol , defined by Equation (S.11):

Axonal stretch
The axonal stretch can be defined from the fourth strain invariant as:

Axonal stretch energy rate
The energy rate from axonal deformation Ψ̇a xon is the time derivative of Ψ axon , defined by taking into account the evolution of the product of the stress tensor with the rate of deformation tensor, projected along the axonal direction: where � = and ̃= are, respectively, the corotational Cauchy stress and the corotational rate of deformation tensor, and where d is the symmetrical part of the total velocity gradient l = ̇− 1 . = − is the rotation tensor, where U is the stretch tensor.

Simplified FEHM
The simplified FEHM was designed as a combination of concentric spheres representing: skull, CSF, brain and ventricles. In addition, another layer was incorporated to introduce a protective shield on the head (see Fig. S7).

Role of acoustic impedance and impulse mitigation for blast protective devices
When a blast wave reaches the head-shield interface, the incident stress σ i is conserved into a reflected stress σ r and a transmitted stress σ t : Similarly, by compatibility at the interface, the transmitted particles velocity v t = σ t �E h ρ h can be obtained from the incident particles velocity and the reflected particles velocity v r : where E h and ρ h are respectively the Young's modulus and the density of the head (assuming the head homogenised).
By combining Equations (S.25) and (S.26), the transmitted stress can be expressed as: where E s is the Young's modulus of the shield and ρ s is the density of the shield. The reduction in the stress amplitude can be thus controlled through the acoustic impedance of the shield material �E s ρ s .
Assuming a perfectly-elastic collision as a first approximation, the transmitted impulse I h can be obtained from the combination of the conservation of momentum and the conservation of energy equations as: where m h and m s are the masses of the head and the shield, respectively, V h is the head velocity after the blast and I o is the impulse imparted to the shield (29). Therefore, impulse mitigation can be reached by increasing the shield mass. Alternatively, energy dissipation can be reached by additional consideration of inelastic deformation in the shield. Fig. S1: Forebrains of Sham and Blast animals at 24 hours post-injury (coinciding with peak urine 3-HPMA elevation) were subdivided into 53 sub-regions for acrolein-lysine adduct immunoblotting analysis. Overall, 8 out of 53 regions demonstrated significant increases in acrolein-lysine adduct immunoreactivity after mild b-TBI when compared to sham-injured animals. Areas with significantly higher oxidative stress relative to Sham included those containing the OFC and agranular insula bilaterally (A4 + A6), the anterior striatum (B5), bilateral somatosensory cortices (C1 + C3), the caudal hypothalamus + rostral midbrain (E8), and areas containing the auditory cortex + temporal association cortex + hippocampal CA2 bilaterally (F4+F6). The anatomical positions of these regions are illustrated in Fig. 1 of the main text. * p<0.05. Unpaired t-tests. Data presented as mean ± SEM. N=4/group.

Fig. S2:
Behavioural facial hypersensitivity following mild bTBI. Mechanical periorbital allodynia was assessed before mild bTBI induction and 8 days post-mild bTBI. Before mild bTBI induction, both sham and blast groups show no difference in facial mechanical threshold; however, following mild bTBI the blast group shows a significant decrease in mechanical threshold compared to sham animals, indicating facial hypersensitivity compared to the sham group. *p<0.05. One-way ANOVA and Tukey test. Data presented as mean ± SEM. N=8/group.

Pre-blast
Mechanical threshold (g)  Protective device Skull Brain CSF and Ventricles