On the level of computational model of a human skull : A comparative study

In this study, different patient-specific computational models of the skull, which are often used in literature, were investigated, analysed and compared. The purpose of this study was to demonstrate the differences in computational model creation and results in case different computational models based on same computed tomography (CT) dataset are used. The selection of computational model directly influences the values of investigated parameters. The effort is to demonstrate, how the selection of the computational model influences the results of biomechanically relevant parameters. The comparison was based on total displacement of the skull and von Mises strain investigated around predefined paths around the skull. The strain values were evaluated according to criterion from literature. The results were obtained using finite element method. The values of the displacement of the skull were higher in case of considering cancellous bone tissue due to its poor material properties or heterogeneous material properties. The same situation occurred during the evaluation of strain. The values were higher in models which include cancellous bone tissue in the structure. c © 2018 University of West Bohemia. All rights reserved.


Introduction
Methodology that has become very popular in implantology-related simulations and production in last few years is the patient specific approach [15].This approach has been used for several years for skull implants as well.To solve problems from the cranial biomechanics field, computational modelling is used very often.Computational modelling allows the assessment of the mechanical interaction and determination of stress-strain states in a bone tissue-implant system.
A valuable analytical tool that offers a solution to the problems outlined above is computational modelling using numerical methods.In recent years, finite element method (FEM) has been introduced as a successful representative of computational modelling [23].It has been implemented in several commercial engineering analysis software, which are commonly used also for solving biomechanical problems.The computational models might be created at different levels depending on the purpose of the model and on the availability of input data.In literature, several types (levels) of computational models can be found.The simplest ones are created using only shell elements [22] or volume elements [20] assuming a simple homogeneous, isotropic, and linearly-elastic material representation.A more complex one contains also heterogeneous material properties [2] which differs due to bone tissue density in different locations.Using of a specific type of model is also related to hardware and time requirements for each problem.
The effects of mechanical interactions are linked to stress-strain states which might be, however significantly affected by the level of computational model.
Despite the recent progress in computer simulation possibilities, even the simplest models are still used in bone-biomechanics simulations.Using such models are often criticized just for their simplicity forgetting the famous aphorism by George Box: "All models are wrong, but some are useful."Cranial biomechanics is not an exception.To the best of authors' knowledge, there is no paper discussing differences among computational models of the same subject but constructed at various modelling levels.
The aim of this study is, therefore, to create various patient-specific computational models of human skull assuming different modelling levels which are most frequently found in literature and to compare outputs from such models.The emphasis was put on frequently used model types but also an approach still not common (such as assuming heterogeneous CT-based distribution of material properties in bone) in cranial biomechanics was considered as well.

Materials and methods
A CT dataset of a human skull (male, 23 years old) was used in this study.CT images were processed in the STL Model Creator software developed by authors [13,14] in Matlab (Matlab 2012, MathWorks, USA) to obtain 3D standard tessellation language (STL) file.The image voxel size was 0.4 mm in x and y coordinate and 0.7 mm in z coordinate.The STL file was processed in software Catia (CATIA V5, 2016, Dassault Systèmes, France).Calculations were done in software Ansys (Ansys 17.2, Ansys Inc., USA).

Model of geometry
In this study, only a part of skull was created (length ∼ 170 mm, width ∼ 140 mm, height ∼ 90 mm and volume ∼ 187 000 mm 3 ).The created part is from the top of the skull near to the supraorbital ridge.The thickness of skull is variable.The minimum thickness is in the temporal bone (∼ 2 mm) and the maximum value is on the back of the occipital bone (∼ 15 mm).
Using automatic segmentation in STL Model Creator, the STL model of the external shape of the skull (cortical bone) was created.The internal structure and the cancellous bone tissue were created using manual segmentation.Surfaces were created from the STL file in the CATIA software and then volume model was created from these surfaces.It was necessary to make space (gaps) in compact bone tissue volume for cancellous bone tissue, which was represented by another volume model.In cortical bone tissue representation, gaps were created using Boolean operations.These gaps were filled by another volume model of the cancellous bone tissue (as shown in Fig. 3).
Four basic geometry models were created as follows (see also Figs. 2 and 3 for model designation): • a shell model, created only from outer surface (consisting of only cortical bone tissue); • volume model (consisting of only cortical bone tissue); • volume model including cortical and cancellous bone tissue; • model, which combines shell elements of cortical bone and volume elements of cancellous bone.

Meshing procedure
The geometry model was meshed in Ansys 17.2.Volume models (A1, A2, B and E) were discretized using quadratic tetrahedral elements SOLID187.The total number of elements was approximately around 2 000 000, which corresponds approximately to 3 000 000 nodes (model A1, A2 and E) and 1 500 000 elements, which corresponds to 2 300 000 nodes (model B).
The shell model including only the cortical bone tissue (model C) was discretized using SHELL181 elements with a constant thickness of 10 mm.This corresponds to an average thickness of the created part of the skull bone.The total number of elements was approximately around 18 000, which corresponds approximately to 17 000 nodes.
Other variant was a combination of solid elements (cancellous bone tissue) surrounded by shell elements (cortical bone tissue).In this study, 2 different thicknesses of cortical bone tissue were assumed and evaluated.In literature, thickness of cortical bone tissue is stated to be approx.2 mm [22].It corresponds to model D2.Based on CT dataset, thickness of cortical bone tissue for another computational model was assumed 1 mm (as demonstrated in Fig. 1), which is represented by model D1.The total number of elements and nodes was 3 600 000 and 4 000 000, respectively.The operations necessary to create the models are summarized in Table 1.The overview of solved variants is shown in Figs. 2 and 3.
* -material properties calculation; ** -defining material properties, FE mesh, defining loads and boundary conditions etc.; *** -creation of shell elements around solid model and meshing them

Materials
Two different levels of material model which are the most frequently encountered in literature in FE simulations of human skull were used as follows.

Homogeneous material
The most commonly used model of material is homogeneous, isotropic, linearly-elastic material, which is defined by Young's modulus (E [MPa]) and Poisson ratio (μ [-]) [27].Values of these parameters used in present study are shown in Table 2.

Heterogeneous material
Experiments proved a correlation between bone tissue density and CT number (pixel intensityin this case, there were 4 096 shades of grey).Similarly, a correlation between bone tissue density and Young's modulus can be found.In previous research, a lot of formulas for determination Young's modulus from bone tissue density [3].In the current study, we used the formula (1) taken from [4,9,11] where E is Young's modulus in GPa, C is Carter-Hayes correlation constant of value 3.79 and ρ a is apparent density in g/cm 3 , which is calculated from formula (2).The formula is valid for bone tissue density from 0.07 to 2 g/cm 3 .Material properties based on CT numbers (as demonstrated in Fig. 3) were mapped into the FE model (model E) using in-house software application CTPixelMapper [2].Thanks to this approach, skull sutures with specific material properties different from those of the surrounding bone can be taken into account.This phenomenon is difficult or impossible to consider in the other types of computational models.

Boundary conditions and loads
Loads which are assumed to act on skull in this study are intracranial pressure (ICP) and pressure loading force (as shown in Fig. 4).The skull was loaded from the inside by ICP of 2 kPa (15 mmHg), which is a physiological value of ICP [7].The static force was applied to an area (∼ 68 mm 2 ) as shown in Fig. 4. The value was assumed 1 000 N. This value is based on preliminary calculations assuming a wide range of loading force (50-6 000 N).The selected value is believed to provide results that best point out on the variability among all studied models.According to literature, this force is lower than the one, which can cause fractures of skull [17,20].All models were fixed in the same way at the bottom side of the modelled part of the skull.Given the purpose of the study and because the evaluation paths introduced further are sufficiently far from the constraints, the simplification of boundary condition is justifiable.

Evaluation of FEA results
In literature, the most commonly used parameters in biomechanical analysis are displacement and von Mises strain [8,20,24].It has been reported that bone tissue is pathologically overloaded at strain values higher than 0.004 [-].In an extreme case (strain value above 0.025 [-]), cracks and fractures may occur [8].
Preliminary calculations were performed to decide whether large deflection effects were significant in the studied cases.Since the differences between small deflection and large deflection options were detected to be up to 20% (total displacement at model D1), the large deflection effects were considered to be significant and, therefore, were taken into account in all calculations.

Results
For a comparison, total displacements and von Mises strains were assessed.

Total displacement
The distribution of total displacement is shown in Fig. 5.In all analysed variants, the highest displacements were located in the region of loading force application.The highest maximum values of total displacement were detected in model D1 (1.04 mm) and E (0.78 mm).In case of model D2, the maximum value of total displacement was 0.37 mm.Maximum values of displacement for the volume models with homogeneous materials were 0.26 mm (model B), 0.13 mm (model A1) and 0.48 mm (model A2).The maximum value of total displacement was the lowest in case of model C (0.03 mm).

Von Mises strain
The distribution of von Mises strain in all analysed variants is shown in Fig. 6.The highest maximum values of von Mises strain were detected in model E (0.048).In models D1 and D2, the strains were 0.021 and 0.010.In case of solid models, the strains were 0.003 (model A1), 0.009 (modelA2) and 0.011 (model B).The lowest maximum value of strain was detected in model C (0.001).For a comprehensive comparison, von Mises strains were evaluated along paths of interest.The paths are curves which were defined as intersections of parallel planes at certain heights (three different heights as shown in Fig. 7) on the outer and inner side of the skull.
In total, there were six different paths (three from inside and three from outside the model).The middle path had considerably higher values of strain (as demonstrated in Fig. 8).External paths had higher values of strain than internal ones at the same height.The internal middle The values of strain in all analysed models along the middle path are shown in Figs. 9 and 10.Along most of the path length, the strains are not significantly different among the studied variants.However, near the location of the loading force (approx.20% of the path length), the strain values as well as the differences among the variants are significantly higher.The values of strain are significantly higher there.
Model C had the lowest values of strain (the maximum was 0.001.Volume model A1 had the highest maximum value 0.002, volume model B 0.007, (which is triple of that of A1), and model A2 0.008.Model D1 and D2 performed as the least stiff ones.Model D1 had higher a maximum value of strain (0.078) than model D2 (0.026).Model E had the highest value of strain 0.018 of all analysed variants.

Discussion
All computational models have the same input, CT dataset.The process of creating geometry of computational model is the most time consuming and significantly differs for various model types.Different computational models also require differently long time to compute the results.From the model-creation point of view, computational models A1, A2, and C took the same amount of time.Both models D (D1 and D2) took longer to create than A-models, due to solid and shell nature of FE mesh.Computational model E took longer to finalize too due to bigger complexity of material properties.The longest process of model creation was associated with computational model B, because it consisted of two volumes, which were segmented separately.Cancellous bone tissue had to be segmented manually, which took much longer than automatic segmentation of A-models.All other required operations were done twice.The difficulty of computational model creation and computation would become even more apparent if analyses of the skull including cranial implant with fixation (mini-plates and micro-screws) would be performed.
The highest values of total displacement were found near the acting force location.The highest maximum value was observed in model D1 and the lowest was in model C. The difference between the maximum values of displacement in those two cases was 97%.This difference and the maximum values for each model show that the level of computational model directly affects observed parameters and it is necessary to select input parameters appropriately, considering the actual purpose of the simulation (e.g.model C for some dynamic simulations, model E or A2 for stress-strain analysis of cranial implant fixed by mini-plates and micro-screws).
Von Mises strain describes state of bone tissue that might be assessed using criteria for bone overloading or damage.The highest value of strain was detected in models E and D1, which were extremely overloaded.On the other hand, the lowest values were in computational models C (∼ 0.000 8 values of model E are the highest on both paths.Almost similar values of strain to model D2 had volume model A2 along both paths.Similar values were also in models A and C.These models are both very stiff because the poor material properties of the cancellous bone tissue (as listed in Table 1) were not considered.Values of strain in model B are higher because the cancellous bone tissue was considered in that case.Despite this fact, values of strain were lower than in shell models D1 or D2.Distribution of von Mises strain is similar in all analysed models except for model E where higher values appeared not only around the loading force, but in skull sutures as well.It can be caused by different material properties (Young's modulus) of the bone tissue around the sutures.This was the only model of all studied cases, which was able to capture this phenomenon.
In previous studies, various computational models were used depending on the purpose of the analysis.The shell model of cortical bone tissue with a thickness of 2 mm (D2) was used very often (e.g. in [1,21,22]).Based on CT data, it was observed in this study that the thickness of cortical bone tissue is lower (D1) than usually used in other studies.Solid models (A1, A2 and B) were used very frequently in [6,12,18,19,25].The computational model, which has the material properties based on CT data (model E), is a modern approach, which is still not so often used in literature for analysing skull defects.
Research which focuses on mechanical interactions with skull [1,2,5,6,12,18,20,25,26] uses all types of models analysed in this study, except for model C. Model C is the stiffest one and values of analysed parameters may be quiet far from reality when used for an inappropriate objective.Solid model A1 is closer to reality because the geometry is similar to a real skull (the thickness is not constant), but it is too stiff.Model A2 considers spongy part in the value of Young's modulus, which is ∼ 3.5 times lower than in model A1.Solid model B takes into account cancellous bone geometry and material properties.If computational model is stiff, values of analysed parameters are lower, which can lead to underestimation of analysed situation.On the other hand, models D1 and D2 consist of shell-like cortical bone tissue with solid-based cancellous bone tissue inside the entire skull bone, which is far from reality as well.The percentage ratio of cortical and cancellous bone tissue is higher in these cases than in reality, which is resulting in higher pliability than in models A1, A2, B and C. It may lead to overestimation as well.
Model E which assumes heterogeneous material properties based on CT data represents a much more complex approach than the other models.Despite the fact that heterogeneous distribution of material characteristics is more realistic in this case compared to all previous models, Carter-Hayes formula was used to calculate Young's modulus should be used with caution as it was not directly for the skull, but for a generic bone (including both, cortical and cancellous bone tissue).This drawback cloud be solved by deriving a formula specifically for human skull.
Model E was the only model, which was able to detect another strain concentrators, such as skull sutures.Results of total displacement and von Mises strain were similar to those of shell model D2, but if more complex computational models (including implant, fixation mini-plates and screws) were investigated, using shell elements is would not be entirely applicable, because the interaction between microstructures (such as micro-screws) and cortical bone (based on shell elements) would be difficult to capture.

Conclusion
In this study, seven computational models on different modelling levels were investigated.When it comes to appropriate model type/level selection, it must be considered what is the actual objective of the study.It is evident that different computational models showed different complexity of creation process, different computational demands and results precision.It is, therefore, important to carry out thorough pre-analysis considerations and to choose appropriate computational model type.Nowadays, interaction of cranial bone with implant, fixation miniplates and micro-screws is a typical example of situation that is very sensitive to the proper selection of the modelling level.
In case of using volume models with homogeneous distribution of material properties, cancellous bone role should be considered for the solution with regards to trade-off between the model construction difficulty and the resulting parameters demand.In case of detailed investigation of strains in skull, model E model was found to be beneficial, because it takes into account heterogeneous structure of bone tissue, which can lead to possibility of detection of another risk locations and stress concentrators.

Fig. 1 .
Fig. 1.Cancellous bone tissue distribution within the skull, upper view (left).CT-image of cranial bone detail including the estimation of cortical bone thickness (right)

Fig. 2 .Fig. 3 .
Fig. 2. Model variants: (A1 andA2) Volume model of skull with cortical bone material properties (without cancellous bone); (B) Volume model of skull including cortical and cancellous bone (separate volumes); (C) Shell model of skull with constant thickness 10 mm; (D1) Combination of volume model (cancellous bone) and shell model (cortical bone) -1 mm thickness; (D2) is the same as D1 but with a thickness of 2 mm

Fig. 4 .
Fig. 4. Model of geometry and loading forces -Intracranial pressure (ICP) and static force F

Fig. 9 .
Fig. 9. Values of von Mises strain of each analysed model in location most influenced by loading force on external middle path.The curves are distinguished by several marks.Note that the marks do not highlight any particular results but they serve only for better clarity [-]) and A1 (∼ 0.002 4 [-]).The highest values of strain on paths were in model D1 on the external path and in model E on the internal path.The lowest values were in model C.Both shell models (D1 and D2) had similar values of strain along the path to the heterogeneous model E. Values of model E are between the values of both shell models in the location of loading force along the external path.For the internal path, values of model E are the highest.Except for the location of loading force,

Table 1 .
Operations necessary to create computational model (CM)

Table 2 .
Material properties