Abstract
This paper describes a new method for simulating tissue growth which can handle anisotropic changes in volume. The method takes advantage of the total Lagrangian formulation which allows the computation of nodal forces for each element in a finite element mesh based on a theoretical stress-free configuration, obtained by considering the unconstrained anisotropic growth of the considered element. The method allows the modelling of shrinking (atrophy), swelling, or tissue growth and the computation of the resulting mechanical stresses in the surrounding tissue. The steady-state solution is obtained using an explicit integration method and dynamic relaxation. The proposed method allows the coupling of continuum mechanical simulations with underlying growth mechanisms, offering a tool for the multiscale study of tissue growth.
1. Introduction
The majority of studies using finite element models to simulate growth in soft tissues have used implicit time integration schemes to solve the equations of motion. We propose an explicit method for modelling anisotropic growth of soft tissues that uses dynamic relaxation to obtain a steady-state solution. Explicit time integration has been used to model soft tissue deformations for neurosurgery [1, 2], tissue swelling [3], and cardiac remodelling [4]. Explicit methods have the advantage of decoupling the equations of motion, therefore being very well suited for parallel implementation [5, 6].
Many growth processes are anisotropic whereby the growth rates are not equal in all directions. Examples include cardiac growth and remodelling, skin growth and expansion, tendon and muscle fibre growth, and tumour growth. Most tissue can be modelled using hyperelastic constitutive laws [7].
Approaches that model growth by modifying the Jacobian, as done in the commercial finite element software Abaqus for temperature induced volumetric changes [8], for example, cannot handle anisotropy. Therefore we apply the multiplicative decomposition of the deformation gradient introduced by Rodriguez [9] to separate the deformation into growth and elastic components. The use of the total Lagrangian formulation allows us to define the nodal forces based on the known initial configuration and greatly simplify the solution procedure.
2. Materials and Methods
2.1. Anisotropic Tissue Growth Modelling
Prior to growth the tissue is found in a physiological state considered to be the stress-free configuration. Anisotropic growth introduces changes in the stress state of the tissue, even if the motion is not constrained, as different tissue regions can grow with different rates and in different directions.
We consider the deformation from an initial configuration to a final configuration as in Figure 1. The stress-free configuration after growth can be obtained for each element if we consider the growth is unrestricted and uniform over the element. The stress-free configuration is a hypothetical configuration for each element, and in these configurations elements are disconnected, as they may grow at different rates and therefore compatibility between elements will be lost. Growth can be expressed as a deformation gradient describing the changes in geometry of the element with coordinates and in the configurations and , respectively. Based on the definition of the deformation gradient and applying the chain rule of deformation we arrive at the multiplicative split introduced by Rodriguez [9]: To satisfy the requirement that zero stress corresponds to zero strain the constitutive law must be written in terms of the elastic deformation gradient The total deformation gradient can be expressed as where is the displacement from the initial to the final configurations, is the vector of unknown nodal displacements from the initial to the final configurations, and is the vector of shape functions associated with the nodal displacements in the initial configuration.
Hyperelastic material laws are defined in terms of a strain energy function . The second Piola-Kirchhoff stress associated with the stress-free configuration can then be computed as where is the elastic Green-Lagrange strain and is the elastic right Cauchy-Green deformation tensor.
The nodal forces for an element, created by the internal stresses, are computed as where are the coordinates in the configuration and are the element shape functions.
We consider the transform , for an individual element. By performing a change of variables in (6) we can express the nodal forces based on the known initial configuration: If we select the shape functions in the configuration as the shape function derivatives can be computed as By combining (3), (7), and (9) the nodal forces can be calculated as We notice that (10) is written only based on the deformation gradient from the initial to the final configuration, which is a function of the unknown displacements, the growth deformation gradient, and the material law; also the integral is defined over the known initial configuration. By solving the equilibrium between the nodal forces given by (10) and the externally applied forces, the unknown displacements can be computed.
2.2. Solution Method
The timescale of soft tissue growth is on the order of hours or days; therefore we are interested in static solution. We use explicit integration in the time domain with dynamic relaxation to solve the model. This approach was developed for soft tissue simulation during neurosurgery on commodity hardware [1, 10] and has recently been applied to the modelling of soft tissue swelling and shrinking [3].
The solution algorithm is based on the total Lagrangian formulation and uses explicit integration to solve an artificial dynamic equation [1, 11–13]. Dynamic relaxation is applied to obtain the steady-state solution [10, 13]. Dynamic relaxation involves the addition of mass proportional damping to the equation of motion, which damps down the free oscillations leading the solution towards its steady state. We use a lumped (diagonal) mass matrix which enables the decoupling of the finite element equations and the use of explicit integration, leading to a simple solution procedure which does not involve the solution of large systems of equations. Mass and mass proportional damping have no influence on the final steady state; therefore they are chosen to guarantee stability and fast convergence toward the steady-state solution.
Explicit integration has many advantages over the alternative implicit integration methods; it does not require storage of large stiffness matrices, simple handling of material, and geometric nonlinearities and can handle very large problems; stepping in time requires only vectorial operations making it easy to implement and well suited for parallel implementation. Nevertheless, explicit integration is only conditionally stable, and therefore the time step used for integration cannot be larger than a maximum allowable time step in order for the solution to be stable (the Courant stability condition).
The proposed method is not the only alternative for finding the solution. The equilibrium equation between the nodal forces and the externally applied forces can be solved by assembling the system’s stiffness matrix (the derivative of the nodal forces with regard to nodal displacements) and using Newton-Raphson iterations. A detailed example of how to compute the consistent stiffness matrix for an equation similar to (10) is given in [14].
2.3. Numerical Experiments
To demonstrate its capacity to model anisotropic growth we performed a simple 2-dimensional numerical experiment showing the growth of a rectangular tissue sample rotated by 45 degrees which doubles its width. Therefore, the deformation gradient that describes the growth can be defined as with a rotation matrix (which aligns the local and global coordinate systems) and a growth matrix. By setting we expect no growth in the local direction and doubling of the width in the local direction.
We also performed a constrained shrinking experiment, to demonstrate the capability of the method to handle both growth and shrinkage. By setting we expect no growth in the local direction and halving of the width in the local direction.
The model geometry consists of a rectangle with a width of 4 cm and a height of 10 cm. The material behaviour is characterised using a hyperelastic neo-Hookean constitutive law. The material represents brain tissue with Young’s modulus of 3000 kPa and Poisson ratio 0.49 [15].
Explicit solution methods are only conditionally stable, with the maximum allowable time step defined for each element in the mesh by the material properties and the element geometry (characteristic length). Because the growth deformation gradient is expressed as a combination of rotations and stretches and given the fact that rotations do not modify the geometry of the element, only the stretches must be taken into account. Therefore we used the following procedure to determine the stable time step: we considered the characteristic length after growth equal to the characteristic length before growth multiplied with the minimum of stretches.
3. Results and Discussion
The results of the first numerical experiment show that the sample grows as expected (Figure 2), doubling in size in the local direction, when the movement is not constrained.
By constraining the nodes on one end of the sample the growth is restricted and causes stress in the material (Figure 3). The stresses are higher at the constrained end of the sample and almost zero at the free end. Both growth and shrinkage of the sample are properly handled by the method.
The integral defining the nodal forces, given by (10), involves a nonlinear, nonpolynomial integrand. Exact integration can only be obtained for uniform strain elements (such as the linear triangle/tetrahedron elements). Any higher order elements will require complex integration schemes to ensure integration accuracy, such as adaptive integration [16].
The stability conditions for explicit time integration, used in our solution method, are derived based on the theory of linear elasticity. For highly nonlinear material behaviour, involving material stiffening, the computed allowable time step may be overestimated, leading to the loss of solution stability. This has been generally handled by slightly reducing the computed allowable time step by using a safety factor. Other possible solutions include accurate (but time consuming) computing of the stable time step for each element based on the eigenvalue of the element’s stiffness matrix or the use of a different solution method altogether, such as implicit Newton-Raphson iterations.
4. Conclusions
We developed a method for modelling anisotropic growth of soft tissues. The method is based on the deformation gradient decomposition and the use of total Lagrangian formulation. The novelty of the method consists in the rewriting of the nodal forces based on the known, fixed initial configuration, instead of the unknown stress-free configuration. This is achieved by using a mapping of the shape functions between the two configurations and a change of variable in the total Lagrangian definition of nodal forces. The resulting nodal forces definition can be easily implemented into an existing total Lagrangian explicit dynamics code or can be used with implicit Newton-Raphson procedures. The implementation does not require any changes to existing material laws. Because the initial configuration is fixed, the recomputation of shape functions or redefinition of integration points and weights are avoided during the computation of the solution.
Numerical experiments were performed which demonstrated the capacity of the method to compute anisotropic growth/shrinkage of a tissue sample. By using definitions of growth which lead to known deformed configuration we were able to verify that the algorithm behaves as expected.
Due to the explicit nature of the solution method used, the algorithm has the potential to be implemented on graphics processing units (GPUs) to considerably increase the computational speed. The algorithm can be applied to model soft tissue growth and remodelling. The introduction of stress or strain dependent growth within the model could be pursued in the future.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The financial support of the Australian Research Council (discovery Grant no. DP120100402) and the National Health and Medical Research Council (Grant no. APP1063986) is gratefully acknowledged.