A Nonlinear Reduced-Order Model of the Corpus Callosum Under Planar Coronal Excitation

Approximately 2.8 million people in the U.S. fall victim to traumatic brain injury (TBI) annually [1–3]. Such a high number of incidents have become a major socioeconomic concern. This calls for a deep and detailed understanding of the mechanisms and physics behind TBI. To achieve this goal, some researchers have tried to view brain injuries as a result of brain movements and deformations. In 1943, observing that the brain behaves as an elastic body, Holbourn proposed a mechanical system to represent and model the mechanics of brain tissue damage through the theory of elasticity. He postulated that large-strain deformations in this tissue are more dominant due to rotational rather than translational acceleration excitations [4,5]. Following the path of examining representative reduced mechanical models, Kornhouser studied brain deformations and head motion as second-order mechanical systems over a wide range of impact conditions and suggested isodisplacement curves, also known as sensitivity curves, for relative brain displacement as an injury classifier [6].

Many scientists have studied the dynamical behavior of the brain-skull system to quantify the severity and predict the occurrence of TBI. In doing so, kinematics-based criteria, such as the head injury criterion [7,8] and brain finite element (FE) model-based criteria such as cumulative strain damage measure [9], have been developed. Others have taken advantage of the recent advances in imaging techniques and hypothesized the axonal injuries to be one of the main contributors of mild TBI (mTBI) [10,11]. Furthermore, large strains in regions of the brain with high density of axonal fibers such as the periventricular region have often been associated with concussions [12,13]. Studying such fiber tracts, the corpus callosum (CC) has been implicated in TBI pathology. Using a FE model, Giordano et al. determined an axonal strain threshold of 0.07 for CC and 0.15 for brainstem as mTBI thresholds, which is in agreement with experimentally derived injury thresholds [14]. Furthermore, they suggested that this mTBI threshold predicts the occurrence of concussion more accurately than other criteria such as head injury criterion and cumulative strain damage measure. In a later study, Hernandez et al. developed a six degree-of-freedom (DOF) model to evaluate and predict injury risks. With this model, they studied 513 cases of head impact measurements, two of which caused concussion, and showed that peak principal strain in the CC followed by peak rotational acceleration magnitude and the head impact power can be potentially used as the strongest predictor of concussion [15]. To understand how an external head impact is converted into regional brain responses and affects the biomechanics of substructures such as the CC, FE models of the human head play a crucial role. These models are often validated based on displacement fields acquired from cadaver head in situ under excitation [16]. Others use ex vivo samples to further improve the accuracy of FE models [17,18]. An alternative solution to measure displacement fields and validate FE models is using novel imaging techniques such as magnetic resonance elastography [19,20], amplified magnetic resonance imaging [21], and tagged magnetic resonance [22]. Kleiven et al. implemented a KTH brain model and calculated strain, strain rate, and stress in the CC as potential indicators of concussion [12]. Other FE models such as the Worcester Head Injury Model (WHIM), a FE brain injury model in which the brain tissue is modeled as an isotropic, homogeneous material using a second-order Ogden hyperelastic model) [23,24], and the Global Human Body Models Consortium (GHBMC) [25] have also been used to locally calculate tissue-level mechanical responses of the brain. Studying the dynamics of the brain under head impact using a WHIM FE model, Abderezaei et al. found local nonlinear effects near the CC and suggested a correlation between this phenomenon and TBI pathology [26]. Distinct modal behavior of the CC at low frequencies was also reported to be causing high strain levels in the brain [27]. Other studies have used imaging techniques and identified the CC as an important region in TBI pathology, where abnormalities in the fractional anisotropy of white matter (WM) tracts, especially at the CC, thalamus and hippocampus, have been reported in patients suffering from TBI [28].

Although they provide accurate and thorough simulations of the response of the human brain under mechanical loads, these brain FE models have limitations since they incur a significant computational cost to simulate a single head impact [9,29,30], therefore rendering these methods impractical in situations where a large number of head impacts are to be simulated, such as impacts occurring over the span of a sports season [31,32].

In order to avoid the high computational cost of using FE models for on-site simulations, several reduced order models have been developed to study the dynamics of the brain by modeling brain tissues as mechanical systems such as viscoelastically supported rigid masses. The reduced order model (ROM) proposed by Alem [33] was a 12DOF linear system of oscillators coupled by springs and viscous dampers. He hypothesized that the damage to the brain is due to its translational and rotational motions relative to the skull. Based on his model, Alem predicted high relative rotational strains in impact cases compared to relative axial strains. Low and Stalnaker [34] developed a linear 3DOF system validated by experimental data to model the dynamics of the brain under coronal or angular accelerations with the aim of providing a way to establish injury criteria. More recently, Laksari et al. showed that overall dynamics of the brain-skull system can be modeled as a single DOF under-damped linear oscillatory system to represent the motion of the brain under low-impact scenarios [35]. Their analysis showed that typical head exposures are able to drive the brain close to its natural frequency which can potentially lead to high-amplitude brain-skull relative motion. Also recently, Gabler et al. [36] developed a new injury metric, diffuse axonal multiaxis general evaluation (DAMAGE), based on a 3DOF linear oscillatory system, predicting maximum brain strain. Neither of these models, however, focuses exclusively on the deep white matter region of the brain, often implicated in TBI [27,37], thereby missing the local, potentially important dynamics within these regions. A potential alternative approach to avoid the high computational cost of using FE models is to utilize precomputed brain response atlases (pcBRA) [32,38]. These are an attempt to simplify and speed up the process of studying different impact profiles and the response of human brain, while maintaining the sophistication of the FE models.

In this paper, inspired by the configuration of the human brain in the coronal plane, we develop a geometrically nonlinear ROM that focuses on the dynamics of the deep WM region of the brain and more specifically on the CC of the brain, by modeling the dynamics of each of the brain hemispheres and the CC as coupled nonlinear spring-mass system of oscillators. First, we correlate the configuration of the human brain in the coronal plane with that of the ROM, and then express the dynamics of the ROM in terms of an initial-value problem. Next, we test the accuracy of the model by comparing its responses to the responses of the advanced FE model of the human brain proposed by Ji et al., and the Worcester Head Injury Model (WHIM) [23,24] which were reported in Ref. [26]. This is followed by comparison of concussion thresholds, i.e., isostrain curves, obtained from the ROM with those from the literature. Next, by using numerical modal analysis techniques, we show that the nonlinear behavior of the proposed ROM is in close agreement with that of the FE model. Finally, we inspect the dynamics of the CC exclusively under excitations of different intensity to determine qualitative changes in its behavior.

Our purpose in this paper is to present a nonlinear dynamical reduced-order model, which enables us to study and predict the behavior of the deep WM region of the human brain in the coronal plane under impact-like excitation (Fig. 1). We chose the coronal plane since many studies have shown that rotations in this direction impose high strains on the CC [26,37,39]. In the proposed ROM, we neglected the dynamics within the hemispheres and considered them to be lumped masses located at their respective centers of mass. These representative point-masses were coupled to each other through two pinned viscoelastic links, representing the CC. Figure 1 illustrates the schematics of the ROM superimposed on a cross section of the deep WM region of the human brain in the coronal plane and its detailed configuration.

where the vectors $y=[rc,θc,rpr,θpr,rpl,θpl,rr,rl]T$ and $p=[M,m,mp,κg,λg,κp,λp,kl,dl,kr,dr]T$ contain the state variables and parameters, respectively, with F representing the half-sine moments applied to the hemispheres corresponding to rotational acceleration with amplitude of $αM$ and duration of $td$ as shown in Fig. 1. Table 1 lists all the state variables and parameters and explains what each represents.

The procedure of obtaining the governing equations of motions is explained in detail in the Supplementary Information (SI), available in the Supplemental Materials on the ASME Digital Collection, portion of this paper. Moreover, each element in the parameter vector, p, was identified either through a series of system identification analyses or obtained directly from the configuration of the human brain finite element model [23,24]. A system identification method—in this case the restoring force surface method—was employed to determine the unknown parameters of the proposed ROM considering the datasets reported in Ref. [26]. These datasets were obtained for the case where a half-sine acceleration was applied to the center-of-mass of the head in the coronal plane. The dataset corresponding to the half-sine acceleration excitation with $αM=5 krad/s2$ and $td=5 ms$ was used for system identification analyses. After obtaining the unknown parameters, the responses of the left and right hemispheres along with the response of the CC were reproduced using the ROM and compared against those from the FE model simulations. The quality of the comparison between the response of the ROM and that of the FE model was quantified using the coefficient of determination, R². Detailed descriptions of the system identification analysis are provided in the SI discussion, available in the Supplemental Materials on the ASME Digital Collection, as well.

An important feature of the proposed ROM in this work is its nonlinear behavior due to geometric and kinematic effects. The nonlinearity in nonlinear systems of oscillators causes the nonlinear normal modes to interact, e.g., exchange energy. Note that, in linear time-invariant systems of oscillators, no two governing normal modes can interact and exchange energy. Modal energy exchange is one of the many ways that nonlinearity in dynamical systems can manifest itself. In order to determine the contribution of the governing normal modes of the system in the dynamics of our ROM, for each case of simulation of the ROM excited by a single half-sine pulse with duration $td∈[5,30] ms$ and amplitude $αM∈[0.01,15] krad/s2$⁠, we used the discrete proper orthogonal decomposition (POD) method. Discrete POD is the singular value decomposition of a matrix, which aims to reduce a high-dimensional set of data to one of lower dimension by simplifying the process of revealing the phenomena of interest that are originally embedded in the data-set [26,35,40–42]. By applying singular value decomposition to a displacement data-set recorded for a certain time duration—in the case of periodic or quasi-periodic responses of a system at least half the period of the mode with the lowest frequency—we are able to identify an orthogonal set of eigenvectors spanning the vector-space where the displacement at any time-instant is represented by a vector. In oscillatory dynamical systems the eigenvectors are directly related to the governing normal modes of vibration, along with their associated singular values, quantifying their participation in the dynamics.

Results presented in this section are divided into two categories. In the first section, we verify the accuracy of the identified parameters from the system identification analyses mentioned in Sec. 2.1, as well as the accuracy of the ROM. This is done by direct comparison of the time-responses of the two hemispheres and their corresponding frequency content, obtained from our ROM and from the FE model [23,24,26].

In Sec. 2, we investigate the nonlinear behavior of our proposed ROM by computing the modal energy distribution obtained through application of POD to the time-responses of the ROM for different values of $αM$ and $td$⁠. Next, we show the isostrain curves calculated by our proposed ROM and compare them against those obtained by Margulies and Thibault [43] and the experimental TBI data reported by Hernandez et al. [15]. Finally, we study the nonlinear behavior of the CC.

Following the system identification procedure described in Sec. 2.1, we reproduced the response of the ROM and compared it against the response of the FE model. Figure 2 shows the comparison between the response of the left hemisphere obtained by the ROM and the response of the center of mass of the left hemisphere determined by the FE model, with coefficients of determination of $R2=0.70$ and $R2=0.45$ in the x and y directions, respectively. It should be noted that, since the response of the right hemisphere of the brain is not qualitatively different from that of the left hemisphere, here we only show the response of the left hemisphere. More comparisons to ensure the accuracy of the identification process are reported in the SI available in the Supplemental Materials on the ASME Digital Collection.

The brain as a biomechanical system ranks among the most complex natural systems. This is due to its mechanical properties, complex geometry, and interaction with its surrounding skull. From a dynamical point of view, the aforementioned properties lead the brain to exhibit strongly nonlinear behavior originating from its nonlinear and nonhomogeneous interactions with the skull and other nearby organs, its inhomogeneous and nonlinear material composition, and its geometry [24,35]; therefore, when studying the dynamics of the brain, it is critical to take such factors into account. However, there is a tradeoff between taking these properties into account and the simplicity of the representative dynamical ROM.

As noted earlier, the nonlinearity of an oscillatory dynamical system can cause modal interactions, i.e., energy exchanges among governing normal modes of the system. In order to demonstrate this, the percentage of the input energy allocated within the first two proper orthogonal modes (POMs) governing the response of the FE model and the ROM for different amplitudes and durations of excitation, was evaluated to analyze the capability of our model to capture the intrinsic nonlinearity of the brain. It was observed that as the intensity of the excitation increases the first POM begins exchanging energy with higher POMs, a typical behavior observed in nonlinear dynamical systems (Fig. 3). More specifically, in the extreme case of $αM=15 krad/s2$ and $td=30 ms$⁠, one can see that the first POM possesses only 60% of the input energy to the system, and 30% of the input energy is allocated within the second POM. On the other hand, for small $αM$⁠, i.e., $αM<2 krad/s2$⁠, almost all the input energy is captured by the first POM. Comparing the two extreme cases and how energy is distributed between the two POMs, it becomes obvious that the nonlinear behavior in the system becomes more pronounced as the input energy increases.

Next, we analyzed the possible correlation between the energy content of the first POM corresponding to various excitation levels with isostrain curves obtained from the ROM as well as those reported in the literature ([26]; Fig. 4). Additionally, the experimental head impact data of 15 athletes (6 injury, 9 noninjury cases, [15]) was compared against the energy contribution of the dominant POM and isostrain curves. These results hint on the connection between TBI and the appearance of nonlinear effects in the deep WM region of the human brain; this can be inferred also from the contour plot in Fig. 4, and the experimental cases where TBI has occurred. This connection is inferred due to the accumulation of the experimental TBI cases close to the region where the level of the energy ratio contour of the first POM starts to change, i.e., in precisely the regime where the geometric nonlinearity begins to strongly affect the dynamics of the system.

Considering the possible implications of this nonlinearity on the damage mechanism and the importance of the CC in the TBI mechanism, it is essential to characterize how the dynamics of the CC changes as the system makes the transition from linear to nonlinear. Here, the time response and its corresponding frequency content of the CC in the ROM, between low (⁠$αM=100 rad/s2$⁠, $td=5 ms$⁠) and high (⁠$αM=23 krad/s2$⁠, $td=25 ms$⁠) energy inputs, is compared—cf., Fig. 5. Considering the response of the CC for the low input rotational acceleration—cf., Figs. 5(b) and 5(e)—as the linear baseline for the behavior of the CC, a clear change in the behavior of the CC as the input energy increases is observed compared to its baseline behavior—cf., Figs. 5(c) and 5(f). The qualitative difference is more discernible in the wavelet transforms (WTs) of the time-series, which represent the temporal evolution of the frequency content (dark regions in the contours of Fig. 5) of their corresponding time-series and show the generation of high frequency harmonics in the system at approximately 45 Hz.

Due to its biological importance and susceptibility to injury, it is crucial to study the dynamics of deep regions of the brain exclusively. With this in mind, we developed a geometrically nonlinear dynamical ROM of the CC to closely follow the dynamics of the deep WM region of the human brain. We first compared the time response of the center of mass of each hemisphere obtained from simulating the FE model [26] and the ROM, for the case where the head is excited by an ideal half-sine rotational acceleration impulse. The quality of the reproduced responses by the proposed ROM compared to the FE model was found to be $R2=0.70$ and $R2=0.45$⁠, in the x and y directions, respectively. Such relatively low values of $R2$ between the models in the y-direction are due to the apparent phase difference between the two time-series. However, the frequency content of both responses correlates closely when compared against that of the simulated response of the FE model of the brain. These close correlations confirmed and verified the accuracy of our proposed model and the system identification analysis. It should be noted that the ideal half-sine excitation is rather unrealistic and less likely to occur in real-world situations; however, due to its simple and yet broadband nature, it is suitable to employ when examining and studying the transient responses of dynamical systems.

In order to further verify the accuracy of the proposed ROM, we compared the percentage of input energy distributed between the first two POMs for different amplitudes and durations of excitation, determined by the ROM and the FE model [26]. As demonstrated by Fig. 3, the ROM is able to accurately replicate the distribution of the input energy percentages allocated in each of the POMs calculated from the FE model. It was observed that, by increasing the amplitude and duration of the excitation, the POMs begin exchanging energy. This is a clear indication of nonlinear behavior of the system [26]. It should be noted that, quantitatively, the two sets of surfaces in Figs. 3(a) and 3(b) do not match closely. The reason for this mismatch is that the energy allocation percentage surfaces from the FE model were obtained by taking into account the energy that is not localized in the deep white matter region, whereas this is not the case when computing those surfaces for the ROM. Accounting for the energy of the system outside the deep white matter region of the brain when computing the modal energy percentages for the FE model cannot be avoided because of the nature of the POD method. The qualitative similarity of Figs. 3(a) and 3(b), in terms of amplitude-duration dependence of the modal energy ratios, computed for the ROM and the FE model, suggests that similar nonlinearity governs the dynamics of both the ROM and the FE model. As mentioned in Sec. 2.1, the ROM consists only of linear mechanical elements, i.e., linear translational and torsional springs and viscous dampers; hence, from all the possible sources of nonlinearity, such as the nonlinear behavior arising from the geometry, material properties, and boundary effects, the dominant source in this model is the geometric nonlinearity. This leads us to conclude that the dominant nonlinear effect observed in the deep WM region in the ROM and the FE model appear to be due its special geometry. This observation is in line with the reported nonlinearity in the simulated response with the WHIM model, where a shift in the mode shape of the falx, which was perpendicular to the brain's deformation is suggested to be the source of geometric nonlinearity near the CC [19]. Such nonlinear behavior has been reported in other studies, as well. Sabet et al. analyzed the brain's deformation under a rotational acceleration and observed disrupted strain fields near dura matter and CC through tagged magnetic resonance imaging, which suggested local nonlinearity in those regions [44]. There was also evidence of localized modes in the deep regions of the brain where nonlinearity was observed [27].

Having acknowledged the importance of nonlinearity in the ROM, it is of interest to see how this phenomenon affects the isostrain curves—curves in $αM$–$td$ parameter-space that can roughly serve as mTBI thresholds. The obtained curves were compared against the existing isostrain curves in the literature [43] and were viewed along with the energy of the POMs—which implied the strength of the nonlinearity of the ROM for different excitation levels (Fig. 4). An interesting observation drawn from Fig. 4 was that the isostrain curves computed by the ROM or by Margulies and Thibault [43], along with several experimental TBI cases [15], were located in close proximity of where CC begins behaving nonlinearly [26]. This suggests a possible link between the intensity of the nonlinear effects and the likelihood of TBI occurrence. This conjecture is similar to that made by Abderezaei et al. [26], where they found local nonlinearity in the deep WM and reported possible correlation between these concussion thresholds and the volume of nonlinearity in those regions. Furthermore, it was observed that, similar to the findings in Ref. [43], for low-amplitude excitation levels, $αM<500 rad/s2$⁠, the nonlinearity is not dominant in the dynamics regardless of the duration of brain exposure to the excitation. This conclusion originates from the fact that, according to Fig. 4, for peak acceleration levels less than $500 rad/s2$ the POMs do not interact, i.e., the color of the contour which represent the energy localized in the first POM barely changes for $0≤αM≤500 rad/s2$⁠; hence, the nonlinearity in the system can be disregarded. Similar behavior holds for exposure durations less than 5 ms, where the brain is exposed to small amplitude impact energies, and its dynamics is not affected significantly by the nonlinearity, regardless of the amplitude of excitation.

After reporting the existence of nonlinear behavior in the deep WM by studying energy redistribution among the governing POMs and considering the importance of the CC in the TBI mechanism, we investigated the effects of the nonlinearity on the dynamics of this region as the energy of the system changed with increasing input amplitude and duration. Surprisingly, we observed that in concussion-level impact simulations, which according to Fig. 4 results in 5% to 10% strain in the brain hemispheres, there was generation of higher harmonics in the response of the CC as a direct result of the geometric nonlinear effects discussed earlier. Higher harmonics generated due to nonlinearity in an oscillatory system are associated with high frequency oscillations. Hence, for high-intensity excitations the CC oscillates at higher frequencies compared to its surrounding regions. Additionally, high frequency harmonics as high-frequency carriers of energy enable the POMs to exchange energy and possibly cause localized strain concentrations in deep regions of the brain. Such behavior of the CC has previously been reported by Laksari et al. [27] who showed dependence of tissue-level strain patterns on the brain's frequency response.

While the presented reduced-order model is capable of capturing and reproducing the general nonlinear dynamics governing the FE model of the human brain, it possesses certain limitations. For instance, the ROM is modeling the dynamics of the human brain, only in the coronal plane; as such, it neglects the dynamic coupling between the motions of the brain in the coronal plane and those in the sagittal and axial planes. Additionally, due to the nature of reduced-order modeling, the elastic and viscous elements are merely representatives of the effective continuum viscoelasticity of the deep WM tissues. Although such limitations are unavoidable when creating reduced-order dynamical models, they do not diminish the ability of those models to replicate the behavior of their original (intricate) systems and reproduce their responses accurately.

In summary, motivated by the evidence showing that the deep WM of the brain has been implicated in TBI pathologies, and that employing sophisticated FE models to study the mechanism of TBI in a large cohort of patients is intractable due to their large computational cost, we created a geometrically nonlinear dynamical (ROM) of the deep WM region of human brain. The accuracy of the ROM was verified by comparing its performance with that of the WHIM brain FE model, in both time and frequency domains. Next, using numerical modal analysis techniques, we calculated the percentage of the input energy allocated within the first two normal modes governing the dynamics of the ROM and observed good agreement when compared with those computed with the FE model. This observation led to two conclusions: first, the ROM is accurate enough to represent the local dynamics of the deep WM region of the brain, and second, owing to the fact that the only type of nonlinearity built into the ROM is geometric (by construction), this phenomenon originates from the specific local configuration of that region, i.e., geometric sources of nonlinearity. Finally, to investigate possible injury mechanisms, we compared the response of the CC in both time and frequency domains for noninjury (low-energy) and injury (high-energy) cases. It was revealed that in injury cases new high-frequency harmonics were generated. This may imply the possibility of correlation between TBI pathology and existence of high-frequency harmonics as a result of the brain's inherent geometrical nonlinearity.

No benefits in any form have been or will be received from a commercial party related directly or indirectly to the subject of this paper.

Any opinions, findings, and conclusions or recommendations expressed in this work are those of the authors and do not necessarily reflect the views of the National Science Foundation.

We thank professor Sonbai Ji and professor Wei Zhao for providing us with numerical data from their FE model of the human brain.

• •

  National Science Foundation (Grant Nos. CMMI-17-1727761 and CMMI-17-1728186; Funder ID: 10.13039/100000001).