Simulation of mTBI Utilizing White Matter Properties from MRE

Tissues of the brain, especially white matter, are extremely heterogeneous with constitutive response varying spatially. In this paper, we implement a high-resolution Finite Element (FE) head model where heterogeneities of white matter structures are introduced through Magnetic Resonance Elastography (MRE) experiments. Displacement of white matter under shear wave excitation is captured and the material properties determined though an inversion algorithm are directly used in the FE model. This approach is found to improve model predictions when compared to experimental results. In the first place, responses in the cerebrum near stiff structures such as the corpus callosum and corona radiata are markedly different compared with a homogenized material model. Additionally, the heterogeneities introduce additional attenuation of the shear wave due to wave scattering within the cerebrum.


Introduction
With approximately 2.8 million cases reported in the United States [1], traumatic brain injury (TBI) -commonly caused by a direct blow or impulse to the head -remains a pressing concern for study. Clinical results show that damage to the brain in blunt head injuries has a tendency to occur in regions with highly organized axon tracts such as the brainstem and corpus callosum [2][3][4][5][6][7]. These regions serve as pathways to other parts of the brain, meaning damage to them is potentially more dangerous. In this work we introduce a heterogeneous material description of white matter structures to our high-resolution Finite Element (FE) model to account for the local differences in mechanical response between different regions of the brain.
The FE method is commonly used to determine the mechanical response of brain tissue in order to develop improved diagnostic tools and protective measures to reduce the prevalence of TBI. The Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 17 July 2020 doi:10.20944/preprints202007.0366.v1 accuracy of these FE models is highly dependent on the accuracy of the material model. To date, most finite element models have utilized mechanical properties averaged over large portions of the brain, as reviewed in [8], thus ignoring potentially significant effects of local structures. In actuality, the tissues of the brain are heterogeneous, their constitutive response varying from location to location. This is most noticeable for white matter due to the presence of axons with diverse orientations. Overall, the shear stiffness of white matter is 1.2-2.6 times higher than that of gray matter [9]. Locally, white matter tracts with highly oriented fibers such as the corpus callosum and corona radiata have material properties vastly different from other regions. Johnson et al. [10] determined that global white matter was softer on average than either the corpus callosum or the corona radiata. This can be explained by considering the structure of these regions. The corpus callosum provides more structural rigidity than the superficial white matter due to the presence of highly aligned fibers. The fan-like structure of the corona radiata exhibits similar behavior, but to a lesser degree as the fibers are not as highly aligned.
As such the corpus callosum was found to be approximately 11% stiffer than the corona radiata.
The brainstem is another structure with a high level of heterogeneity. Arbogast and Margulies [7] investigated the prevalence of trauma observed in the brainstem after head injuries. They determined that the brainstem was 80-100% stiffer than the cerebrum and concluded this regional stiffness variation is one reason for the selective vulnerability of this region to rotational motion. FE models that utilize homogenized white matter material properties have no way to resolve these local features.
A very useful tool to measure heterogeneity in-vivo is Magnetic Resonance Elastography (MRE) [11] where the head is excited with shear waves to quantitatively assess the local mechanical properties of brain tissue. This technique has been successfully applied to a variety of different applications such as decrease in whole-brain stiffness with age [12] and as a result of certain neurodegenerative diseases  [10]. In order to properly resolve the high-frequency waves generated during such vibrations, the finite element mesh must be properly refined, as argued in [8].
The mechanical properties used in this work are reconstructed at the same spatial resolution as the displacement data captured during the MRE process. This, in turn, generates a finite element mesh that has more than sufficient resolution to accurately capture the dynamic shear wave propagation during impacts. To the best of our knowledge, this is the first attempt to include heterogeneity of brain tissue via MRE in a high-resolution FE model.

MRE Aquisiton and Inversion
The heterogeneous properties of the white matter are taken from the work by Johnson et al. [10].
A brief overview of the acquisition and inversion is presented here. A more detailed discussion of this by iteratively updating the material property description, given by θ. Here, u m i is the measured displacement amplitude at the location i; u c i (θ) is the computed displacement amplitude sampled at the same point; and * indicates the complex conjugate.
Following the development in [22], a Rayleigh damping model is used to represent the material response of brain tissue under time-harmonic conditions. The motion amplitude field u, is calculated from Navier's equation for an inhomogeneous, locally isotropic linear elastic medium where λ is the first Lamé parameter; G is the second Lamé parameter, or shear modulus; ρ is the density; and ω is the activation frequency. The Rayleigh damping model introduces the complex-valued shear modulus and density to account for attenuation related to both elastic and inertial forces -where the imaginary shear modulus includes damping effects due to inertial forces, while the imaginary part of Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 17 July 2020 doi:10.20944/preprints202007.0366.v1 the density includes damping related to inertial forces. Including inertial damping effects -something that commonly used viscoelastic models do not include -allows for better characterization of material response when performing the inversion. We use the notation G and G to denote the storage and loss shear moduli respectively, which are the real and imaginary part of the complex valued shear modulus i being the imaginary unit. Due to the nearly incompressible nature of brain tissue, we assume that λ is very large compared to the shear modulus. The damping ratio can be determined as well, which physically describes the level of motion attenuation in the tissue.
The distribution of storage and loss moduli within the white matter in the imaging resolution (2×2×2 mm 3 ) is presented in Figure 1 for the coronal, sagittal and horizontal planes. The relative stiffness of both the corpus callosum and the corona radiata can be clearly observed. The average values and the standard deviation for the global white and gray matter is presented in Table 1 as well as those for the corpus callosum and corona radiata.

Finite Element Mesh Generation
In addition to the MRE image acquisition, T1-weighted anatomical images with a resolution of 1 × 1 × 1mm 3 are generated for the purposes of mesh segmentation. Image voxels are directly converted to eight-noded hexahedral elements through a custom code, thus preserving the same spatial resolution as the MRI scans. The resulting mesh consists of approximately 2.2 million elements.
The mesh is segmented into the following tissue types: scalp, skull, cerebrospinal fluid (CSF), gray matter and white matter. Each element is assigned a different material definition based on the results of the segmentation, as detailed in section 2.3. Features that are below the imaging resolution such as membranes, blood vessels, bridging veins, and draining sinuses are excluded from the model.
Since T2-weighted scans were not collected, the segmentation of the interfaces between tissue types is negatively affected. We perform a mesh smoothing operation at these interfaces to minimize this effect.
We utilize a volume-preserving Laplacian smoothing algorithm as outlined in [17]. Smoothing has the added benefit of eliminating numerical artifacts that may manifest from jagged edges along interfaces.
We ensure traction and displacement continuity at all material interfaces.
Our MRI voxel-based approach produces meshes which more realistically model the complicated folding structure of the cerebral cortex (i.e. gyri and sulci) as well as the differentiation of gray and white matter of brain tissues. Additionally, in order to accurately resolve the shear wave motion within the brain, the mesh must be sufficiently refined [8]. For example, using the typical speeds of shear wave propagation in brain tissue, c T ≈ 5m/s, and the frequency of the vibration as 50Hz, we arrive at a minimum element size of 10mm (using a conservative choice of 10 elements per wavelength [23]).
For vibrations of higher frequency or more transient loading (such as the cases for impact loads), this requirement is even more stringent.

Material Properties
As discussed above, we assign different material definitions to each voxel based on the result of the mesh segmentation. Due to the presence of highly oriented tracts of myelin-sheathed axonal fibers in white matter, significant heterogeneity exist in this region, especially within the corpus callosum (CC) and corona radiata (CR). On the other hand, gray matter is made up of cell bodies and supporting vascular networks that can be assumed to be isotropic and homogeneous [10,24]. This assumption Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 17 July 2020 doi:10.20944/preprints202007.0366.v1 For homogeneous tissues, the choice of material properties is determined from data used commonly in the literature, presented in Table 2, and is the same as that in our earlier works [17,19].
The skull is assumed to be linear elastic and modelled as a single-layer structure with homogenized properties given in [25]. We assume a linear viscoelastic material model for the shear response of the brain tissues. The standard linear solid model is chosen for its simplicity, with the shear relaxation modulus given by Here, G 0 is the short-term shear modulus, G ∞ is the long-term shear modulus and β is the decay factor. .94m/s. The impact was along the specimen's mid-sagittal plane in an anterior-posterior direction.
The skull was rotated as depicted in Figure 2a. The resultant input force time history is given in Figure   2b. Intracranial pressure history was recorded during the impact event.

Simulation of Impact
The recorded impact force from the frontal cadaveric impact experiment conducted by Nahum, as discussed in the preceding section, is directly applied to the model in the form of a distributed load with the peak pressure input of 4.37 MPa on the frontal region of the skull. The impact pulse lasts about 9 ms and the simulation is run for 15 ms. Rotation of the model is permitted as we the base of the skull is not constrained. We find that this free boundary condition gives better correlation of coup pressure to experimental results [17]. The FE simulations are carried out using Abaqus/Explicit. The pressure-time history for the Nahum loading is presented in Figure 3 at the coup impact site.
We plot the results for both the homogeneous and heterogeneous model. The model predicts tensile pressure for nearly the whole duration of impact. The peak pressure predicted by the heterogeneous Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 17 July 2020 doi:10.20944/preprints202007.0366.v1 material model more closely matches that from Nahum's experiments than the homogeneous model.
In both cases the peak pressure occurs roughly at the same time -indicating very similar wave speeds. We plot the contours of the von Mises stress distribution on the sagittal plane due to the frontal impact in Figure 6. The spherically convergent structure of the shear wave is clearly observed. The wave attenuates as it travels inwards and eventually dissipates. Reflections of wave due to scattering from heterogeneous white matter structures can also be observed at later times.
Next, we consider three distinct points along the sagittal plane, as depicted in Figure 7. The points are chosen within regions of strong heterogenities due to the presence of highly aligned axon tracts, such as the corpus callosum and corona radiata. The differences in mechanical properties of these regions are given in Table 3. We see that the material phases at these points are relatively stiffer than   Figure 7 is affected accordingly. We find that the difference in peak pressure response is proportional to the difference in shear stiffness between the homogeneous and heterogeneous models. However, the time at which these events occur is not significantly affected. This indicates that the pressures in regions of high stiffness within the brain are over-estimated in the homogeneous models. In summary, relative to the MRI-based model, the new MRE-based heterogeneous model more accurately predicts the local response within the white matter by taking into account the differences in tissue stiffness of local white matter structures.
A few points are in order regarding the qualitative differences between the shear modulus of different regions in our model. Globally, the white matter is found to be approximately 32% stiffer than the gray matter. In general, the white matter properties in local regions differ significantly from the average ones. For instance, the storage modulus G is significantly lower in the rest of the white matter than within the corpus callosum and the corona radiata. This is quite logical given the fact that the corpus callosum contains highly oriented, tightly packed axon tracts. The corpus callosum, in turn, is stiffer than the corona radiata, again evident from the composition of the corona radiata which contains axon fibers that fan out and are not as highly aligned as the corpus callosum.  Table 3. Difference of material properties and peak pressure response for three distinct points (indicated in Figure 7) along the sagittal plane within the white matter. and DTI measures correlate well with each other within the corpus callosum -not surprising since they are both sensitive to the underlying tissue microstructure. They hypothesize that these measures are highly dependent on axon diameter since larger axons provide greater structural rigidity to the tissue [38]. Within the corona radiata, however, the correlation is not as significant. The corona radiata comprises fiber tracts that fan towards the cortex and contain numerous crossings [33] which are not captured well by DTI [39]. This has been hypothesized as the reason for the poor correlation within the corona radiata. More work is needed to determine the differences between these two methods when used within FE models.

Stochastic Wave Propagation
The highly heterogeneous structure of the brain tissue introduces wave scattering that competes with wave amplification due to spherically convergent implosion. Following the development in [40], we investigate this effect by considering the theory of wavefronts. For the case of one-dimensional wave motion, we assume that a compressive load produces a shock wavefront that propagates from a disturbed domain to an undisturbed one with a speed c T . The initial conditions can be given as where H is the Heaviside function.
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 17 July 2020 doi:10.20944/preprints202007.0366.v1 Assuming a plane wave in a homogeneous medium, we have the dynamic compatibility condition denotes the discontinuity in a function across the boundary of two materials, σ is the shear stress, ρ is the mass density, u is the displacement normal to the direction of wave motion, and c T is the transverse wave speed. The linear viscoelastic stress-strain relation for a process that started at time t = t + 0 is where ε is the shear strain. We can derive the relationship for the wavespeed as c T = G(0)/ρ, where 212 G(0) is the glassy modulus. Following the derivation in [40], we obtain the equation governing the 213 evolution of the discontinuity of τ at the wavefront as: On account of the initial conditions above, the solution of equation 3.2 is Given that G ,t (0) ≤ 0, and G(0) > 0, the stress jump exhibits exponential attenuation and has a 216 tendency for blow-up as r → 0. As our simulations here and in [19] demonstrate, the attenuation is sufficiently strong so that the imploding waves generated from transient impacts do not blow up into a singularity at the head center.
The impact results not only in a fast pressure wave, but also in a slower shear wave. Due to its relatively low shear modulus, brain tissue deforms much more easily in shear than in dilatation mode.
Thus, the shear wave is potentially more damaging. Recall that the spherically convergent shear wave patterns are observed even in the case of homogeneous material description, c.f. [17]. The attenuations of pressure along the sagittal plane for both the homogeneous and heterogeneous models are presented in Figure 8. It is clear that the attenuation is greater in the heterogeneous model as predicted. This is consistent with studies of transient wave propagation in elastodynamics of random media [41,42].

Limitations
This section discusses a few limitations of our methodology, some of which are likely to be improved in future works. First, features that are below the imaging resolution such as membranes (in particular, the falx cerebri and tentorium), bridging veins and blood vessels are excluded from the current model. However, we argue that the inclusion of heterogeneous MRE-derived parameters in white matter tissue alleviates some of the drawbacks of excluding these features. Additionally, since the MRE inversion was performed assuming locally linear viscoelastic material model, we have chosen to utilize linear viscoelasticity for white and grey matter. Research has shown [43] that the the shear modulus of the brain tissue is much lower than that used in this work. We aim to perform an optimization study to tune the values chosen here to better match experimental results.  [10]. While many FE models employ homogenized or averaged, mechanical properties to approximate constitutive relations of brain tissues, our approach allows us to investigate response due to local structures within the white matter. We note that both the corpus callosum and corona radiata are significantly stiffer than overall white matter, with the corpus callosum exhibiting greater stiffness and less viscous damping than the corona radiata.