Brain Tissue Response Analysis Based on Several Hyperelastic Models, for Traumatic Brain Injury Assessment

Numerous geometrically simplified models may be found in literature for simulation of the traumatic brain injuries due to the increased intracranial pressure caused by sever translational accelerations of the brain inside the cranium following the impact waves. Some researchers have used more accurate models but employed specific hyperelastic material models. No research has presented a comprehensive comparison among results of various geometric and hyperelasticity models, so far. In the present research, two distinct finite element models and four hyperelastic constitutive models (i.e., polynomial, Yeoh, Arruda-Boyce, and Ogden models) are employed to accomplish the mentioned task. Therefore, the motivation is checking accuracy of the modeling procedure and discussing the results according the traumatic brain injury criteria. In this regard, a realistic skull-brain model is reconstructed in CATIA software based on the MRI scans and employed for optimized mesh generation in HYPERMESH finite element software. Influence of the contact and nonlinear characteristics of the brain tissue are considered in simulation of the relative motions in LS-DYNA finite element code. Time histories of the accelerations and the pressures (von Mises stresses) are derived from ANSYS finite element analysis code. Finally, the responses are discussed based on the available traumatic brain injury criteria and tolerances. Comparisons made with the available experimental results for the four hyperelastic constitutive equations confirm that employing Arruda-Boyce or Ogden models may lead to inaccurate or even erroneous results. On the other hand, the polynomial model is the most accurate model but underestimates the injury probability and may be used with care.


Introduction
The relative translational brain motions due to severe accelerations and decelerations following the impact waves, lead to increased intracranial pressure gradients that are exerted on the brain, cerebrospinal fluid (CSF), and the brain's blood supply vessels. These pressures can lead to traumatic brain injuries, blood vessels damages, and restricting blood flow to the brain and brain trauma and consequently, to fatal conditions if they exceed 40 mmHg in adult persons [1]. Some of the frequent origins of the damages are the relative motions of the brain with respect to the skull, striking and bouncing of the parenchyma against inner skull protrusions, cavitation phenomena induced by negative pressures, and rupture of the bridging veins, axonal fibers and vascular tissue [2,3]. Yue [4,5] studied dynamic characteristics of the human skull-dura mater system, using a hollow sphere model and finite element codes. The model was used to determine the resulting deformations of the human skull due to variations of the intracranial pressure. The thin-walled skull was simulated by composite shell elements. The viscoelasticity of human skull-dura mater system was studied and analyzed by Maxwell's model. The important mechanical characteristic of the cancellous bone and dura mater is viscoelasticity [6,7]. Willinger et al. [8] and Ding et al. [9] determined mechanical properties and specifications of the compact bone, cancellous bone, and dura mater.
Early finite element models used a simple fluid-filled spherical shell to represent the skull/brain complex but they have been developed gradually. Zhang et al. [10], Willinger et al. [11], Kleiven and von Holst [12], and Horgan and Gilchrist [13] have recently developed more complicated 3D models based on anatomical drawings and medical images such as computer tomography (CT) and magnetic resonance imaging (MRI) scans. Liu et al. [14] presented a transparent physical head model with air bubbles to simulate the brain cavitation phenomenon in head decelerating impacts. The transparent skull model was generated based on a real human skull through the turnover formwork technique, and a transparent gel was used to substitute the brain tissue. Recently, Chen and Ostoja-Starzewski [15] presented a 3D finite element model of the human head that accounts for important geometric characteristics of the various components within the human head through an efficie1nt magnetic resonance imaging voxel-based mesh generation method. El Sayed et al. [2] presented a biomechanical model for simulation of the traumatic brain injury and soft tissues damages to reproduce axonal damage and cavitation injury through inelastic deformations, examining frontal and oblique head impacts with external objects. The material response was divided into elastoplastic and viscoelastic components, including rate effects, shear and porous plasticity, and finite viscoelasticity.
Some researches considered simple hyperelastic or finite viscoelastic models [16][17][18][19][20][21][22][23]. Elastic and hyperelastic properties of the brain gray and white matters were characterized by Kaster et al. [24], based on the forcedisplacement data of the tissues for 25 different brain samples, using an indentation apparatus. Young's modulus and the hyperelastic parameters corresponding to the commonly used Polynomial, Yeoh, Arruda-Boyce, and Ogden models were obtained. Recently, Post et al. [25] studied effects of loading curves with identical area under the curves, on the time history of the resulting von Mises and maximum principal strains as measures of the brain tissue damage.
In the present paper, two finite element models are presented for simulation of the traumatic brain injuries due to the increased intracranial pressure caused by the sever translational accelerations and decelerations following the shock wave and interaction between the skull and brain soft tissue, taking into account the hyperelastic nature of the brain tissue. In this regard, a realistic skull-brain model based on the MRI is created in CATIA modeling software and employed for mesh generation on the basis of an optimized accuracy algorithm, in HYPERMESH finite element code. Influence of the contact and nonlinear characteristics of the brain tissue are considered in simulation of the relative motions in LS-DYNA code. Finally, time histories of the accelerations and the von Mises stresses are derived from ANSYS finite element analysis software. Results are compares with those extracted based on the traditional spherical shell model and discussed in detail. Furthermore, various hyperelastic constitutive models are considered and their results are compared with each other and with the available experimental results.

Geometric and Finite Element Modeling
The skull (cranium) as a container, serves partly as a bony capsule for the brain and the central nervous system. It is composed of a brain case (neurocranium) and a facial skeleton (splanchnocranium, viscerocranium), both sharing the base of the skull, which descends obliquely backward (Fig. 1). All the bones that form the skull are connected together by some homogenous tissues and don't have any relative motions. The finite element model of the human's skull is constructed by considering the geometric complexities that are shown in Fig.3. Two models are used for the cranium cavity, based on the brain model: i. Spherical cavity, by simplifying the real model for easier performing of the dynamic analyses, as the models employed by Yue [4,5] and Chen and Ostoja-Starzewski [15]. ii.
Realistic MRI-based model. The cerebrum is the highest and biggest region of the brain that covers majority of the portions of the brain. Almost all of the skull space is filled by cerebrum. The cerebrum is composed of two main parts, namely, right cerebral hemisphere and left cerebral hemisphere, and several bilateral gray nuclei (basal ganglia). These two parts are separated at all sections but in some inside regions they have little connections by some white fibers.
In brain sections in the horizontal or coronal planes, the cerebral hemisphere can be seen to consist of an outer gray cortical layer; the cerebral cortex with about 2-5 mm thick and mainly composed of cell bodies and an inner white core composed mainly of nerve fibers (myelinated axons). On the other hand, the brain is composed of other parts such as: brainstem, cerebellum, temporal lobe, occipital lobe, parietal lobe, central fissure, frontal lobe and lateral sulcus. The cerebrospinal fluid (CSF) is a clear fluid, with a volume of about 130-150 ml, which fills the whole subarachnoid space. It acts as a protective fluid cushion around the brain and the spinal cord, damping shocks caused by blow or falls of the central of the spinal cord.
Mid-section of the 3D finite element brain model Universal Journal of Biomedical Engineering 4(2): [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]2016 13 associated with the first skull model is shown in Fig. 2. Yue [4,5] have assumed that the external diameter of the cranial cavity is about 200 mm; a diameter that corresponds to the circumscribing sphere of the cavity. Chen & Ostoja-Starzewski [15] have chosen this diameter as 140 mm that may be equal to the inscribing sphere. As mentioned before, similar models that have been employed by other researchers have proven to lead to somewhat satisfactory results for simulation of the brain trauma [4,5,15]. The realistic brain model of the present research is constructed based on the magnetic resonance images (MRI) of the brain and cranium which some of them are shown in Fig. 3 [26]. The whole scan set consists of 55 parallel sections. The three-dimensional geometric and finite element models were reconstructed based on the axial magnetic resonance images available in ''The Whole Brain Atlas" of Harvard Medical School [27]. In the present research, the geometric model is constructed in CATIA software whereas the finite element model (Fig. 4) is created in HYPERMESH finite element code to obtain elements with enhanced (optimized) topologies. The mesh includes the: (1) cerebrospinal fluid (CSF) in the form of a 3-mm thick layer, (2) gray matter, (3) white matter, (4) cerebellum, (5) corpus callosum, (6) telencephalic nuclei, (7) brain stem, and (8) ventricles. The whole assembly is illustrated in Fig. 5.
In both models 10-node second-order tetrahedral composite elements whose material properties are identical to that of the Wayne State Brain Injury Model [26], are employed. In both models, the 3mm gap between the inner surface of the cranium and the outer surface of the brain is filled by the cerebrospinal fluid (CSF). Due to complexity of the geometries, the general contact constraint is defined between these three media and the dynamic simulation is performed in LS-DYNA code.

Description of the Material Properties
No medical model may be found that can actually reproduce the normal characteristics or may be regarded as a realistic model for in vivo modeling and molecular characterization. All models are approximate; especially those relevant to the cellular evolution, aging processes, and degradation or changes in microscopic and even macroscopic properties of the cells or tissues. The most that can be expected from any model is supplying a useful approximation to reality. In this regard, the following two real difficulties in proposing a model for the in vivo functions of cell cycle, e.g., the cerebral tissues, may be emphasized: • Inside the skull cavity, there are many compartments separated by tissue of different elasticity, capacity, compliance, and resistance. These qualities should be taken into account in detail, since they substantially affect the changes that may occur during the TCE-s. • The patient's age in a specific way changes the characteristics of the intracranial cavity, such as cerebral mass, the content of CSF, and tissues that divide the compartments.
Hence, there is no model that can cover all variants, but it is advantageous to discuss procedures that may lead to development of an ideal experimental model. In this regard, the authors have previously attempted to discuss the problem in detail and published papers on techniques for characterization of the human tissues. Since the brain tissue is extremely soft, the traditional techniques cannot be employed for determinations of the material characteristics. The authors have already proposed and implemented special techniques to undertaken the mentioned task; e.g., they used the micro-indentation technique for the periodontal ligaments, using indenters with proper configurations [28]. However, due to the nature of the brain tissue, special indenters have to be used (the flat or semi-spherical indenters may be more efficient). Since our research is in progress, more developments may be reported in future.
Elastic behaviors of the human skull and the cerebrospinal fluid are considered to be linear and the relevant material properties are listed in Table 1. Studies performed on deformations of the brain tissue under impact loading confirm that this tissue mainly exhibits a hyperelastic behavior [2,[16][17][18][19][20][21][22][23][24]. Hyperelasticity refers to a constitutive response that is derivable from an elastic potential function W (e) (strain energy density function) and is typically used for nearly incompressible materials that experience nonlinear large elastic deformation (i.e., materials with nonlinear one-dimensional stress-strain curves in the elastic region) such as rubber and some biological materials. When the strain energy density function per unit undeformed volume is defined, one may determine the stress expressions from: where ij S , ij E , and ij C are components of the second Piola-Kirchhoff stress tensor, components of the Lagrangian strain tensor, and components of the right Cauchy-Green deformation tensor, respectively. Several forms of the strain-energy potential have been proposed for simulation of the incompressible or nearly incompressible hyperelastic materials, among them: Arruda-Boyce, Blatz-Ko, Gent, Mooney-Rivlin, neo-Hookean, Ogden potential, polynomial, and Yeoh models.
Evaluation of efficiency and accuracy of these models for the behavior and injury simulation of the skull-brain assembly is an important issue. In this regard, a comparative study is accomplished in the present research among results of four distinct hyperelasticity models and the previously reported experimental results. Mathematical bases of the mentioned four models are presented briefly in the following subsections.

Polynomial Model
The polynomial form of the strain energy potential is: where ij C are material constants that have the unit of force per unit area. 1

I and 2
I are the first and second invariants of the strain tensor, J is the elastic volume strain, and k D is the material incompressibility parameter that tends to zero for incompressible tissues. The neo-Hookean model can be obtained by setting N=1 and C 01 =0. We adopted the second-order polynomial strain energy density function in the present study where N=2 (the five-parameter Mooney-Rivlin model) and determined the required five parameters through solving the corresponding equations. In this model, the initial shear and bulk moduli are defined as:

Yeoh Model
The Yeoh [29] model is a reduced and revised form of the polynomial model (for N=3) wherein the strain energy potential depends on the first strain invariant only: The employed notations are identical to those of the previous model. The required equations may be solved to determine the three parameters of the model. In this model, the initial shear and bulk moduli are defined as: 10 1

Arruda-Boyce Model
Arruda-Boyce model just depends on 1 , the first invariant of the strain tensor. The strain energy density function of Arruda-Boyce model is: 3 11 1050 = C , 4 19 7000 = C , 5 519 673750 = C and l is the locking stretch (limiting network stretch) and is the initial shear modulus [30].
l and  are measureable parameters and: As the parameter l goes to infinity, the model reduces to the neo-Hookean model.

Ogden Model
The Ogden form of the strain-energy potential density is dependent on the principal stretches of left-Cauchy strain tensor: a a  l l l a (8) where i a are non-dimensional exponents and 1 2 3 , , l l l are parameters of the deviatoric stretches. As the polynomial model, a value of N > 3 leads to difficulties and is not recommended. The first-order Ogden model reduces the unknown variables to two parameters where this two parameters are the initial shear modulus  and 1 a .
Generally, theinitial shear and bulk moduli may be found from: For N=1 and 1 2 = a , Ogden potential is equivalent to the neo-Hookean potential. For N=2, 1 2 = a and 2 2 = -a , Ogden potential can be converted to the 2 parameter Mooney-Rivlin model.

Comparing Origins, Advantages, and Disadvantages of the Popular Hyperelastic Models
Each hyperelasticity model has been proposed based on some assumptions that are valid for specific conditions. Therefore, it may be adequate for specific materials or deformation ranges. For instance, Mihia et al. [31] claimed that the Fung and Gent models, which have typically been used to model soft tissues, are inadequate for the modelling of the brain or fat under combined stretch and shear, and so are the classical neo-Hookean and Mooney-Rivlin models used for elastomers. However, they found that a subclass of Ogden hyperelastic models is in excellent agreement with the experiments.
In particular, the shear modulus of normal brain can be increased nearly four times by compressive stresses. During experimental tests, tissue samples may exhibit a predominantly isotropic behavior and their volume may remain virtually constant [46]. An efficient hyperelastic model should have the following four main properties [47]:  It should have the ability to exactly reproduce the whole 'S' shaped stress-strain curves.  If the model operates sufficiently in the uniaxial tension, it must also be exact in simple shear or in a biaxial tension.  The number of fitting material parameters should be small, in order to decrease the number of the required experimental tests.  The mathematical implementation should be as simple as possible.
Based on the presented guidelines for adopting the proper brain hyperelasticty model, origins, advantages, and disadvantages of the most common hyperelastic models are compared in Table 2. -Suitable for small strains.
-It can provide rough predictions only.
-It is based on explanation of physics of the molecular chains and statistical methods.
-The strain energy is sum of the strain energies of chains oriented randomly in space [37].
-The complete Arruda-Boyce model is more accurate than the Neo-Hookean and Yeoh models.
-Suitable for small strains only.
-At high strains, the predicted tangent stiffness tends to infinity at finite strains.
-Generally less accurate than the polynomial and Ogden models.  [40,41]. -Ogden series, may give good fits in one deformation mode (namely uniaxial tension, uniaxial compression or pure shear) and very unrealistic results in another mode. -Its six material parameters necessitate a large experimental database to be fitted. -It can lead to errors when used out of the deformation range in which its parameters are identified.

Boundary Conditions
In the present research, possibility of the brain injury occurrence is investigated through simulation of an impact (collision) with a rigid plate. The initial velocity of the skull is assumed to be 7.5 m/s (27 km/h) and the front region of the skull hits the rigid plate. Value of the initial velocity is chosen identical to that of an experiment performed by Nahum et al. [48] and its results are available in literature.

The Available Criteria for the Traumatic Head Impact Injury Tolerances
Over the past decades, the biomechanical researchers provided criteria for tissues against impact which led to safety standards. While the Abbreviated Injury Scale (AIS) is used as a general qualitative measure, some quantitate criteria have been proposed for the head impact injury tolerances. The experimentally derived Head Injury Criterion (HIC) that has mainly been proposed for road traffic accidents may be interpreted in terms of the heat acceleration as follows [44]: a t dt t t t t (10) where a(t) is the head acceleration in g and t is time. The HIC has to be less than 1000 when a 36ms time duration is considered. Usually, t 1 and t 2 are so determined that the HIC is maximized by using a 15 ms limit. Studies accomplished by Prasad and Mertz [45] show that an HIC of 1000 is equivalent to an 18% probability of a severe (AIS 4) head injury, a 55% probability of a serious (AIS 3) injury, and a 90% probability of a moderate (AIS 2) head injury, to the average adult.
Following recent research [46,47], the injury risk may also be interpreted as a function of the mean acceleration: (11) in which critical a is the lower limit of the mean acceleration a for which the risk equals 1.0. The empirical constants critical a and n are derived by regression of real world injury risk data. Evans [48] investigated the injury and fatality probability for belted and un-belted drivers according to numerous car crashes and derived the associated mean accelerations for each risk level, using Monte Carlo methods. However, since this criterion does not consider the impact duration time that is an important parameter, it is not employed here.
The statistical correlation studies between the calculated head loadings and the observed injuries have led to the identification of other injury criteria. Concerning the brain neurological lesions, Newman et al. [49] and Zhou et al. [50] suggested that the mild traumatic brain injuries occur at a brain von Mises stress of 20kPa. Baumgartner et al. [51] observed that a brain von Mises stress of 18 kPa generates moderate neurological lesions which become severe from 38 kPa.
Viano [52] correlated the concussion with the translational acceleration. He observed that the nominal tolerance for concussion was HIC = 250. The conventional measures of head injury risk were effective in assessing concussion risks. He stated that no concussion occurs with average accelerations of 68 ± 15 g and HIC of 143 ± 37, while concussion may occur with 94 ± 28 g acceleration and HIC= 345 ± 181.

Specifications of the Finite Element Models and the Employed Material Properties
In the present section, a comparative study is carried out for a better judgment about the usefulness of the hyprelasticity models mentioned in section 3, for behavior simulation of the skull-brain system. Such comparative study has not been performed so far. The imposed boundary conditions are defined in section 4. Material properties associated with the adopted hyperelasticity models are listed in Table 3. To compare results of the traditional simplified spherical brain model against results of the realistic model (as well as the experimental results), both finite element models are considered and two sets of results are derived. Quadratic pyramid elements are used to discretize the models. Convergence of the created finite element meshes is checked through comparing results of successive refinement of the meshes; so that finally, ignorable changes in the results may be noticed by further refinement of the mesh. The performed convergence study whose results are not included here to save space, has revealed that choosing about 20000 elements is vital for achieving convergent results.

Stress Time Histories Predicted by Various Geometric and Hyperelastic Models
The experimental results include time histories of pressures of the coup and countercoup regions (Fig. 6 [53]) reported by Nahum et al. [43] and the time history of acceleration of the center of gravity of the brain measured by Trosseille et al. [54]. Results of the polynomial, Yeoh, Arruda-Boyce, and Ogden hyperelasticity models are illustrated in Figs. 7 to 10, respectively and compared with the experimental results. Figs. 7 and 8 correspond to the spherical finite element model (Fig. 2) whereas Figs. 9 and 10 are associated with the realistic finite element model (Fig. 4).
Deviations of results of the peak pressures of the employed finite element models from those of the experimental results are given in Table 4, for various hyperelasticity and finite element models.
Results illustrated in Figs. 7 to 10 reveal that maximum magnitudes of the coup and countercoup pressures predicted by the finite element models associated with the polynomial constitutive material models are less than those of the experimental (real) results. However, based on results listed in Table 4, the discrepancies between the finite element and the experimental results are lowest for the polynomial method. Moreover, curvatures of the time history curves of the finite element and experimental results show its most conformability. So that effects of the vacuums occurred in the coup and counter coup regions in the cerebrospinal fluid are simulated more accurately. Therefore, accuracy of predictions of this model is more reliable.
On the other hand, since the polynomial hyperelasticity model has underestimated the coup and countercoup peak pressures, its results are not in the safe side of the clinical predictions. In both finite element models, the peaks are slightly postponed in comparison with the experimental results. As it may be deduced from Figs. 7 and 8, the superimposed influences of the fundamental and higher natural frequencies on the responses of the spherical finite element model are detectable. This phenomenon has occurred because mass and stiffness of the spherical model are greater than those of the real one and consequently, natural frequencies of the mentioned model are different from those of the real one.
Results of the finite element Yeoh hyperelasticity model (Figs. 7 to 10) are slightly greater than the real (experimental) results. Thus, they are in the safe side of the clinical predictions. On the other hand, accuracy of these results is located in the second rank (after the polynomial hyperelasticity model). The predicted behaviors for the vacuum regions are not reliable, especially for the spherical finite element model. The peaks appeared in this region of response of the spherical finite element model are not reasonable. The absolute peak pressures occur slightly before the real occurrence times. As the previous and other models, discrepancies of results of the spherical finite element model relative to the experimental results are higher in comparison with results of the realistic finite element model.
In the Arruda-Boyce model, times of occurrence of the peak pressures conform to a great extent to those of the experimental ones. But as results of Table 4 confirm, discrepancies between the finite element and the experimental results are significant. Magnitudes of the peak pressures predicted by the finite element models are remarkably higher than those of the experimental ones. A much oscillatory response (in comparison to the previous hyperelasticity models) is observed in times following the peak pressures. However, responses of this region of the time history curve are more consistent with the experimental results (in comparison with Yeoh results). Therefore, it may be concluded that this model does not simulate the hyperelasticity nature of the brain tissues accurately.
The last hyperelasticity model, i.e. Ogden model significantly overestimates the peak pressures of both the coup and countercoup regions of the brain. Results presented in Figs. 7 to 10 and results reported in Table 4 reveal that in spite of using this model by some researchers [2], results of this model are not reliable at all. Results of the Ogden finite element model show the greatest discrepancies with the experimental results. Furthermore, redundant oscillations are appeared in the predicted responses shown in Figs. 7 to 10.
The absolute peak of each response is the most important parameter for injury assessment of the brain tissue and the blood vessels that are in the neighborhood of the coup and countercoup regions. Therefore, the minor oscillations can be ignored in the responses. In the realistic finite element model, the areas of the coup and countercoup regions are smaller than those of the spherical model. Therefore, it is expected that the peak pressures of the realistic model be greater than those of the spherical finite element model. But the mass and consequently, the inertia forces (and the linear momentum) of the spherical finite element model are greater; a fact that leads to prediction of higher peak pressures by the spherical finite element model (Figs. 7 to 10).
Comparing results shown in Figs. 7 to 10 with the experimental results reported by Trosseille et al. [54][55][56], reveals that predictions of the polynomial model are even closer to those of the real ones (if the peak pressures are of concern). Figs. 7 to 10 reveal that the coup peak pressures are greater than those of the countercoup region. The distributions predicted for the coup pressure by the mentioned four hyperelasticity models are shown in Figs. 11 to 14, for both the spherical and the realistic finite element models.

Time Histories of the Resulting Translational Accelerations
One of the general discrepancy sources is the shortcoming in accurately modeling the inertial and volume characteristics. For this reason, time histories of the acceleration of the center of gravity of the brain are determined by LS-DYNA software based on the various 22 Brain Tissue Response Analysis Based on Several Hyperelastic Models, for Traumatic Brain Injury Assessment hyperelasticity models and plotted in Fig. 15, for both the spherical and realistic finite element models. Results of the realistic model extracted based on the polynomial hyperelasticity model show a better agreement with results of the time history reported by Trosselie et al. [54] regarding both the peak acceleration and the relevant time of this peak. Again, Ogden model has led to the worst results. Results of Yeoh model are still in the second rank.
With exception of the polynomial hyperelasticity model, acceleration results of the spherical finite element are close to those of the realistic model. Therefore, both acceleration of the center of gravity of the brain and the pressure results lead to an identical conclusion that the models may be ordered with respect to the accuracy as: the polynomial model, Yeoh model, Arruda-Boyce model, and Ogden model.

Discussion of the Results from the Clinical Point Of View
Various criteria were explained in Sec.5; among them, the HIC, brain von Mises stress, and brain peak acceleration criteria. The HIC values are copmuted here for various hyperelastic models of the spherical and realistic model, numerically, using Simpson's numerical integration method and given in Table 5. With one exception, the HIC values predicted by the realistic model are higher. However, these values are less than the nominal tolerance for concussion (345 ± 181) proposed by Viano [52]. Although the HIC criterion does not predict severe brain injuries, predictions of Ogden and Arruda-Boyce are close to the concussion onset state. As may be noted, the ploynomial hyperelastic model, predicts lower injuries. The acceleration criterion of Viano [52] for concussion occurrence (a=94 ± 28 g) holds for all of the models (Fig.  15). Therefore, concusion is possible according to this criterion.
Regarding the brain von Mises stress critrion developed and extended by Newman et al. [49], Zhou et al. [50], and Baumgartner et al. [51] for neurological lesions, results of Figs. 7 to 10 may be used. Although results of the polynomial model underestimate the von Mises stresses of the brain (and this model has to be used with care from the safety point of view), magnitudes of all the coup and countercoup stresses are above 38 kPa. Therefore, all of the considered finite element and hyperelastic models predicts severe neurological lesions based on the resulting brain von Mises stresses as well.

Conclusions
In the present research, a comparison is made among results of various (polynomial, Yeoh, Arruda-Boyce, and Ogden) hyperelasticity models that may be candidate for simulation of the traumatic brain injuries, for the first time. A simultaneous comparison is made between results of the traditional spherical finite element model and results of the realistic MRI-based finite element model reconstructed by the authors. CATIA geometry modeling, HYPERMESH and ANSYS finite element modeling and solvers, and LS-DYNA dynamic modeling code are employed in the present research. A pair of experimental results is employed to evaluate accuracy of the peak pressures and accelerations of center of gravity of the brain predicted by the mentioned four hyperelasticity models for the two finite element models. Comparing with both experimental results confirms that employing Arruda-Boyce or Ogden models may lead to inaccurate or even erroneous results. Comparisons made with the experimental results of accelerations of the center of gravity of the brain and pressure results of the coup and countercoup regions lead to an identical conclusion that the models may be ordered with respect to the accuracy as: the polynomial model, Yeoh model, Arruda-Boyce model, and Ogden model. The traumatic head impact injury criteria reveal that the polynomial model has to be used with care due to underestimation of the injuries level and the acceleration and that the pressure criteria predict severe injuries for the considered simulation data.