Coupling Finite Element and Mesh-free Methods for Modelling Brain Deformation in Response to Tumour Growth

,

In biomechanics of soft tissues, it is common to encounter extreme deformations that cannot be handled by traditional modelling methods, such as the Finite Element method.An example of this is brain tumour growth.Very little is known within this field and such a model should be useful for medical use, particularly in quantifying tumour aggressiveness, therapy planning, as well as surgical planning and simulation.The three dimensional mechanical response of the brain is highly non-linear, involving extremely complex constitutive models and geometry, which is very time consuming to model using public Finite Element (FE) software (Miller, Taylor, et al., 2005).Furthermore the Finite Element method on its own will be inaccurate and problematic for modelling the brain deformation response to tumour growth, since the mesh surrounding the tumour is easily distorted, consequently destroying elements.Complicated re-meshing can combat this, however it is extremely time consuming.Alternatively a solely mesh-free Element Free Galerkin (EFG) model will be capable of handling larger deformations and topology changes (Li and Liu, 2004).Despite this the method is known to be less reliable than Finite Element analysis for moderate deformations and suffers from Dirichlet boundary difficulties (Fries and Matthies, 2003).A coupled Finite Element / Element Free Galerkin approach is proposed to overcome the shortcomings of each individual method, by placing a mesh-free domain around the tumour affected location, with the remaining brain tissue modelled as a hexahedral mesh.

Background Theory
The Finite Element method is a numerical approach for solving systems of partial differential equations, by discretising the domain into small volumes (elements) and estimating the solution in each of the elements via shape functions.The estimated solutions are then substituted into integral differential equations of the weak form with the residuals minimised (Bathe, 1996).The Element Free Galerkin method (Horton, 2006) conducts the same process, without requiring the connectivity of elements.Shape functions are not within elements but small neighbourhoods of nodes, called support domains, each of which is associated to an integration point (Belytschko, Krongauz, et al., 1996).The Moving Least Squares formulation is used to minimise residuals within the EFG method.

Coupled Finite Element / Element Free Galerkin Method
Mixed-mesh coupling is achieved by constructing interface support domains in between FE and EFG boundaries, as shown in Figure 2.1.Interface support domains are created by allowing the EFG nodal support domains to extend into the FE region, consuming nodes.They follow the same numerical approach as the EFG method.FE nodes that exist within a nodal support domain are considered by both the EFG and FE methods with their nodal forces summed together.

Total Lagrangian Formulation
The Total Lagrangian (TL) formulation is a general convergent method for dealing with materially nonlinear effects and large deformations.It is derived from the Principle of Virtual Displacements.The TL formulation references all static and kinematic variables to the initial configuration of the system, which means all derivatives with respect to spacial coordinates can be precomputed.Hence, the number of mathematical operations performed in each time step of the TL algorithm is reduced, and is therefore favourable for surgical simulation (Bathe, 1996, andMiller, Joldes, et al., 2007).

Explicit Integration
The method of integration for the TL formulation is very important for the performance of the algorithm.Computing the deformation field and internal forces within a soft tissue organ, such as the brain, requires the utilisation of an efficient integration method for integrating static or dynamic equations within the time domain.Both implicit and explicit integration methods are suitable for this problem.The Central Difference Method of explicit integration is chosen over implicit integration as it does not require the solution of simultaneous non-linear equations and performs much faster on materials of very low stiffness, such as the brain (Miller, Joldes, et al., 2007).
The Central Difference Method is derived from Newton's second law as shown, where, M represents the nodal mass, Δt the simulation time-step, F i n  t and R i  n t the net force and reaction force on the node.The previous, current, and future nodal displacements are described by Unfortunately two issues arise with the use of explicit integration.Firstly explicit methods are conditionally stable, requiring time-steps below a critical value (Joldes, 2006).Furthermore explicit methods cannot obtain static solutions unless damping is used.Quasi-static simulations can be conducted such that the load is applied very slowly over a long period of time.The longer the simulation time, the closer the solution is to a static one.

Total Lagrangian Explicit Dynamics Algorithm
The Total Lagrangian Explicit Dynamics (TLED) algorithm combines the TL formulation with explicit integration.The TLED algorithm has been successfully utilised for non-linear soft tissue deformations in FE and EFG numerical approaches individually (Horton, 2006, andMiller, Joldes, et al., 2007).Both numerical methods have provided reliable and accurate results using the TLED algorithm, thus it is desirable for use in the coupled FE/EFG approach.The following psuedocode describes the TLED algorithm in relation to solving mixed-mesh problems:

Pre-process:
1 . Combine local reaction forces to obtain net nodal reaction forces at time, t. 10.Explicitly determine displacements, u t+Δt , using Central Difference formula: ( t Ri -Total nodal reaction force at time, t)

Program Implementation
The implementation of the tailor made localised soft tissue deformation simulator is divided into three main sections: • Preprocessor -Reads in mixed mesh model and constraints, pre-computing all initial configuration stationary properties.• Analysis Solver -Executes the main time loop performing calculations in accordance with the Total Lagrangian Explicit Dynamics Algorithm.• Postprocessor -Uses series of visualisation tools to view and identify implications of analysis.
The script language, MATLAB, was chosen for the preprocessing and postprocessing stages of the simulator.MATLAB is a very powerful, high level, language, containing many built in functions, which are of particular use for dealing with matrices.This was beneficial for the preprocessing stage, which requires many large matrix operations.Furthermore MATLAB holds significant advantages for the postprocessing phase as it has excellent visualisation tools, allowing for advanced analysis of the results.Unfortunately MATLAB performs much slower than compiled languages, hence it was not feasible for use in the implementation of the analysis solver.The functional programming language, C, was chosen as it is very fast, with inbuilt optimisation compilation abilities.

Preprocessor
The implementation of the preprocessor phase can be broken down into a number of smaller subsections.It should be noted that the computational performance of this phase is less important than the analysis solver, since for coarse mixed-meshes, the run time is negligible in comparison.

Simulation Properties
All of the major simulation properties are user defined and must be set prior to running the preprocessor.This includes the maximum displacement, time-step, and the total simulation time.The simulation properties provide enough information to setup the deformation loading curve, for applying incremental displacements at each time-step.The default deformation loading curve is defined in (4.1), where T is the total simulation time.

Mixed-Mesh Reader
The nodes, elements, and boundary conditions of the mixed-mesh are read in from ABAQUS output files.Two ABAQUS output files are required, both containing nodal positions and boundary conditions for each method, with the FE output file containing additional information about element composition.The data format of the information read from the FE and EFG ABAQUS files is displayed in Table 4.1.A coupled list of all nodes, Xcoupled, is then formed by combining XFE and XEFG in that order.

Constrained FE Nodes
Displaced FE Nodes XFE (Nnodes-FE x 3) matrix of all FE node locations.Each row of the matrix corresponds to the FE node number.
EFE (Nelements-FE x 8) matrix of all FE elements.Each row of the matrix corresponds to the element number, containing eight node numbers forming a hexahedral.
FE_node_fix (Nnodes-FE x 3) binary matrix of all FE nodes, 1 represents if the node is fixed for that dimension.
FE_node_disp (Nnodes-FE x 3) binary matrix of all FE nodes, 1 represents if the node is displaced for that dimension.

EFG Nodes
Constrained EFG Nodes Displaced EFG Nodes XEFG (Nnodes-EFG x 3) matrix of all EFG node locations.Each row of the matrix corresponds to the EFG node number.
EFG_node_fix (Nnodes-EFG x 3) binary matrix of all EFG nodes, 1 represents if the node is fixed for that dimension.
EFG_node_disp (Nnodes-EFG x 3) binary matrix of all EFG nodes, 1 represents if the node is displaced for that dimension.Table 4.1: Data format of information read in from the ABAQUS output files.
An integration point grid for the EFG and interface region should also be read in.This is just a matrix holding the three dimensional coordinates of the integration point locations.It is ideal to have a regular grid, such that each integration point can be assigned the same volume.

Material Model
The material model information is to be set, requiring Young's modulus, Poisson's ratio, and the density for each material identified within the mixed-mesh.From this data the lame` material constants can be calculated as follows.
The material model is based upon the Neo-Hooken model (Bathe, 1996).

Support Domains
The construction of nodal support domains for the EFG and interface region is quite simple following on from Horton (2007).The method requires a fixed number, n, of nodes per support domain, which is user defined.Support domains are then constructed by finding the n closest nodes to each integration point.A limit on the number of Finite Element nodes allowed within a single support domain removes the possibility of an entire element being consumed by a support domain, which would have no hourglass control measures.Having a fixed number of nodes per support domain is faster and more robust than typical support domain constructions, which rely on defining a fixed local volume with a varying number of nodes.

Hexahedral Shape Functions
The hexahedral shape functions and derivatives are determined from a series of calculations.The matrix of hexahedral shape function natural derivatives is defined as, using the node numbering convention as described in Bathe (1996).Each elemental Jacobian, J, is then calculated based on the element nodal position vector,  x .J =∂hr T  x (4.5) Using the elemental Jacobians, hexahedral shape function derivatives, ∂ h , are then computed by, ∂ h=∂ hr J −1  T (4.6)

Moving Least Squares Shape Functions
Calculations of the Moving Least Squares shape functions for each nodal support domain of size n, are derived from Fries and Matthies (2003).Consider a three dimensional space vector of monomial basis functions, p, of length m.
The nodal displacement approximation, u h (x), is calculated with respect to the coefficient vector, a(x). i.e.
The formulation of a(x) is determined by minimising the weighted residual function, J, where, In (4.8) W(di) represents a weight function with di being the distance between the node, xi, and integration point, x.Minimising J is done by considering which leads to the following linear relationship, A x a  x =B x U (4.10)In (4.10)A is an m x m matrix known as the moment matrix defined by, B is an m x n matrix given by, and U is the vector of length n as shown, (4.15) where the shape function vector Φ of length n at the i th node in the support domain is given by, The length m of p is user defined and should be chosen such that shape functions are all interpolated in a similar fashion in each dimension.There is a trade-off between the total number of integration points and the size of m due to the limitations on computational speed.Single point integration is well suited for low order interpolations, hence a lower value of m is chosen, while using a larger number of support domains.More support domains relieve the emphasis on stress calculations at any integration point.It has been found in Horton (2007) that setting m = 4 and using 8 nodes per support domain (n = 8) is substantial enough for deformation to be transferred between support domains.In addition it has been noted that by using very small support domain sizes, the weighting of each node can be considered equal, without having a negative impact on the accuracy of the solution.This reduces the risk of generating singular matrix A.

Mass Allocation
Initially a matrix, MFE, is setup for handling the mass of all nodes within the FE domain.The mass of each node within an element for MFE is calculated using the determinant of each elemental Jacobian from (4.5) and the material density, ρ.
(4.17)The nodal contributions for all elements are then summed up to give MFE.
A coupled mass matrix, Mcoupled, is then created for allocating masses to all EFG and interface nodes involved in support domains.Each integration point is allocated a volume and consequently a mass based upon the materials density.This mass is equally divided amongst the number of nodes within the support domain, where, n represents the number of nodes per support domain, and V (g) is the volume of the specific integration point, g.The mass of each FE node in MFE is then added to Mcoupled, giving the entire nodal mass of the system.This is a very effective method of distributing mass throughout the EFG/interface region since nodes that appear in more support domains will receive more forces.One concern, however, is that nodes not included in many support domains will have a low mass, which can result in unbalanced forces and high accelerations.This is undesired, often leading to unstable simulations.It can be avoided by involving each node in at least two or three support domains as suggested in Horton (2007).A further measure is implemented so that any node that manages to escape support domain allocation is removed to prevent massless nodes entering the analysis.

Analysis Solver
The analysis solver is the most computationally intensive phase of the simulator.It consists of the main time loop described in the Total Lagrangian Explicit Dynamics Algorithm with a few additional considerations.Efficient programming is very important to minimise the number of calculations required in the main time loop, substantially increasing the performance of the algorithm.

Main Time Loop
Both the Finite Element and Element Free Galerkin methods follow the same calculations for the main time loop of the TLED algorithm, making it quite easy to implement the coupling as treating the entire domain as a single method.-Inverse Right Cauchy-Green deformation tensor

Hourglass Control
One of the biggest disadvantages to using single-point integration for hexahedral elements is the requirement for controlling zero energy modes, known as hour-glassing (Hallquist, 2006).In order to control hour-glassing within the Finite Element domain, resistance providing artificial stiffness is implemented, which has a negligible effect on stable global modes.This is an efficient method following on from modifications of Flanagan and Belytschko (1984) perturbation method (Joldes, Wittek, et al., 2007).An additional hourglass control force is added to the total force of the system, based on the hourglass resistance and displacement.Hourglass base matrix is setup as: The k th elemental hourglass control force, F hg , can then be calculated by the following series of equations: where ∂ h  k  hg is the k th elemental hourglass shape function derivative, u ij hg represents the elemental hourglass displacement derivatives, and Rhg is the hourglass resistance constant.The hourglass force is then added to the elemental force calculated without hourglass control as shown in (4.23).
A good value for the hourglass resistance, Rhg, was found to be Rhg = 0.04/9.

Mixed-Mesh Generation
One of the main considerations of utilising coupled FE/EFG methods is the generation of mixed-meshes.Meshless methods may have much more freedom with node placement in comparison to hexahedral FE methods, although very regular node placements can lead to singular shape functions within one plane (Horton, 2007).Regular node and integration arrangements for meshless methods have been successful, however this becomes equivalent to a hexahedral mesh, suffering from hour-glassing and furthermore is quite complicated to generate over irregular domains, such as the brain (Li, and Liu, 2004).A way around this is to model the meshless region using a well refined tetrahedral mesh, which helps to maintain a roughly even density, while not conforming to very regular node placements.The tetrahedral element information is omitted, since meshless methods do not require the connectivity of nodes.It is possible to convert sections of a hexahedral mesh into tetrahedrals in commercial meshing software such as Hypermesh, which easily gives rise to the generation of mixed-meshes.Unfortunately this procedure relies on the existence of a Finite Element mesh, which, depending on the complexity of the shape, may still be quite complicated and time consuming to create.
Another aspect of mixed-mesh generation is the requirement for an integration point domain.As found in Horton (2006), a background grid of integration points for the EFG domain leads to greater accuracy and stability, while still performing efficiently.The background grid can be extended up to the boundary between the EFG and FE region in order to maximise coupling as mentioned in Chapter 5.2.

Coupling Integration Point Distribution
An investigation into the depth of coupling and its effects on the accuracy of the method was conducted.The level of coupling is indicated by the number of FE nodes involved in interface support domains.This is strongly dependent on the distribution of integration points around the boundaries of each region.The best results occur when the coupling integration point grid is mapped up to the boundary of the FE region.This maximises the number of FE nodes within interface support domain.
To avoid unnecessarily using too many integration points around the coupling region the findings of Horton ( 2007) should be employed with respect to interface support domains.That is, by using twice as many interface support domains to interface nodes, the Moving Least Squares approximation will still give very accurate results.

Interface Support Domain Composition
An area of interest in the coupling method is the number of allowable FE nodes to be considered within a single interface support domain.A support domain should not be allowed to contain only Finite Element nodes, as there is the possibility that it has coincided with an entire element, which would not have any hourglass control measures.For this reason a limit on the maximum number of FE nodes within an interface support domain must be defined, based on the number of nodes within a hexahedral element.This condition would then be checked when allocating support domains within the preprocessor.
It was found that by reducing the amount of allowable FE nodes per interface support domain causes the level of coupling to decrease, which produces results with poorer accuracy.Allowing up to 7 FE nodes per interface support domain, one node less than a hexahedral element, maximises the coupling, consequently leading to increased accuracy of solutions.

Validation
The coupling method was validated by a series of quasi-static deformation tests, using the material properties of healthy brain tissue.Initially the algorithm was trialled and compared against a Finite Element solution using commercial software (ABAQUS) for a homogeneous cylinder undergoing compression, extension, and shear deformations.The mixed-mesh contained an outer FE region with an inner EFG core.The results were further compared against validated Element Free Galerkin software, showing that the coupling method performs slightly better than a solely EFG method and is still very close to the FE solution.Comparisons against a FE solution is preferable over an analytical solution due to the complicated mathematical nature of analytical methods and general loss of accuracy for non-linear problems.The results shown in Table 6.1 reflect the maximum nodal displacement error in comparison to the FE ABAQUS solution.
The final validation test involved a partially constrained ellipsoid undergoing indentation on the surface.The mixed mesh gave highly accurate results in comparison with a FE mesh simulated in ABAQUS.

Deformation Model
Maximum Error (mm) It is evident from Table 6.1, Figure 6.2, and Figure 6.3 that the maximum error in all cases falls within the allowable 0.85 mm tolerance for surgical accuracy (Bourgeois, Magnin, et al., 1999).An additional investigation has been conducted comparing tumour growth on an ellipsoid using a stand alone FE mesh and a mixed-mesh, containing an EFG region of high density surrounding the proposed area of localised high deformation.The accuracy of the FE results became questionable as tumours grew larger than 523.6 mm 3 .Figure 6.4 demonstrates the localised deformation to the FE mesh with a series of increasing tumour growths.The onset of hour-glassing is present during the very early stages of tumour growth, despite control measures in place to prevent this.Table 6.2 compares the maximum nodal displacements, surrounding the localised high deformation region, using the stand alone FE mesh and the mixed-mesh.It is apparent that initially both methods give quite similar results, however they begin to differ significantly as the distortion to the FE mesh increases.Discrepancies between the two methods are observed for tumour growths greater than 523.6 mm 3 .The FE mesh fails as the tumour reaches a volume of 3053.6 mm 3 .Significant distortion to the FE mesh is present in Figure 6.5, with the hexahedral elements compressing up to 70%, well beyond the reliable limits as discussed in Wittek, Dutta-Roy, et al. (2008).The mixed-mesh deformation for the same tumour growth volume is shown in Figure 6.6.

Tumour Growth Analysis
A mixed-mesh of a brain was created allowing for tumour growth to occur behind the left ventricle, mimicking an MRI scan of a tumour affected brain in Urbach, Binder, et al. (2007).An EFG nodal domain of high density surrounds the proposed tumour region allowing for large deformation.The tumour was grown as an ideal sphere, of which the analytical equations are well defined.
For healthy brain tissue and tumour we assume that Young's modulus, E, is 3000 Pa, and Poisson's ratio, v, is 0.49 (Miller, Chinzei, et al., 2000, Miller, 2002).The ventricles contains cerebro-spinal fluid (CSF), which has very similar material properties to water, hence they are modelled as a soft elastic compressible solid, with E = 10 Pa and a low Poisson's ratio, v = 0.1.A low Poisson's ratio allows to simulate leakage of the cerebro-spinal fluid which may occur under static deformation conditions.(Wittek, Miller, et al., 2006).
A number of different tumour growth sizes were investigated, with deformation volume change in the ventricles, from an initial volume of 57.1 ml, noted in Given the location of the tumour, the ventricular deformation and associated volume loss, displayed in Table 6.3, is likely to correspond to the leakage of CSF.In reality CSF may leak between the left and right ventricles, however it is also known to leak out of the ventricles completely, particularly under static deformation (Rando and Fishman, 1992).Furthermore, large tissue deformation is apparent, particularly for the 4118 mm 3 tumour, as shown in Figure 6.7, with local displacements of up to 9.66 mm.For this example 14.5% of the brain mesh experienced displacements greater than 5 mm.
This follows on from Clatz, Bondiau, et al. (2004), which declared volume variation within the ventricles and large tissue deformation in response to brain tumour growth mass effect.It should be noted that simulated tumour growths larger than 14000 mm 3 on the given mixed-mesh reduced the reliability of the method as the resulting deformation to the Finite Element region became too large.A greater EFG domain would be required surrounding the tumour affected region in order to simulate larger growths accurately.

Conclusion
A new coupling method has been proposed to combine the Finite Element and Element Free Galerkin methods for modelling the non-linear soft tissue deformation of the brain in response to tumour growth.The method was verified against FE commercial software and a validated EFG simulator on a number of different mixed meshes.All results were very accurate, easily falling within the 0.85 mm error tolerance, corresponding to the working resolution of an MRI scan.Simple analytical tumour growths were conducted on a comprehensive brain mesh.The tumour's close proximity to the ventricles caused observable volume changes, which may involve leakage of CSF.Furthermore large tissue displacements were noted, with a significant portion of the brain undergoing moderate deformation.In reality this may have a detrimental effect on the cell metabolism and function of the brain, altering the stress distribution and blood flow.Further investigation into realistic tumour growth models and implementation of a brainskull contact algorithm would increase the reliability of the results.Ultimately this would become beneficial for both clinical prognosis and operation planning as well as for simulated training applications.

Figure 2
Figure 2.1: Coupled FE/EFG domain with interface region and support domains highlighted.

Figure 6
Figure 6.3: Cross-section comparison of deformed boundary for indented ellipsoid.
λ -Lame` parameters; J -Jacobian determinant) 4.13)By finding the inverse of A equation (4.14) can be solved,

X 0 t =u i , j I 3x3
Psuedocode below presents the implementation of the main time loop.
t, were defined.

Table 6 .
1: Maximum displacement errors in coupling method compared against FE ABAQUS solutions.