Development of a Finite Element Model of Decompressive Craniectomy

Decompressive craniectomy (DC), an operation whereby part of the skull is removed, is used in the management of patients with brain swelling. While the aim of DC is to reduce intracranial pressure, there is the risk that brain deformation and mechanical strain associated with the operation could damage the brain tissue. The nature and extent of the resulting strain regime is poorly understood at present. Finite element (FE) models of DC can provide insight into this applied strain and hence assist in deciding on the best surgical procedures. However there is uncertainty about how well these models match experimental data, which are difficult to obtain clinically. Hence there is a need to validate any modelling approach outside the clinical setting. This paper develops an axisymmetric FE model of an idealised DC to assess the key features of such an FE model which are needed for an accurate simulation of DC. The FE models are compared with an experimental model using gelatin hydrogel, which has similar poro-viscoelastic material property characteristics to brain tissue. Strain on a central plane of the FE model and the front face of the experimental model, deformation and load relaxation curves are compared between experiment and FE. Results show good agreement between the FE and experimental models, providing confidence in applying the proposed FE modelling approach to DC. Such a model should use material properties appropriate for brain tissue and include a more realistic whole head geometry.


Introduction
Uncontrolled brain swelling and raised intra-cranial pressure can lead to death or poor functional outcome in patients suffering from severe traumatic brain injury, ischaemic stroke and other type of brain insults [1,2]. Decompressive craniectomy (DC) is a surgical procedure which involves the following steps: opening of the skin (scalp incision), removal of a large piece of skull (bone flap), opening of the outermost membrane covering the brain (dura mater) and closure of the skin incision. The creation of an opening in the skull essentially turns the cranial cavity from noncompliant into compliant [3]. As a result, intra-cranial pressure is reduced after a DC [4,5]. A number of clinical studies have been launched recently in an attempt to collate high-quality evidence regarding the clinical effectiveness of this operation in patients suffering from traumatic brain injury and ischaemic stroke [6][7][8]. While evidence is still accumulating, it is important that further work is undertaken in an attempt to better characterise the effects of DC on the brain. This is especially important for clinicians, as the optimal parameters of craniectomy (particularly size and location) remain controversial [9].
In the context of brain injury and swelling, DC will induce mechanical strain in the brain tissue which may cause damage [10]. Engineering models of DC can provide predictions of the induced strain and hence identify regions of the brain in which damage may occur. Although there have been various finite element (FE) models of the brain during surgery (e.g. [11]), there is only one model in the literature that focuses on decompressive craniectomy [12]. Gao and Ang [12] present a three dimensional model of the brain and skull with material properties taken from Cheng and Bilston [13]. However their model does not consider how details of the material model or geometric features affect the accuracy of the predicted strain fields in the brain, nor do they provide experimental validation. Indeed clinical data to support such models are difficult to obtain, reinforcing the need for a careful validation study of the FE approach before applying it to simulate DC. Strains extracted from patient CT scans (see [14] for methods) could be used to compare full skull geometry models. However post-op CT scans generally take place a number of days after surgery; therefore early tissue response to the craniectomy would be lost.
One aspect that needs to be considered carefully is the type of material model required. In decompressive craniectomy there is a hybrid of confined and unconfined loading, which means that conclusions drawn on appropriate material models from other brain deformation studies may not be valid for DC. Kyriacou et al. [15] suggest that a poro-viscoelastic (PVE) model of the brain combines the advantages of both poroelastic models at low strain rates and viscoelastic models for high strain rates. PVE parameters for the brain are available from Cheng and Bilston [13] and Franceschini [16].
This paper aims to identify an appropriate FE model for DC. Experimental tests on an idealised DC geometry, using gelatin hydrogel to represent brain tissue, are compared with corresponding FE models with a PVE material model. Rather than adopting the complex geometry of the surgical situation, the approach used in this paper is to consider a simplified axisymmetric model of DC. Validation of the model is based on load response and strain comparisons between FE and experimental models. Good correlation between these experiments and FE models can provide confidence in adopting the proposed FE model, using material models representative of brain tissue and realistic head geometries, to analyse the effects of DC on the brain.

Methods
Typically a DC has an approximately circular or elliptical (for unilaterial DC) geometry [9]. This suggest an idealised model of DC using a cylindrical geometry for the 'skull' and a circular craniectomy opening. The experimental model and FE approach are described below and results for the two approaches are then compared.

Experimental model
The experimental model of the idealised craniectomy with semicylindrical geometry of radius 30 mm and height 20 mm is illustrated in Figure 1. Gelatin was used to simulate the brain tissue, while a semi-circular aluminium mould modelled the constraint of the skull. A flat perspex front plate was in contact with the gelatin to constrain the movement while keeping the front face visible. It is known that gelatin hydrogels exhibit PVE behaviour [17] and gelatin has also been used as a surrogate for brain in experimental brain models [18,19]. While not an exact match for brain tissue, gelatin provides insight into how a PVE material would respond under similar loading situations to DC. In clinical practice intra-cranial pressure causes deformation of the brain out through the craniectomy opening. In the experimental simulation this deformation is instead associated with a reduction in the size of the cavity, induced by movement of a top platen containing the semicircular craniectomy opening. There was a fillet radius r of 0:5 mm around the edge of the craniectomy opening and craniectomy radii of 5, 10 and 15 mm were used.
Manufacture. The experimental model used 8% 180 g bloom porcine gelatin hydrogels (Sigma Aldrich), made in a similar manner to the method described by Galli and Oyen [20]. Samples were created by casting molten gelatin into silicone moulds, acetate was then placed on the top surface preventing a meniscus forming. The gels were refrigerated at 4uC overnight for ease of extraction from the moulds.
After extraction, the gels were marked with dots of Indian Ink (Winsor and Newton) on the front face, see Figure 1. The application of these markers on the flat face allowed image analysis of the strain on a plane in the centre of the simulated craniectomy for comparison with the FE model.
The gels were left overnight in a sealed container for the ink to dry whilst minimising sample dehydration. The samples were tested 43 hours +1 hour after the gels were originally cast.

Materials characterisation
Gelatin hydrogels can be produced with varying concentrations of gelatin and water producing large variations in material parameters. Therefore it is necessary to characterise the specific gelatin hydrogels used. For this type of gelatin, confined compression tests approximately isolate the poro-elastic behaviour Gels were manufactured in the same manner as for experimental tests. Cylindrical confined and unconfined compression samples were cast with radii of 5 mm and heights of 5 mm and 10 mm respectively. Unconfined compression tests were carried out in displacement control. Strain was applied over a 10 s ramp to a peak strain of 10%, followed by a hold time of 600 s. Confined compression tests were load controlled. After initial contact, load was increased over a 2:5 s ramp to a peak load of 2 N which was held for 10,000 s. The PVE material model used in this paper is a well-established approach to considering time-dependent material behaviour. The following paragraphs in this section provide an overview of the model and detail the procedures used to identify the material parameters. References are given for a more detailed description of the model.
The governing equations for a PVE material are as follows: the poro-elastic component of the response is characterised by the instantaneous Young's modulus E of the solid, the Poisson's ratio n and a permeability k k. The time domain visco-elasticity of the solid is introduced by a relaxation function which is applied to the solid component of the material, as follows [21]: The g i and k i terms represent dimensionless shear and bulk relaxation moduli respectively and t g i are corresponding time constants, with K t ð Þ and G t ð Þ being time-dependent bulk and shear moduli. The PVE model used in this paper assumes g i~ki and t k i~t g i for all time constants. This is the same model as the ''BPVE2'' model used by Suh and DiSilvestro [22]. It is also assumed that the relative relaxation rates of visco-elasticity and poro-elasticity differ significantly in unconfined compression, so that the visco-elastic response can be determined from unconfined compression tests.
The unconfined compression tests were analysed using a standard three parameter Prony series form of the relaxation function to solve the Boltzmann hereditary integral [23], as per [24] by Equation 2.
where C i and t i are the Prony parameters, h max the maximum displacement, z the sample height, D the sample diameter, t the time and RCF is a ramp correction factor defined as: Matlab (The Mathworks, Natick, MA, US) was used to identify the parameters C 0 , C i and t i which best fit the measured responses using Equation 2. The dimensionless shear relaxation moduli g i can be determined from the C i terms directly as g i~C Confined compression tests were analysed by curve-fitting the consolidation curves assuming a Terzaghi solution in a similar manner to Franceschini [16]. The Terzaghi solution takes the form: where f tÃ ð Þ is: following Detournay and Cheng [25]. Here the consolidation parameter u Ã is a non-dimensional settling of the sample under representing a non-dimensional time, n is the Poisson's ratio with n u the undrained Poisson's ratio, L is the layer thickness and G is the shear modulus. The constant c (assuming incompressible constituents) is given by: Again Matlab was used to identify the material parameters E, n and k k which best fit the measured response in the above equation.   Table 1 details the values of the poro-elastic and visco-elastic material parameters derived from these confined and unconfined compression tests.
Test protocol. This section details the protocol used for the deformation of the idealised DC model illustrated in Figure 1. The aluminium platen was fixed to an Instron 5544 universal testing machine (Instron, UK) fitted with a 500 N load cell. A compressive pre-load of 0:5 N was applied to the platen in order to locate the surface of the sample and fully confine the sample in the rig. The need for a preload was highlighted by Cheng and Bilston [13] (for unconfined compression of a PVE material). A total displacement of 0:25 mm was then applied at a speed of 0:025 mm s 21 , which is a similar speed to that used by Miller and Chinzei [26] and representative of surgical strain rates. This was followed by a hold period of 600 s at the final displacement of 0:25 mm. This hold period provides substantial visco-elastic and poro-elastic relaxation of the gelatin without dehydrating the sample. All tests were carried out for a minimum of 3 repeats per opening size.
Imaging. Images of the front face of the gelatin were captured at 0:5 s intervals during testing using a camera (Pixelink, Ottawa, Canada) and analysed using Matlab.
The markers on the front face of the gelatin were tracked throughout the experiment using methods developed by Crocker [27]. The centres of the markers in each frame determine the path of the markers in successive images given constraints on the maximum distance each particle can travel.
Markers which were tracked for less than 20% of the frames were discarded as noise. The remaining paths were fitted with a second-order polynomial in both the x and y directions. This both smooths the output and interpolates the location of the markers in any frame where they are missing.
The resultant fitted paths of the markers were used to determine the applied strain on the front face. First, Delaunay triangulation [28] was used to define a triangular mesh which described the marker positions in the first time-step. This mesh was then used to find the strain in each element from the displacement of each marker using the methods of Screen and Evans [29]. Figure 2 illustrates the axisymmetric FE model of an idealised craniectomy which was developed. The model has a similar cylindrical geometry and circular craniectomy opening as for the experimental model. The craniectomy opening a and the fillet radius r were defined as parameters to simplify a parametric study. The craniectomy opening a was varied between 8 and 15 mm in 1 mm increments with r~0:5 mm and the fillet radius r was varied from 0:3 mm to 1 mm in 0:1 mm increments with a~10 mm. The FE models with a smaller than 8 mm failed to converge due to excessive strains under the craniectomy edge. Abaqus Standard [21] with non-linear geometry was used in all simulations.

Finite element model
Materials. The measured material properties for gelatin given in Table 1 were used in a baseline set of FE calculations for comparison with the experiments. In addition, the influence of the poro-elastic and the visco-elastic components was analysed by varying material parameters from these base-line values. A linear PVE model was used as preliminary simulations showed that using a non-linear PVE model had little effect on the model response to loading, presumably as regions of high strains were comparatively small.
Mesh. Axi-symmetric quadratic elements with reduced integration (Abaqus element CAX8RPH) were used in all FE models. A fine mesh was created directly under the craniectomy edge with an average element length of 0:005 mm. This was graded to the edges where elements had a side length of 0:6 mm. The mesh used  Loading and constraints. Displacement boundary conditions of zero x-displacement on the vertical edges and zero ydisplacement on the bottom horizontal surface were applied, see Figure 2. The platen underwent a 0:25 mm y displacement at a speed of 0:025 mm s 21 and was constrained to displace only in the vertical direction with no rotation. The platen was held at this position for 600 s. It was assumed that the interface between the PVE material and the platen was a frictionless, hard contact.
Boundaries were considered to be permeable in all simulations; this assumption is discussed in the results section. Load response curves were normalised by peak loads in order to analyse the effect of a, independent of the elastic response, as follows:

Finite element results
This normalisation, by eliminating the effect of craniectomy opening a on peak load, highlights the shape of the response, see Figure 3B. Figure 4 compares the load-time response for visco-elastic, poroelastic and poro-viscoelastic FE models having the same material properties as the gelatin PVE model. There are two distinct phases of relaxation. The first phase before a time of 50 s has a shape characteristic of visco-elasticity, followed by a second poro-elastic consolidation phase after 50 s. A similar form of load-time response was seen in Figure 3B for the effect of the craniectomy radius a. The curves of Figures 3B and 4 initially fall on a master curve but diverge with increasing time. The divergence occurs when the poro-elastic relaxation dominates the visco-elastic relaxation for tw50 s. For these later times, an increase in craniectomy size a slows the relaxation of the sample.   Figure 5A shows that, as the fillet radius r on the craniectomy increases, the force decreases slightly. The shapes of the curves in Figure 5A are similar and the normalisation of the force collapses the responses to a single curve, see Figure 5B. Hence changing r has little effect on the relaxation response of the material. However, the peak load and peak shear strain do depend on r. The peak strains E 12 under the craniectomy decrease from {0:82 to {0:49 for r increasing from 0:3 to 0:9 mm. Figure 6A shows the load versus time response for individual tests at three values of a tested in the experimental model. There is a decrease in load as a increases, with the same two-phase relaxation with the visco-elastic and poro-elastic portions. Figure 6B shows the result of the normalisation of Equation 7 for the experimental results. The curves collapse in a similar manner to the FE results, with divergence from this single curve at later times.

Comparison between FE and experimental models
Load comparison. Figure 7 compares FE load response curves with the average experimental curves for a~10 mm. The dashed lines in Figure 7 represent the average experimental results + the standard deviation. The FE simulations follow a similar path to the experimental data, within the range of experimental results.
Material parameters were varied in the FE model to match the variation in the parameters identified by compression tests, see Table 2. This variation in parameters produces a spread of FE curves which approximately reproduces that in the experiments ( Figure 8A), suggesting that the experimental variation could be largely due to inherent variation in material parameters for gelatin. Figure 8B shows the effect of varying the dimensionless shear modulus g 1 on the predicted load response. An increase in g 1 increases early relaxation (including during the loading phase) as seen in Figure 8B. The relaxation shape of the experiments can be matched by varying g 1 . However, an increase in Young's modulus E would be necessary in order to fully match the response if a larger g 1 term were chosen.
Strain and deformation comparison. Figures 9A and B compare the experimental and FE shear strains on the front face (experiment) and central plane (FE) at the end of the loading step for a craniectomy radius a~10 mm. Figures 9C and D plot the strain variation along the lines x~10 mm and y~10 mm for both the experimental and FE cases.
The width, depth and overall shear strain values are comparable between the experimental and FE models. The FE model has smaller features of high shear strain in the region directly below the craniectomy edge (yw15 mm in Figure 9C). Results (not shown here) are similar for the larger 15 mm opening. Figure 10A illustrates the vertical y-displacements sampled on a vertical centreline on the front face of the gelatin. Figures 10B and  C show the y-displacements along horizontal lines with y~10 mm and 20 mm, i.e. the middle and the top surface of the front face. Along all sampled lines there is close agreement between the FE model with permeable boundaries and experiment.
Boundary conditions. The FE model with impermeable boundary conditions predicts deformations that are much larger in the bulge region, see Figure 10, and elevated loads (12 N as against 8 N for permeable condition and experiment). The results with a permeable condition are used in all other comparisons as this provided much closer agreement between the FE and experiments.

Discussion
This paper develops an idealised FE model of DC, validating this with experimental measurements on a model with similar geometry. Whilst this study is a step towards understanding DC and the resultant strain, specific clinical inferences should not be taken from the model as the materials and dimensions of the model differ from the brain geometry. Rather the aim of the modelling work is to determine an appropriate modelling approach to be applied in the clinical application.
Complexities of using a PVE model are explored with constraints similar to DC. An area of particular interest is the mechanical strain applied to the brain under different conditions. This is especially important for clinicians, as the optimal parameters of craniectomy (particularly size and location) remain controversial.
Finite element predictions of the load-time response (Figure 7) match the experiment within experimental variability. The characteristic two-phase relaxation seen in experiments is also seen in the FE model. When only the peak load is required, Figure 4 shows that a visco-elastic model is adequate. From a clinical perspective the conditions at peak load are likely to  Table 2. Extreme PVE parameters from unconfined and confined compression testing of gelatin, as used in Figure 8A.  Figure 4). In this case it is likely that physiological factors would significantly modify the time-dependent material response, so that modifications to the PVE model would be required to capture these factors. Significantly smaller displace-ments have been applied in this model in comparison with those likely in clinical conditions. However the relative importance of viscoelastic and poro-viscoelastic deformations is expected to be similar at larger strains, so that the conclusions of this study can be applied to clinical models of the procedure. Shear strains on the front face of the gelatin experimental model derived using image analysis were compared with FE predictions. There is a reasonable agreement between strains in the FE and experiment, although the spacing of the markers in the experi-  mental tests precluded the measurement of the highest strain region as this region is smaller than the spacing between markers. It should be noted here that a linear PVE material has been used throughout the FE modelling, but peak shear strains in some simulations are higher than a linear limit. Comparison with a neo-Hookean model of the gelatin at the initial and infinite times (where n is known) revealed that, despite this limitation, the peak shear strains for the hyper-elastic case were within 5% of the equivalent linear elastic model ({0:793 vs. {0:761). This highlights the need for preliminary FE models as presented in this paper to examine the modelling assumptions, showing in this case that a linear elastic model is adequate to capture the details of the strain response in DC. While the material properties of this study are not those appropriate to brain tissue, the strains observed in our tests are comparable with strains expected to cause microstructural damage in brain tissue (of the order of 0.05-0.35 [30]), confirming that the study explores a regime of deformation which is clinically relevant.
The modelling work assumed that the material was uniform and isotropic, matching the experimental behaviour of the gelatin. In practice brain mechanical properties will depend on position and are likely to exhibit an anisotropic dependence of damage on strain, associated with the anatomical and axonal structure of the brain. Wright and Ramesh [30] have shown in a study of traumatic brain injury how such factors can be incorporated into a model of brain damage associated with deformation in DC. A challenge in such work would be identifying the appropriate damage parameters associated with the timescales of deformation, although the micromechanical modelling approach of Cloots et al [31] provides an attractive way of helping determine these parameters.
It was shown that the permeable boundary condition matched the experimental data better than the impermeable condition. Flow on the boundary could be caused by a local dehydration of the gelatin. Although samples were kept in a sealed container, there may have been surface dehydration prior to testing. This surface dehydration could lead to a region of the gel which would allow the fluid from the interior to flow to the exterior in the same manner as a permeable boundary condition. Inspection of the boundary nodes in the FE models with permeable boundary conditions showed that, for a typical craniectomy size, the modelling predicts a fluid flow of 1.5 ml over the first 100 s of relaxation. This equates to a fluid build up of 0.15 mm over the surface of the model during that 100 s interval. In the experimental situation, the sample is kept hydrated and lubricated by applying a thin layer of water on the surface of the model prior to loading. It is therefore suggested that the small volume of fluid due to poro-elastic fluid flow could have been accommodated in the test rig despite the impermeable nature of the test rig walls. The permeable boundary may also replicate the surface of the brain. Fluid collections often develop on the surface of the brain near the craniectomy between the brain and the scalp [32] which may suggest a similar permeable surface behaviour. In any case results confirm that careful attention needs to be paid to this aspect of the model in the clinical application.
Parametric studies have highlighted the significant effects of craniectomy size and fillet radius on the deformation and force response for the assumed model material properties. Hence a study using brain-specific material properties should include these effects to provide guidance on optimising the clinical procedure for DC.

Conclusions
The reasonable correlation between the FE model of the idealised craniectomy and the experimental model using gelatin confirms that the proposed FE model has the potential to provide useful predictions for the brain tissue strains developed during decompressive craniectomy procedures. The work has provided clear guidance on the requirements of such a model as applied to the clinical situation, in addition to the prerequisites of appropriate material and geometric properties. In order to obtain the peak response and critical strain during DC a visco-elastic approach is adequate; only to model the subsequent transient decay of strains is a poro-viscoelastic model required. In this case it is likely that physiological factors would significantly modify the time-dependent material response, so that modifications to the PVE model would be required to capture these factors. A linear elastic model is adequate to capture the peak strains during the procedure, but close attention should be paid to the boundary conditions at the brain/skull interface. In addition a clinical study should examine the effect of craniectomy opening and fillet radius on the developed strains.