A detailed and personalizable head model with axons for brain injury prediction: Sensitivity of brain strain to inter-subject variability of the brain and white matter tract morphology

8 Finite element head (FE) models are important numerical tools to study head injuries and develop 9 protection systems. The generation of anatomically accurate and subject-specific head models with 10 conforming hexahedral meshes remains a significant challenge. The focus of this study is to present two 11 developmental work: First, an anatomically detailed FE head model with conforming hexahedral meshes 12 that has smooth interfaces between the brain and the cerebrospinal fluid, embedded with white matter 13 (WM) fiber tracts; Second, a morphing approach for subject-specific head model generation via a new 14 hierarchical image registration pipeline integrating Demons and Dramms deformable registration 15 algorithms. The performance of the head model is evaluated by comparing model predictions with 16 experimental data of brain-skull relative motion, brain strain, and intracranial pressure. To demonstrate 17 the applicability of the head model and the pipeline, six subject-specific head models of largely varying 18 intracranial volume and shape are generated, incorporated with subject-specific WM fiber tracts. DICE 19 similarity coefficients for cranial, brain mask, local brain regions, and lateral ventricles are calculated to 20 evaluate personalization accuracy, demonstrating the efficiency of the pipeline in generating detailed 21 subject-specific head models achieving satisfactory element quality without further mesh repairing. The 22 six head models are then subjected to the same concussive loading to study sensitivity of brain strain to 23 inter-subject variability of the brain and WM fiber morphology. The simulation results show significant 24 differences in maximum principal strain (MPS) and axonal strain (MAS) in local brain regions (one-way 25 ANOVA test, p<0.001), as well as their locations also vary among the subjects, demonstrating the 26 necessity to use subject-specific models for studying individual brain responses. The techniques 27 developed in this study may contribute to better evaluation of individual brain injury and development 28 of individualized head protection systems in the future. This study also contains general aspects the 29 research community may find useful: on the use of experimental brain strain close to or at injury level 30 for head model validation; the hierarchical image registration pipeline can be used to morph other head 31 models, such as smoothed-voxel models. 32

1 Introduction 1 Traumatic brain injury (TBI) is a leading cause of injury-related death and disability, with a devastating 2 impact on the patients and their families (Maas et al. 2017). TBI causes a substantial threat to global 3 public health and enormous economic burden for the society, with an estimated number of 50-60 million 4 new TBI cases occur annually (Feigin et al. 2013). TBI influences all age groups, including children, 5 adolescents, and the elderly, which can happen due to traffic accidents, sports injuries, and falls. Studies 6 also suggest that TBI might represent an important modifiable risk factor for epilepsy, stroke, and late-7 life neurodegenerative diseases such as dementia and Parkinson's disease (Maas et al. 2017). 8 Biomechanical studies, including experimental and computational studies, have long been carried out to 9 understand brain injury mechanisms due to external mechanical input (Meaney et al. 2014). In particular, 10 finite element (FE) head models have emerged as valuable numerical tools to study head injuries and aid 11 the development of protection systems ( diffusion tensor imaging (DTI), brain-skull interface improvement, as well as mesh refinement, have led 22 to updated versions compared to the original. Especially recent efforts on mesh refinement led to average 23 brain element size of about 1.8 mm in the WHIM V1.5 model (Zhao and Ji 2019; Zhao and Ji 2020), 24 while the element sizes in GHBMC (Mao et al. 2013a) and the refined THUMS (Atsumi et al. 2016) are 25 also on the order of millimeter, being 2 mm, and 1.2~5 mm, respectively. However, the brains in the 26 classical models are simplified by smoothing out sulci and gyri, accompanied by a homogenous layer of 27 outer cerebrospinal fluid (CSF). Further, the brain ventricles often lack anatomical details and also have 28 jagged interface connecting with neighboring brain elements; some models or earlier versions have no 29 ventricles. Mesh simplification as such is partially attributed to the challenges for current meshing 30 techniques, e.g., blocking technique (Mao et al. 2013b) to capture anatomical details, while is also a 31 reasonable trade-off for computational efficiency. However, lacking these anatomical details hinders a 32 model's capacity for studying certain localized injuries such as at sulci, gyri and surrounding ventricles 33 (further discussion found below). Nevertheless, the classical models with high computational efficiency 34 have played important roles in improving our understanding of TBIs; some have found wide applications 35 for improved vehicle safety and helmet design. Modeling techniques learned from the classical head 36 models also pave the way for future models with higher anatomical accuracy. 37 To address anatomical accuracy, voxel-based approach has been used to generate head models including 38 detailed sulci, gyri, and ventricles (Ho and Kleiven 2009; Chen and Ostoja-Starzewski 2010; Miller et al. 39 2016; Ghajari et al. 2017). Voxel-based approach by converting voxels to hexahedral elements directly 40 or with various smoothing algorithms is efficient and has been used widely for FE analysis of bone 41 structures. However, a known concern is the less accurate peak strain/stress predicted from such models, 42 especially on the surfaces due to jaggedness (Viceconti et al. 1998; Samani et al. 2001). Although the 43 jaggedness could be reduced with various smoothing algorithms, e.g., (Camacho et al. 1997; Boyd and 44 Müller 2006), with a larger smoothing factor, which, however, is at the expense of decreased element 45 quality. Similarly, brain strains predicted from voxel-based head models may also have accuracy issues 46 at jagged surfaces of outer CSF-brain and ventricle-brain interfaces, but the accuracy level is unknown 1 and yet to be studied. Nevertheless, careful choice of result analysis, e.g., evaluating overall regional 2 brain strains or strain distributions, allows such models to provide valuable insights attributed to its 3 anatomical accuracy compared with classical head models, such as high strains at sulci depth (Ho and  4 Kleiven 2009; Ghajari et al. 2017), in line with an earlier experimental study (Lauret et al. 2009). 5 Integrating neuroimaging with model-predicted brain strains has provided a possible association between 6 mechanical response and chronic traumatic encephalopathy (CTE) (Ghajari et al. 2017). However, when 7 brain strains at the jagged interfaces are of primary interest, models with conforming meshes capturing 8 sulci, gyri, and brain ventricles are preferred. The jagged interfaces also hinder a reliable implementation 9 of sliding or fluid-structure interaction (FSI). Lastly, falx and tentorium also need to be manually 10 generated (Ho and Kleiven 2009;Miller et al. 2016), affecting its subject-specific efficiency; some 11 smooth-voxel model chose not to include falx/tentorium (Chen and Ostoja-Starzewski 2010). 12 Another technique for efficient generation of subject-specific models is by mesh morphing (also called 13 warping). The concept has been used extensively in many biomechanics fields on different organs 14 (Couteau et  geometrical difference between the subject and baseline mesh is obtained. Next, the displacement field 20 is applied to morph the baseline mesh, resulting in a personalized mesh with updated nodal coordinates 21 while remaining element connections. In general, the computed displacement field should comply with 22 continuum mechanics conditions on motion, requiring diffeomorphic, non-folding, and one-to-one 23 correspondence to avoid excessive element distortions (Bucki et al. 2010). Otherwise, without such 24 reasonable element quality, not only an FE analysis is prevented from being carried out, also numerical 25 accuracy is influenced. Morphing an anatomically detailed head model that has refined mesh sizes poses 26 a higher requirement on smoothness (associated with Jacobian) of the computed displacement field, 27 meanwhile provides an opportunity attributing to FE model's direct correspondence with neuroimaging 28 and allows utilizing the advanced registration algorithms developed within the neuroimaging field. 29 Therefore, although mesh morphing is efficient, one major challenge for using it to generate detailed 30 subject-specific FE head model is how to design an image registration pipeline that leads to a 31 displacement field representing well the inter-subject difference of local brain structures, meanwhile not 32 cause excessive element distortions. 33 Note that anatomical accuracy for a head model is only one among the several factors influencing its 34 biofidelity; other factors include material properties, representation of interactions between various 35 intracranial components (e.g., brain-skull interface). Nevertheless, an anatomically accurate model 36 provides a prerequisite for capturing local strains in areas of interest. Generation of anatomically accurate 37 and subject-specific head models with conforming hexahedral meshes remains a significant challenge 38 based on above literature review. Though conforming tetrahedral meshes are relatively easier to generate, 39 it's not preferred in head models intended for studying TBIs due to known unfavorable characteristics, 40 such as overstiffening and volumetric locking especially with incompressible material using first order 41 tetrahedral element; though second order could alleviate but may lead to a larger computational cost than 42 hexahedral meshes (Samani et al. 2001). 43 Despite promising progress, detailed TBI mechanisms remain largely unknown, reflected by not able to 44 predict clinical symptoms, and individual-specific injury tolerances may explain to some extent (Rowson 45 et al. 2018). Subject-specific models for more detailed mechanics of brain injuries are needed. Yet, how 46 brain morphology and WM fiber tract morphology differences among individuals may influence brain 47 injuries remain unclear, although a previous study investigated influences of brain sizes by global scaling 1 (Kleiven &  and Dramms deformable registration leads to personalized (i.e., subject-specific) models with 18 satisfactory element quality without further mesh repairing. The uniqueness of the ADAPT head model 19 is the equipped pipeline that allows fast generation of detailed subject-specific models with large 20 variations in head size/shape as well as local brain regions and lateral ventricles with competitive 21 personalization accuracy. The research community may find the hierarchical image registration pipeline 22 useful to morph other head models as well, such as smoothed-voxel head models. 23 This study is organized as below: Firstly, the development and validation of the ADAPT head model is 24 presented. Secondly, the hierarchical image registration pipeline for personalization is described and its 25 capacity is exemplified by generating six subject-specific head models with largely varying intracranial 26 volumes (ICVs) and brain shapes. Personalization accuracy is quantified by DICE similarity coefficients. 27 Lastly, we use the six subject-specific head models to study the influences of brain size/shape on brain 28 strain response under the same concussive impact. We hypothesize that a large variation in brain strain 29 and location may exist among subjects and could be revealed by anatomically detailed subject-specific 30 models. 31 software Hexotic to generate all hexahedral elements using an Octree-algorithm (Maréchal 2009). The 40 head model includes the brain, skull (compact and diploe porous bone), meninges (pia, dura, falx, and 41 tentorium), CSF, and superior sagittal sinus (SSS) (Fig. 1). The brain is divided into primary structures 42 of cerebral gray matter (GM) (i.e., cerebral cortex), cerebral white matter (WM), corpus callosum (CC), 43 brain stem (BS), cerebellum GM and WM, thalamus, and hippocampus. The cerebrum is further divided 44 into frontal, frontal, parietal, temporal and occipital lobes; CC is divided into midsagittal CC, and side 45 CC and each into four parts including the genu (anterior portion), the body (middle portion), isthmus 1 and the splenium (posterior portion); CSF are divided into outer CSF and ventricular system including 2 lateral ventricles, 3 rd and 4 th ventricles connected by cerebral aqueduct. Continuous mesh is used 3 throughout the model, with all meshes node-connected from the brain, pia, CSF and dura to the inner 4 skull, including all interfaces between the outer CSF and the brain near sulci and gyri. Note as the same 5 material property is used for the entire brain, dividing it into subcomponents is mainly for post-6 processing purposes. Further, a maximum level of recursive partitioning on the initial octree cube is set 7 to eight in the software Hexotic during meshing to allow capturing the complex structures of sulci, gyri, 8 and ventricular system, resulting in smallest element size of about 0.5 mm at these areas, and transits to 9 larger-sized elements to 1 mm and largest about 2.5 mm at inner brain areas as shown in Fig. 2b. A 10 smooth element size transition is ensured by the balancing rule implemented in Hexotic, with details 11 found in Maréchal (2009). The total number of elements in the head model is 4.4 million hexahedral and 12 0.54 million quad elements. The minimum Jacobian in the brain is 0.45. All simulations are conducted 13 with LS-Dyna 971 R11 using an explicit dynamic solving method. 14 15 16 Fig. 1 The ADAPT head model with major components illustrated (upper), embedded with WM fiber tracts (lower 17 left), and with connecting ventricular system to the outer CSF (lower right).

18
The brain is modeled as hyper-viscoelastic material to account for large deformations of the tissue, with 19 additional linear viscoelastic terms to account for the rate dependence. Material properties presented by 20 Kleiven (2007)  spatial coordinates, and the FA and 1 st principal eigenvector for this voxel is linked to each FE element. 16 The resultant FA mapped at FE brain resolution is shown in Fig. 2b. The final extracted WM tracts in 17 the whole brain contain polylines aligned with the FE head model is shown in Fig. 2c (left), from which 18 the CC and BS fiber tracts are enlarged (Fig. 2c right)   The performance of the head model is evaluated by comparing experimental data of brain-skull relative 2 motion, brain strain, as well as intracranial pressure close to or at injury level. For all the selected 3 validation experiments mentioned above, the model is scaled to match the anthropometric measurement 4 of the cadaveric heads. To further evaluate whether the model could predict brain response under 5 noninjuries level in living subjects, brain-skull relative motion and brain strain are compared with the 6 experimental displacements and strains measured in a human volunteer using tagged MRI during mild 7 frontal impact presented in Feng et al. (2010). All details of the validation setup are presented in 8 Supplementary Material, with a brief description provided below. 9 For brain-skull relative motion validation, neutral density target (NDT) displacement curves from seven 10 representative cases from Hardy et al. (2007) are selected, including one sagittal (C288-T3), one 11 horizontal (C380-T2), and five coronal impacts (C380-T1, C380-T3, C380-T4, C380-T6, and C393-T3). 12 The recalculated cluster brain strains of these seven selected cases presented in Zhou et al. (2019b) are 13 used as experimental strain data to evaluate brain strain performance of the model. Cluster brain strains 14 from the above seven cases are chosen, with further motivation provided in Discussion. Intracranial 15 pressure response of the head model is compared with recordings from experiment No. 37 conducted by 16 Nahum et al. (1977). 17 CORA (CORrelation and Analysis, version 3.6.1) scores are calculated to assess the level of correlation 18 between a pair of time history curves using a sub-method included in CORA, i.e., the cross-correlation 19 method. CORA score (CS) reported in this study is calculated as (V+G+P)/3 in terms of shape (V), size 20 (G), and phase (P), meaning equal weights for the three parts. CSs range from 0 to 1 with 1 indicating a 21 perfect match. Note another sub-method included in CORA, i.e., the corridor method is excluded, and 22 recommended settings from (Giordano and Kleiven 2016) are adopted in this study with details provided 23 therein. 24

Hierarchical image registration pipeline for mesh morphing 25
The personalization approach for subject-specific head model generation is based on Demons and 26 Dramms deformable registrations ( Fig. 3a-i). First, the diffeomorphic Demons registration (Vercauteren 27 et al. 2009) implemented in the open-source software Slicer 3D is performed between the segmented 28 cranial masks of the baseline (corresponding to the baseline FE mesh) and the subject after being rigidly 29 aligned using a 6 degree-of-freedom rigid registration available in Slicer 3D. Further details for the 30 Demons registration steps can be found in a previous study (von Holst and Li 2013). Afterward, Dramms 31 registration algorithm (Ou et al. 2011) implemented as open-source code by the authors (Dramms 32 version 1.5.1, 2018) is performed on the skull stripped T1W images inherited from the Demons step. 33 The resultant displacement field from the two-step registrations is then applied to morph the baseline 34 mesh of the ADAPT model, obtaining a subject-specific model. The subject's own DTI is then mapped 35 to the personalize model using the same procedure described in Sec. 2.2 resulting in a subject-specific 36 model incorporated with the subject's own WM fiber tracts (Fig. 3j). 37 1 Fig. 3 The workflow of the proposed hierarchical image registration pipeline for subject-specific head model 2 generation by morphing demonstrated with the results from the smallest female. Baseline T1W image (i.e., the 3 T1W image corresponding to the baseline FE mesh) and the subject's T1W are segmented to obtain the cranial 4 mask (a, b), which are used as input for Demons registration from which displacement field #1 is obtained as 5 indicated by the arrows (c). Displacement field #1 is further applied to the baseline T1W image (d), which is then 6 skull stripped (e) and afterward together with the subject's skull stripped T1W (f) as input to Dramms registration 7 (g). The obtained displacement field #2 from Dramms registration is applied to the baseline T1W and obtain the 8 personalized T1W (i.e., the T1W image corresponding to the subject-specific mesh), which is compared with 9 subject's T1W to evaluate personalization accuracy. Finally, the two displacement fields add up to morph the 10 baseline mesh (i), obtaining the subject-specific head model, including both the mesh and WM fiber tracts (j).

Subject-specific models with axonal fibers 12
The capacity of the personalization approach is demonstrated by generating six models with largely 13 varying ICV. First, six subjects are identified by analyzing the ICVs from the WU-Minn Human 14 Connectome Project (WUM HCP) database (WU-Minn 2017). For all the six subjects, both T1W and 15 DWI images are all openly accessible. T1W image is used for subject-specific mesh generation, and 16 DWIs are mapped to the personalized model, with details presented below. 17 18 1  Out of the 1200 subjects from the WUM HCP database, 3T structural scans of T1W are available for  2 1113 subjects, and the data of ICV are already processed using the software FreeSurfer (WU-Minn 2017). 3 As listed in Table 2, six subjects covering a wide range of ICV are selected, including the smallest head 4 (turns out to be a female), the largest head (turns out to be a male), and four heads with ICVs in between. 5

WM fiber tracts extracted from DWIs 7
The DWIs of the six subjects with a resolution of 1.25 mm are processed to extract WM fibers. The 8 downloaded DWI dataset has already been preprocessed, including corrections for gradient non-linearity, 9 motion-correction, and eddy-currents (WU-Minn 2017), which are further processed to extract diffusion 10 tensors and WM fiber direction (i.e., 1 st eigenvector) in each voxel using FSL v6.0.2 DTIFit with a 11 weighted linear least squares option. The WM fiber directions are then mapped directly to the subject-12 specific mesh using the approach described in Sec.2.2, based on which axonal strains are calculated. 13 14 The baseline T1W image is morphed to a personalized T1W image for each subject (see Appendix 1) 15 using the procedure presented in Sec. 2.4. DICE is then calculated to quantify personalization accuracy, 16 i.e., how well the personalized T1W image (corresponding to the subject-specific mesh) reflects the 17 subject's T1W as the ground truth. DICE is a single metric commonly used in neuroimaging field 18 ( implies no overlap at all between both, whereas a DICE coefficient of 1 indicates perfect overlap. 21

Personalization accuracy evaluation
To calculate DICE, automated segmentation is performed using the software FreeSurfer with the default 22 brain segmentation pipeline (recon-all) for both the personalized and subjects' T1W images without 23 additional manual editing. Segmented labels are combined to the following components for DICE 24 evaluation: cranial mask, brain, local brain regions of cerebral GM & WM, CC, BS, hippocampus, 25 thalamus, and cerebellum, as exemplified with two subjects (Fig. 4). Note the segmented labels are only 26 used during DICE calculation. Thus, the quality of the automatic segmentation has no influence on the 27 subject-specific mesh development process.

Loading conditions and brain strain evaluation 9
The six subject-specific head models are loaded with the same impact kinetics measured in a collegiate 10 American football player resulting in loss of consciousness reported earlier (Hernandez et al. 2015). 11 Translational accelerations and rotational accelerations (Fig. 5) are imposed on the center of gravity 12 (C.G) of the head models. 13 Maximum of the 1 st principal Green-Lagrange strain (MPS), and maximum axonal strain (MAS) (i.e., 14 strains along the WM fiber direction as defined in Eqn.1) during the entire impact, as well as the locations 15 of both metrics, are analyzed and compared between the six subjects.   8 The CORA scores (CSs) for the ADAPT head model on brain motion are presented in Table 3 in  9 comparison with head models previously developed at the same research group, including the original 10 KTH head model (Kleiven 2007) and its updated version with FSI for brain-skull interface (Zhou et al. 11 2019a) interface (referred to as KTH-FSI model), as well as KTH detailed head model (Zhou et al. 12 2019b). Note that, to make the validation performance of different models amenable a direct comparison, 13 the CSs for previous models reported here are either newly calculated or recalculated using exactly the 14 same approach and CORA settings described in this study. Notes: 1 For each case, one CORA score is reported, which is the mean of CSs for all the evaluted NDTs as plotted in Appendix 2. 2 Simulations run for all the seven cases using the KTH head model described in Kleiven (2007) and CSs newly calculated. 3 CSs recalculated based on curves reported in the original study where a different CORA calculation was used. N/A: cases not run in the original study; 4 CSs for the first six cases newly calculated based on curves reported in the original study as CSs were not reported. For case C393-T3, simulation is run for this study with CS newly calculated.

Brain strain 4
The mean CSs for the ADAPT head model on principal and shear strain for the seven evaluated clusters 5 are 0.763 and 0.776, respectively ( The predicted brain strain curves from the ADAPT model are compared with the experimental brain 10 strain data presented in Zhou et al. (2019b) (Fig. 6). Despite a large difference in the peak between the 11 model predicted and experimental data (except for C380-T1 C1 which is closer), the shape and phase 12 show good match, reflected by the high values of V and P; close to 1 in some cases (C380-T3 C1, C393-13 T3 C1 principal strain) ( Table 4). Further, the simulated brain strains are consistently lower than the 14 experimental strain in all the seven evaluated clusters except for C380-T1 C1. More clusters are to be 15 studied to see if the same trend holds true for all the 15 cluster brain strains presented in Zhou et al. 16 (2019b). Further discussion on brain strain validation performance and the implications are found in 17 Discussion. 18 The mean CSs for brain strain are higher than brain-skull relative motion, attributed to high scores in 19 shape and phase compensating the low values in size (G). Using the brain-strain based CSs, the ADAPT 20 head model would be rated as "good" according to the same rating scale used earlier (Zhao and Ji 2020). 21

Intracranial pressure and in-vivo strain comparison 7
The CSs for the ADAPT head model on intracranial pressure are presented in to highlight the need and serves as a basis for future investigation. 20

Personalization accuracy evaluation of DICE for the six subject-specific models 1
The boxplot of the DICE values of the cranial, the brain, and local brain regions are presented in Fig. 7  2 (the DICE values are listed in Appendix 1), in general showing quite good results even for CC with 3 large variations between the baseline and subjects. Especially an average DICE of 0.96, 0.93, and 0.72 4 are achieved for the cranial mask, the brain, and hippocampus, comparable or even higher than some 5 algorithms used in neuroimaging field (Ou et al. 2014) for capturing inter-subject differences. DICE 6 values for local brain regions, as well as lateral ventricles, are all above 0.6, indicating the internal brain 7 structures of the subject-specific head model reflect the subject to an acceptable level. Cranial mask 8 difference defined as (A-B)/B is also calculated and shows an average of 1.6% for the six subjects.

Subject-specific head model mesh quality 14
In the six subject-specific head models generated, most brain elements (95.5% ± 1.2% on average for 15 the six subjects) have a Jacobian over 0.5, and almost all elements (99.9%±0.1%) have a Jacobian over 16 0.45. The minimum Jacobian in the six head models is all above 0.2. In this study, the mesh quality is 17 considered to be satisfactory when at least 95% of the elements have a Jacobian over 0.5. A summary of 18 brain mesh element qualities is listed in Table A2 (Appendix 1). 19

Magnitude and location of MPS in the cerebral cortex 20
The time-history curve of MPS in the cerebral cortex is extracted for each subject (Fig. 8a). Interestingly, 21 the smallest female shows MPS occurring at 36 ms, slightly different from other subjects except for the 22 largest male occurring at 56 ms (Fig. 8a). Note the delay in the peak between both curves is mainly 23 caused by the location difference where maximum strain occurs during the entire impact. The MPS in 24 the cerebral cortex for all subjects is located at the sulci regions, as exemplified with results from two 25 subjects (Fig. 8b). Additional animations of the brain strain response during the entire impact for both 26 subjects are provided as Supplementary Videos. 27 The results show, in general, smaller heads tend to have a lower brain strain under the same impact  The locations of MPS for the six subjects are notably different, shifting  2 between frontal, parietal, and temporal lobes (Fig. 8c).

Brain regional response analysis: MPS and MAS 9
The time-history curves of the 1 st G-L principal strain from all elements of cerebral WM, CC, BS, 10 hippocampus, and thalamus are evaluated and exemplified with the results from the smallest and largest 11 heads (Fig. 9). The color shaded shape is formed by strain-time history curves from all elements, of 12 which the curve from the element with maximum strain is plotted with black color line for each brain 13 region. Of interest note that MPS typically occurs at a similar time between both subjects in most brain 14 regions (i.e., similar phase in black color line between two models), except for CC and thalamus (Fig. 9, 15 2 nd and 5 th row), shifting to a later time in the largest head compared with the smallest head. The 16 difference in the shaded shape indicate element-wise different MPS response between the two models, 17 also indicate a different location of MPS. For example, the shade shape for CC and thalamus are notably 18 different between the two subjects, indicating the largest strain likely occur at a different element 19 (location) between the two models at CC and thalamus (similarly as observed for strain at cerebral cortex 20 as shown in Fig. 8b). The results for axonal strain are plotted in Fig.10 with time-history curves of axonal strain in all elements 6 of cerebral WM, CC, BS plotted. Note axonal strain for thalamus and hippocampus not evaluated (grey 7 matter region with less anisotropy), thus, is not plotted. Similarly, the color shaded shape varies between 8 the two subjects among brain regions, indicating element-wise different MAS response between the two 9 models. The MAS in cerebral WM occurs much later in the smallest female than in the largest male ( Fig.  10 10, 1 st row), a different trend than observed for MPS with maximum strain occurs at similar time for 11 WM (Fig. 9, upper row).  The time-history curves of all elements for only two subjects are presented above, and for all other 17 subjects, only the maximum values (MPS and MAS) are presented in the boxplot (Fig. 11)

Brain shape influence 5
The 50 th percentile female and 5 th percentile male have very close ICV, and brain volumes are also 6 similar (5% difference). Thus, it would be interesting to analyze further to understand how brain strains 7 vary between subjects with similar ICV but different shapes. Simulation results from the two subjects 8 are compared, showing MPS in the cerebral cortex differs by 14.2% (Fig. 12a), occurring at a quite 9 different location as shown in Fig. 8c, despite a similar pattern in strain distribution (Fig. 12b). For CC, 10 not only the MPS value differs between the two, more importantly, the occurring time (Fig. 12c), and 11 consequently, the strain distribution pattern when MPS occurs (Fig. 12d). MPS is located at the mid-12 body of CC in the 5 th percentile male and located at the anterior in the 50 th percentile female. The results 13 provide evidence that for heads with similar ICV, the magnitude of MPS, location, and strain pattern can 14 vary significantly due to head shape difference.  This study presents an anatomically detailed and personalizable head model with WM tracts embedded 7 (the ADAPT head model), and a hierarchical image registration pipeline for subject-specific head model 8 generation by mesh morphing. The model is validated against experimental brain motion and brain strain 9 close to or at injury level, as well as intracranial pressure, showing overall CORA score comparable or 10 higher than earlier models. The developed pipeline allows efficient generation of anatomically detailed 11 subject-specific head models with satisfactory element quality. Subject-specific head models generated 12 using this approach are shown to capture well the subjects' head geometry for the six subjects of largely  13 varying ICVs, both on a global level (cranial mask and the brain), and local brain regions as well as 14 lateral ventricles. Brain strains of MPS and MAS show significant differences among the six subjects 15 due to head/brain size & WM morphology variability, motivating the necessity for using subject-specific 16 head models for evaluating brain injuries. 17 The significant difference in both MPS and MAS magnitude and locations among the subjects (Fig. 8 -2 Fig . 11) seems to suggest ICV as a dominant factor influencing brain strain. Further, heads with larger 3 ICV tend to have increased brain strains under the same loading, though not necessarily follow a 4 monotonic increasing trend. Notably, the two subjects with very similar ICV show quite different strain 5 patterns (Fig. 12), indicating head/brain shape may also be an important factor influencing brain strain 6 response. It will be interesting to investigate in the future to clarify whether size or shape is more 7 important than the other influencing brain mechanical response. For example, by using principal 8 component analysis (PCA) (Wold et al. 1987) with more brain images, which may also allow identifying 9 the characteristics of brain shape and WM morphologies that are most vulnerable to impact. 10 Percentile strain (e.g., 95 th percentile MPS) have been used in previous studies both for adults (Miller 11 et uses MPS (i.e., 100 th percentile MPS) to compare between subjects. Using MPS also allows identifying 14 the element with maximum strain and plotting its time-history curves, as presented in Fig. 8a, which is 15 from a brain element on the cortical surface. Visual check shows the strain in the plotted element has no 16 abrupt difference comparing with its neighboring elements at the cerebral cortex (i.e., sulci/gyri) 17 attributed to the conforming mesh, giving some confidence in the plotted curves, as well as the identified 18 varying locations among subjects (Fig. 8c). Fig. 9 with strain curves plotted from all elements at brain 19 regions also shows no curves jumping outside, brings further confidence on the reported MPS. Finally, 20 the time-history strain curves shown in Fig. 10b from one element in the CC, despite the non-smooth 21 interface between CC and neighboring brain elements due to the same brain material used, no abrupt 22 strain differences are observed in the plotted element comparing with its neighboring elements. 23 Influences of head/brain size on brain mechanical response have been studied in the past. In a 2D 24 numerical study by Prange et al. (1999), coronal rotational accelerations were applied to evaluate brain 25 strain response due to brain size differences between adults and children. Kleiven and von Holst (2002) 26 by globally scaling a 3D adult FE head model to six heads of varying dimensions, showed brain response 27 increased almost monotonically from the smallest to the largest head under a linear acceleration. 28 Similarly, as revealed in this current study, a larger ICV (relating to larger brain mass) tends to have a 29 larger strain under the same impact, which is also suggested by Holbourn's scaling principle (Ommaya 30 et al. 1967), and more recent work by Wu et al. (2020) and Panzer et al. (2014). Of interest would be to 31 investigate whether brain strains predicted from FE models follow or can be predicted by the 32 acceleration-mass scaling law, which indeed has been studied by Prange et al. (1999). Their results 33 demonstrated that the mass scaling relationship was not sufficient to produce brain strain distributions. 34 Similarly, the models in the current study incorporating sulci and gyri, poses even greater challenges for 35 such scaling laws to relate MPS with brain masses. The anatomically detailed subject-specific head 36 models representing individual brain structural differences appear to be critical for revealing the new 37 insights on brain size/shape influence ( Fig. 8 and Fig. 12 capacity for handling large shape differences by performing Demons registration as a first step with 10 binary cranial masks as input that allow obtaining a displacement field reflecting well overall cranial 11 shape (mean DICE of 0.96, Fig. 7). Dramms registration is performed in the 2 nd step on the skull stripped 12 T1W images inherited from the 1 st step. The focus of the 2 nd step registration is to align local brain 13 structures and CSF (outer CSF and ventricles) by using brain MRI information without need of 14 segmentation. The choice of Dramms is based on its promising performance due to its advancing in 15 hierarchical attribute matching mutual-saliency mechanism (Ou et al. 2011). It will be interesting to 16 investigate if other nonlinear algorithms would achieve good performance for the detailed brain or not, 17 such as B-spline and Burr's elastic previously used for morphing head models with smooth brains ( Dramms in the 2 nd step may result in better performance or not. Nevertheless, the proposed hierarchical 21 registration pipeline leads to competitive performance in generating detailed subject-specific head 22 models with satisfactory element quality without the need for further element repairing, meanwhile 23 achieves DICE values comparable to or higher than that in the neuroimaging field (Ou et al. 2014). With 24 the established workflow, new advanced registration algorithms developed in the neuroimaging field can 25 be readily implemented to generate subject-specific models reflecting even better the intersubject-26 variability for both brain and WM fiber morphologies. The hierarchical image registration pipeline is 27 not only applicable for this current ADAPT head model but also can be used to morph other head models 28 as baseline, such as smoothed voxel head models. Further, due to its capacity for handling highly 29 nonlinear warping, the pipeline can also be used to generate models with pathologies with brain structural 30 changes such as decompressive craniotomy when the brain is expanding outside the skull (von Holst 31 et al. 2012). 32

Brain-skull relative motion validation performance comparing with previous models 33
Regarding model validation performance on brain-skull relative motion, the ADAPT head model shows 34 consistently higher CSs comparing with the original KTH head model, while has comparable CSs with 35 the KTH-FSI model and the KTH detailed head model. The seven cases for brain-skull relative motion 36 validation selected here are because they consistently have a duration longer than 40 ms. Besides, only 37 reliable NDTs are included in CS calculation, justifying the same weight factor used for all NDTs (i.e., 38 CORA score for each case is a mean of NDT curves). CSs from previous models are either newly 39 calculated or recalculated in this study to ensure all values are directly comparable among models. In 40 fact, CSs for the original KTH head model for C288-T3, C380-T4, C380-T6, C393-T3 have been  41 presented earlier, showing comparable CSs with THUMS, GHBMC, see Table 8 in the study by 42 Giordano and Kleiven (2016). Note despite this study followed the same CORA calculation method 43 (with corridor method excluded) and used recommended global settings same as proposed by Giordano 44 and Kleiven (2016), the same weighting factor is used for all the evaluated NDTs instead of using the 45 proposed weights for each NDTs as listed in Table 5 (Giordano and Kleiven 2016). The reason is that 46 seems no consensus has been reached among the research community on the proposed weighting factors 47 and using equal weight for all NDTs is justified by only include reliable NDTs as plotted to allow easier 1 comparison for future studies. Nevertheless, the CSs for the overlapping cases are close between the two 2 studies, e.g., for C380-T4 5.26 shown in Table 8 in Giordano and Kleiven (2016), close to 0.551 reported 3 here (note differing by a factor of 10 in the calculation). 4

Experimental brain strain for head model validation 5
Regarding experimental brain strain data used for model validation, the cluster brain strain presented in 6 Zhou et al. (2019b) is used, which is recalculated based on the original brain-skull relative motion 7 experimental data from Hardy et al. (2007) using a tetra approach instead of a triad approach used in the 8 original study ). Though it's well recognized (Zou et al. 2007; Zhao and Ji 2020) and 9 recently has been extensively verified (Zhou et al. 2019b;Zhou et al. 2018) that a model validated against 10 brain-skull relative motion may not necessarily guarantee its strain prediction accuracy. Therefore, it's 11 suggested that a head model with an intended use for strain prediction should be validated against 12 experimental brain strain data. However, despite the availability of the strain data presented in Hardy et 13 al. (2007) along with the brain motion data, it's seldomly used for head model validation in contrast to 14 motion data being widely used. This may be partially attributed to concerns on the quality of strain data, 15 especially the initially reported ~ 2-5% peak strains are rather low for an injurious impact (see a detailed 16 discussion by Zhao and Ji (2020) and references therein). To address the quality concerns on strain data 17 originally presented in Hardy et al. (2007), a two-step effort has been undertaken recently. As a first step, 18 Zhou et al. (2018) reanalyzed and updated the brain strain data using the same triad approach and 19 developed three criteria to assess the eligibility of the NDT clusters suitable for strain calculation. As a 20 second step, a tetra approach is used for brain strain calculation (Zhou et al. 2019b), reflecting better the 21 3D experimental brain deformation than the triad approach. Note the tetra approach is just one of another 22 way to estimate brain strain from NDT motion, an alternative approach has been proposed also, e.g., a 23 generalized marker-based strain sampling approach to estimate and compare regional strains (Zhao and  24 Ji 2020). Indeed, brain strains calculated by the tetra approach are much larger than by the triad approach, 25 indicating the original ), also the reanalyzed triad strains (Zhou et al. 2018) largely 26 underestimated the experimental brain strain, thus is recommended not to be used for head model 27 validation. Note that for nearly incompressible material as the brain, shear strains should be close to 28 principal strains, indeed as the ADAPT model predicted for all the seven clusters (see Fig. 6 simulated 29 shear strain and principal strain). However, unexpected large differences are observed between 30 experimental principal strain and experimental shear strain (Fig. 6). However according to a study by 31 Bradshaw (2003), it is quantified that for brain tissue, the principal strain should fall within the range of 32 2/3 to 4/3 of the maximum shear strain based on the theoretical interpretation. While when reanalyze the 33 experimental strain data from Zhou et al. (2019b), of all the 15 clusters, the experimental principal strain 34 consistently falls within the shear strain band, correlating with the theoretical finding by Bradshaw 35 (2003). Such correlation seems to indicate an acceptable quality of the 15 cluster brain strain dataset 36 including the 7 used here for validation. Therefore, the recalculated cluster brain strains using the tetra 37 approach have been verified and justified thoroughly (Zhou et al. 2019b), thus can be taken as 38 experimental brain strain for validating brain strain response of FE head models. 39

Brain strain validation performance 40
Seven cluster strains are chosen (out of 15 in total) for model validation as they consistently have 41 duration over 40 ms. The shape and phase show a good match between ADAPT model simulated and 42 experimental brain strain, but a large discrepancy in magnitude is observed in most clusters (Fig. 6). 43 Further, the simulated brain strains are consistently lower than the experimental strain in all the evaluated 44 clusters except for C380-T1 C1. It will be interesting to study more clusters, especially, experimental 45 strains for C241-T5 C2, C241-T6 C2, C393-T2 C1 have much lower peak strains (all <0.1) (Zhou et al. 46 2019b). It's expected that for some of these clusters, simulated brain strain may be larger than 1 experimental strain, which, however, is yet to be performed. This also reminds one important difference 2 comparing with NDT motion data (many NDT curves even for one case, e.g. 14 NDTs and each has 3 XYZ displacements), that choosing only several or certain cluster strains may lead to a biased impression 4 that the simulated brain strain too low compared with experimental data; this further reminds that caution 5 should be taken when tuning material parameters to satisfy brain strain performance when only use a 6 limited number of clusters (as shown in Fig. 6 for the seven clusters selected in this study consistently 7 shows lower strain in the model, but may not be the case when more cluster are evaluated as discussed 8 above). 9 Therefore, the large peak discrepancy between ADAPT and experimental cluster strain is suggested not 10 to be seen as a concern for the model performance, rather it opens a question not specific to this model: 11 How to use the experimental cluster strain for validation for the models that chose to use this strain data. 12 Especially, how to evaluate model performance when only a few cluster curves are available versus NDT 13 motion curves where are many (note CSs for brain strain is actually higher than brain motion despite 14 large magnitude discrepancy, seems to remind CSs may not be a proper index). How to weight between 15 NDT motion and strain to reach an overall biofidelity ranking for a model. Further, how to extract brain 16 strain for a model with coarse mesh. The recalculated experimental cluster brain strain data is a result of 17 the two-step efforts (Zhou et al. 2018, Zhou et al. 2019b) attempted to best utilize the NDT motions close 18 to or at injury level originally measured by Hardy et al. (2007). The recalculated experimental cluster 19 brain strain data provides a possibility to evaluate head models' brain strain prediction capacity; however, 20 there are challenges, as mentioned above, that need integrated effort among the research community to 21 allow its proper use for head model validation. it may not be appropriate to adopt material properties of the brain obtained from a coarse mesh to a 28 model with a much finer mesh, and vice versa (Zhao and Ji 2019). This is indeed a valid concern 29 especially if brain parameters were obtained by adjusting/optimization to satisfy model validation, 30 resulting model-specific material properties that may not be directly translatable to other models, 31 especially with largely different brain mesh sizes. Keeping this in mind, it seems questionable to adopt 32 brain material properties from the KTH head model to ADAPT. However, note that the brain material 33 properties presented in (Kleiven 2007)  loadings. Indeed, ADAPT model, in general, has a larger brain-skull relative motion for most NDTs (see 43 Fig. A3 to Fig. A5), and has higher CSs for the seven cases than the KTH head model (Table 3). Further, 44 the detailed brain morphology (sulci/gyri) included in the current model may potentially alter brain tissue 45 response, considering the CSF penetrates and pia mater holding the brain at the sulci/gyri comparing 46 with models with smooth brain surface. Nevertheless, a more systematic investigation is required to 47 study the potential prediction differences between ADAPT and the original KTH head model. versus anisotropic) may influence brain strain in an anatomically detailed head model is yet to be studied. 18 19 There are some limitations in the ADAPT model to be noted. The ADAPT model has conforming mesh 20 at all interfaces between the entire brain and CSF (i.e., outer CSF-brain/dura, ventricle-brain). Still, non-21 smooth boundaries exist at brain subregions, including GM and WM interface. The subregions of the 22 brain, and diploe porous skull bone are grouped according to the spatial correspondence with the 23 segmented images via an automatic script, which seems often similarly practiced in classical head 24 models as well, e.g., WM meshes manually picked according to geometry data in the GHBMC (Mao et  25 al. 2013), and the same for the KTH head model (Kleiven 2007). Despite the non-smoothness, strain 26 concentration is often not a concern, especially if using the same material properties for the entire brain 27 as done in this study. Further, a continuous mesh is used throughout the ADAPT model, which can be 28 further improved by including FSI at the brain-ventricle interface as done earlier (Zhou et al. 2020), and 29 further to implement FSI at the brain-skull interface (Zhou et al. 2019a) in the future, despite a technical 30 challenge to implement FSI on the complex sulci and gyri than a smooth brain model. 31

Conclusion 32
This study presents the development of an anatomically detailed head model with conforming hexahedral 33 mesh (the ADAPT head model) equipped with a hierarchical image registration pipeline for efficient 34 generation of subject-specific models by mesh morphing. The model is validated against brain-skull 35 relative, brain strain, and intracranial pressure, showing comparable performance with previous models. 36 The six-subject specific head models generated using the ADAPT model and the pipeline demonstrate 37 the capacity of the ADAPT model and the pipeline for fast generation of anatomically detailed subject-38 specific head models with largely varying brain sizes/shapes with competitive personalization accuracy 39 on capturing individual's brain structures. The simulation results show significant differences in brain 40 strain and axonal strain, motivating the necessity for using subject-specific head models for evaluating 41 brain injuries. The ADAPT model due to its uniqueness of the complete ventricular systems, including 42 3 rd and 4 th ventricles in connection with outer CSF via aqueduct, together with conforming meshed sulci 43 and gyri and subject-specific WM fibers, could potentially provide new insights into TBI mechanisms. 44 The verified performance of the ADAPT head model equipped with the personalization approach 45 addresses the challenges in subject-specific FE model generation, opening an opportunity for studying 46 for studying personalized brain responses, as well as developing personalized head protection systems. 1 The research community may find the hierarchical registration pipeline useful to morph other 2 anatomically detailed head models, such as smoothed-voxel head models. 3

Acknowledgments 4
The authors would like to acknowledge the earlier version of the Python code from Chiara Giordano,5 which was adapted and used in this study. Our special thanks go to Dr. Yangming Ou and Dr. Steve 6 Pieper, who generously spent their time advising Dramms registration, and to Dr. Loïc Maréchal who 7 provided great support for the usage of the software Hexotic. The authors thank the three anonymous 8 reviewers and Prof. Peter J. Hunter for their stimulating comments and valuable suggestions that 9 substantially improved this paper. Brain tissue response of the 1 st principal strain during the simulated entire concussion impact for two 23 subjects: the smallest (female) and largest (male). 24

Appendix 1 Personalization accuracy and mesh quality of the personalized models by morphing 25
The baseline T1W image Fig. A1a) is morphed to a personalized T1W image for each subject (Fig. A1b) 26 using the procedure presented in Sec. 2.4. The personalized T1W image is paired with the subject-27 specific mesh. Thus a comparison between it and the subject's own T1W (Fig. A1c) (the ground truth) 28 allows evaluating how the subject-specific model reflects the subject's head geometry. The results show 29 not only the overall head shape, but also internal brain structures are well correlated (Fig. A1 b, c). The 30 personalized T1W images (Fig. A1b) and subject T1W (Fig. A1c) are then segmented for the cranial 31 mask, brain, and seven local brain regions illustrated in Fig. 4 in the main text for DICE calculation.   Appendix 2 Brain-skull relative motion and brain strain validation comparing with previous head 9 models 10 Brain-skull relative motion for three cases predicted from the ADAPT, and the original KTH head model 11 (Kleiven 2007) are compared for C288-T3 (Fig. A3), C380-T1 (Fig. A4), C380-T2 (Fig. A5), together 12 with experimental data from ). Note that C288-T3 NDT8 has a duration of 39 ms (<40 1 ms). Thus NDTs from 1 to 14 except NDT8 are included in the plot also in the CORA calculation. 2 Besides, for the six impacts delivered to subject C380 (i.e., C380-T1 to T6), the motion of NDT9 failed 3 to exhibit a motion trend toward its starting position (details found in (Zhou et al. 2018). Thus NDT1- 4 14 except NDT9 are included in the plot also in the CORA calculation. Finally, for C393-T3, NDT13, 5 therefore, NDT1-14 are included in the plot also in the CORA calculation. 6 CSs calculated for the KTH detailed head model (Zhou et al. 2019b) on brain strain are presented in 7 Table A3 in comparison with the ADAPT head model presented in Table 4 in the main text. 8

4 5 6
Appendix 3 Brain strains predicted from the six personalized models 1 The MPS and MAS for brain regions are presented in Table A4 &Table A5, corresponding to the 2 boxplot in Fig. 11 in the main text. 3 4 Table A4 MPS for brain regions predicted from the six subject-specific head models.