Simulation of planar soft tissues using a structural constitutive model: Finite element implementation and validation

Computational implementation of physical and physiologically realistic constitutive models is critical for numerical simulation of soft biological tissues in a variety of biomedical applications. It is well established that the highly nonlinear and anisotropic mechanical behaviors of soft tissues are an emergent behavior of the underlying tissue microstructure. In the present study, we have implemented a structural constitutive model into a finite element framework specialized for membrane tissues. We noted that starting with a single element subjected to uniaxial tension, the non-fibrous tissue matrix must be present to prevent unrealistic tissue deformations. Flexural simulations were used to set the non-fibrous matrix modulus because fibers have little effects on tissue deformation under three-point bending. Multiple deformation modes were simulated, including strip biaxial, planar biaxial with two attachment methods, and membrane inflation. Detailed comparisons with experimental data were undertaken to insure faithful simulations of both the macro-level stress–strain insights into adaptations of the fiber architecture under stress, such as fiber reorientation and fiber recruitment. Results indicated a high degree of fidelity and demonstrated interesting microstructural adaptions to stress and the important role of the underlying tissue matrix. Moreover, we apparently resolve a discrepancy in our 1997 study (Billiar and Sacks, 1997. J. Biomech. 30 (7), 753–756) where we observed that under strip biaxial stretch the simulated fiber splay responses were not in good agreement with the experimental results, suggesting non-affine deformations may have occurred. However, by correctly accounting for the isotropic phase of the measured fiber splay, good agreement was obtained. While not the final word, these simulations suggest that affine fiber kinematics for planar collagenous tissues is a reasonable assumption at the macro level. Simulation tools such as these are imperative in the design and simulation of native and engineered tissues.     

Finite element implementation of a new model of slight compressibility for transversely isotropic materials

Modelling transversely isotropic materials in finite strain problems is a complex task in biomechanics, and is usually addressed by using finite element (FE) simulations. The standard method developed to account for the quasi-incompressible nature of soft tissues is to decompose the strain energy function (SEF) into volumetric and deviatoric parts. However, this decomposition is only valid for fully incompressible materials, and its use for slightly compressible materials yields an unphysical response during the simulation of hydrostatic tension/compression of a transversely isotropic material. This paper presents the FE implementation as subroutines of a new volumetric model solving this deficiency in two FE codes: Abaqus and FEBio. This model also has the specificity of restoring the compatibility with small strain theory. The stress and elasticity tensors are first derived for a general SEF. This is followed by a successful convergence check using a particular SEF and a suite of single-element tests showing that this new model does not only correct the hydrostatic deficiency but may also affect stresses during shear tests (Poynting effect) and lateral stretches during uniaxial tests (Poisson's effect). These FE subroutines have numerous applications including the modelling of tendons, ligaments, heart tissue, etc. The biomechanics community should be aware of specificities of the standard model, and the new model should be used when accurate FE results are desired in the case of compressible materials.     

Intervertebral disc creep behavior assessment through an open source finite element solver

Degenerative Disc Disease (DDD) is one of the largest health problems faced worldwide, based on lost working time and associated costs. By means of this motivation, this work aims to evaluate a biomimetic Finite Element (FE) model of the Intervertebral Disc (IVD). Recent studies have emphasized the importance of an accurate biomechanical modeling of the IVD, as it is a highly complex multiphasic medium. Poroelastic models of the disc are mostly implemented in commercial finite element packages with limited access to the algorithms. Therefore, a novel poroelastic formulation implemented on a home-developed open source FE solver is briefly addressed throughout this paper. The combination of this formulation with biphasic osmotic swelling behavior is also taken into account. Numerical simulations were devoted to the analysis of the non-degenerated human lumbar IVD time-dependent behavior. The results of the tests performed for creep assessment were inside the scope of the experimental data, with a remarkable improvement of the numerical accuracy when compared with previously published results obtained with ABAQUS^®. In brief, this in-development open-source FE solver was validated with literature experimental data and aims to be a valuable tool to study the IVD biomechanics and DDD mechanisms.     

Simulated bioprosthetic heart valve deformation under quasi-static loading

For more than 40 years, the replacement of diseased natural heart valves with prosthetic devices has dramatically extended the quality and length of the lives of millions of patients worldwide. However, bioprosthetic heart valves (BHV) continue to fail due to structural failure resulting from poor tissue durability and faulty design. Clearly, an in-depth understanding of the biomechanical behavior of BHV at both the tissue and functional prosthesis levels is essential to improving BHV design and to reduce rates of failure. In this study, we simulated quasi-static BHV leaflet deformation under 40, 80, and 120 mm Hg quasi-static transvalvular pressures. A Fung-elastic material model was used that incorporated material parameters and axes derived from actual leaflet biaxial tests and measured leaflet collagen fiber structure. Rigorous experimental validation of predicted leaflet strain field was used to validate the model results. An overall maximum discrepancy of 2.36% strain between the finite element (FE) results and experiment measurements was obtained, indicating good agreement between computed and measured major principal strains. Parametric studies utilizing the material parameter set from one leaflet for all three leaflets resulted in substantial variations in leaflet stress and strain distributions. This result suggests that utilization of actual leaflet material properties is essential for accurate BHV FE simulations. The present study also underscores the need for rigorous experimentation and accurate constitutive models in simulating BHV function and design.     

Effects of boundary conditions on the estimation of the planar biaxial mechanical properties of soft tissues

Evaluation and simulation of the multiaxial mechanical behavior of native and engineered soft tissues is becoming more prevalent. In spite of this growing use, testing methods have not been standardized and methodologies vary widely. The strong influence of boundary conditions were recently underscored by Waldman et al. [2002, J. Materials Science: Materials in Medicine 13, pp. 933-938] wherein substantially different experimental results were obtained using different sample gripping methods on the same specimens. As it is not possible to experimentally evaluate the effects of different biaxial test boundary conditions on specimen internal stress distributions, we conducted numerical simulations to explore these effects. A nonlinear Fung-elastic constitutive model (Sun et al., 2003, JBME 125, pp. 372-380, which fully incorporated the effects of in-plane shear, was used to simulate soft tissue mechanical behavior. Effects of boundary conditions, including varying the number of suture attachments, different gripping methods, specimen shapes, and material axes orientations were examined. Results demonstrated strong boundary effects with the clamped methods, while suture attachment methods demonstrated minimal boundary effects. Suture-based methods appeared to be best suited for biaxial mechanical tests of biological materials. Moreover, the simulations demonstrated that Saint-Venant's effects depended significantly on the material axes orientation. While not exhaustive, these comprehensive simulations provide experimentalists with additional insight into the stress-strain fields associated with different biaxial testing boundary conditions, and may be used as a rational basis for the design of biaxial testing experiments.     

Computational efficiency of numerical approximations of tangent moduli for finite element implementation of a fiber-reinforced hyperelastic material model

In this study, we evaluated computational efficiency of finite element (FE) simulations when a numerical approximation method was used to obtain the tangent moduli. A fiber-reinforced hyperelastic material model for nearly incompressible soft tissues was implemented for 3D solid elements using both the approximation method and the closed-form analytical method, and validated by comparing the components of the tangent modulus tensor (also referred to as the material Jacobian) between the two methods. The computational efficiency of the approximation method was evaluated with different perturbation parameters and approximation schemes, and quantified by the number of iteration steps and CPU time required to complete these simulations. From the simulation results, it can be seen that the overall accuracy of the approximation method is improved by adopting the central difference approximation scheme compared to the forward Euler approximation scheme. For small-scale simulations with about 10,000 DOFs, the approximation schemes could reduce the CPU time substantially compared to the closed-form solution, due to the fact that fewer calculation steps are needed at each integration point. However, for a large-scale simulation with about 300,000 DOFs, the advantages of the approximation schemes diminish because the factorization of the stiffness matrix will dominate the solution time. Overall, as it is material model independent, the approximation method simplifies the FE implementation of a complex constitutive model with comparable accuracy and computational efficiency to the closed-form solution, which makes it attractive in FE simulations with complex material models.     

Prediction of brain deformations and risk of traumatic brain injury due to closed-head impact: quantitative analysis of the effects of boundary conditions and brain tissue constitutive model

In this study, we investigate the effects of modelling choices for the brain–skull interface (layers of tissues between the brain and skull that determine boundary conditions for the brain) and the constitutive model of brain parenchyma on the brain responses under violent impact as predicted using computational biomechanics model. We used the head/brain model from Total HUman Model for Safety (THUMS)—extensively validated finite element model of the human body that has been applied in numerous injury biomechanics studies. The computations were conducted using a well-established nonlinear explicit dynamics finite element code LS-DYNA. We employed four approaches for modelling the brain–skull interface and four constitutive models for the brain tissue in the numerical simulations of the experiments on post-mortem human subjects exposed to violent impacts reported in the literature. The brain–skull interface models included direct representation of the brain meninges and cerebrospinal fluid, outer brain surface rigidly attached to the skull, frictionless sliding contact between the brain and skull, and a layer of spring-type cohesive elements between the brain and skull. We considered Ogden hyperviscoelastic, Mooney–Rivlin hyperviscoelastic, neo–Hookean hyperviscoelastic and linear viscoelastic constitutive models of the brain tissue. Our study indicates that the predicted deformations within the brain and related brain injury criteria are strongly affected by both the approach of modelling the brain–skull interface and the constitutive model of the brain parenchyma tissues. The results suggest that accurate prediction of deformations within the brain and risk of brain injury due to violent impact using computational biomechanics models may require representation of the meninges and subarachnoidal space with cerebrospinal fluid in the model and application of hyperviscoelastic (preferably Ogden-type) constitutive model for the brain tissue.     

An anatomically based patient-specific finite element model of patella articulation: towards a diagnostic tool

A 3D anatomically based patient-specific finite element (FE) model of patello-femoral (PF) articulation is presented to analyse the main features of patella biomechanics, namely, patella tracking (kinematics), quadriceps extensor forces, surface contact and internal patella stresses. The generic geometries are a subset from the model database of the International Union of Physiological Sciences (IUPS) (http://www.physiome.org.nz) Physiome Project with soft tissue derived from the widely used visible human dataset, and the bones digitised from an anatomically accurate physical model with muscle attachment information. The models are customised to patient magnetic resonance images using a variant of free-form deformation, called 'host-mesh' fitting. The continuum was solved using the governing equation of finite elasticity, with the multibody problem coupled through contact mechanics. Additional constraints such as tissue incompressibility are also imposed. Passive material properties are taken from the literature and implemented for deformable tissue with a non-linear micro-structurally based constitutive law. Bone and cartilage are implemented using a 'St-Venant Kirchoff' model suitable for rigid body rotations. The surface fibre directions have been estimated from anatomy images of cadaver muscle dissections and active muscle contraction was based on a steady-state calcium-tension relation. The 3D continuum model of muscle, tendon and bone is compared with experimental results from the literature, and surgical simulations performed to illustrate its clinical assessment capabilities (a Maquet procedure for reducing patella stresses and a vastus lateralis release for a bipartite patella). Finally, the model limitations, issues and future improvements are discussed.     

Dynamic finite element implementation of nonlinear,anisotropic hyperelastic biological membranes

We present a novel method for the implementation of hyperelastic finite strain, non-linear strain-energy functions for biological membranes in an explicit finite element environment. The technique is implemented in LS-DYNA but may also be implemented in any suitable non-linear explicit code. The constitutive equations are implemented on the foundation of a co-rotational uniformly reduced Hughes-Liu shell. This shell is based on an updated-Lagrangian formulation suitable for relating Cauchy stress to the rate-of-deformation, i.e. hypo-elasticity. To accommodate finite deformation hyper-elastic formulations, a co-rotational deformation gradient is assembled over time, resulting in a formulation suitable for pseudo-hyperelastic constitutive equations that are standard assumptions in biomechanics. Our method was validated by comparison with (1) an analytic solution to a spherically-symmetric dynamic membrane inflation problem, incorporating a Mooney-Rivlin hyperelastic equation and (2) with previously published finite element solutions to a non-linear transversely isotropic inflation problem. Finally, we implemented a transversely isotropic strain-energy function for mitral valve tissue. The method is simple and accurate and is believed to be generally useful for anyone who wishes to model biologic membranes with an experimentally driven strain-energy function.     

Dynamic finite element modeling of poroviscoelastic soft tissue

Clinical evidences relative to biomechanical factors have demonstrated their important contribution to the behaviour of soft tissues. Finite element (FE) analysis is used to study the mechanical behaviour of soft tissue because it can provide numerical solutions to problems that are intractable to analytic solutions. This study focuses on the development of a FE model of a poroelastic biological tissue, which incorporates the viscoelastic material behaviour, finite deformation and inertial effect. The FE formulation is based on the weak form derived from the governing equation, and Newmark-beta method as well as Newton's method is incorporated into the implicit non-linear solutions. One-dimensional analytical solutions were used to verify the theoretical formulation and the numerical implementation of the proposed model. This study was further extended to analyze two-dimensional biomechanical models and the results clearly demonstrate the importance of including finite deformation, viscoelasticity and inertial effects.     

A simple method to estimate the exponential material parameters of heart valve tissue based on analogy between uniaxial tension and micropipette aspiration

Micropipette aspiration (MA) has been widely used to measure the biomechanical properties of cells and biomaterials. To estimate material parameters from MA experimental data, analytical half-space models and inverse finite element (FE) analyses are typically used. The half-space model is easy to implement but cannot account for nonlinear material properties and complex geometrical boundary conditions that are inherent to MA. Inverse FE approaches can account for geometrical and material nonlinearities, but their implementation is resource-intensive and not widely available. Here, by making analogy between an analytical uniaxial tension model and a FE model of MA, we proposed an easily implementable and accurate method to estimate the material parameters of tissues tested by MA. We first adopted a strain invariant-based isotropic exponential constitutive model and implemented it in both the analytical uniaxial tension model and the FE model. The two models were fit to experimental data generated by MA of porcine aortic valve tissue (45 spots on four leaflets) to estimate material parameters. We found no significant differences between the effective moduli estimated by the two models ( $p > 0.39$ ), with the effective moduli estimated by the uniaxial tension model correlating significantly with those estimated by the FE model ( $p < 0.001; R^{2}= 0.96$ ) with a linear regression slope that was not different than unity ( $p = 0.38$ ). Thus, the analytical uniaxial tension model, which avoids solving resource-intensive numerical problems, is as accurate as the FE model in estimating the effective modulus of valve tissue tested by MA.     

The influence of fiber dispersion on the mechanical response of aortic tissues in health and disease: a computational study

Changes in the structural components of aortic tissues have been shown to play a significant role in the pathogenesis of aortic degeneration. Therefore, reliable stress analyses require a suitable and meaningful constitutive model that captures micro-structural changes. As recent data show, in-plane and out-of-plane collagen fiber dispersions vary significantly between healthy and aneurysmatic aortic walls. The aim of this study is to computationally investigate the influence of fiber dispersion on the mechanical response of aortic tissues in health and disease. In particular, the influence of three different fiber dispersions is studied: (i) non-rotationally symmetric dispersion, the most realistic assumption for aortic tissues; (ii) transversely isotropic dispersion, a special case; (iii) perfectly aligned fibers (no dispersion in either plane), another special case. Explicit expressions for the stress and elasticity tensors as needed for the implementation in a finite element code are provided. Three representative numerical examples are studied: planar biaxial extension, inflation of residually stressed and pre-stretched aortic segments and inflation of an idealized abdominal aortic aneurysm (AAA) geometry. For the AAA geometry the case of isotropic dispersion is additionally analyzed. Documented structural and mechanical parameters are taken from human aortas (healthy media/adventitia and AAA). The influence of fiber dispersions upon magnitudes and distributions of stresses and deformations are presented and analyzed. Stresses vary significantly, especially in the AAA case, where material stiffening is significantly influenced by fiber dispersion. The results highlight the need to incorporate the structural differences into finite element simulations to obtain more accurate stress predictions. Additionally, results show the capability of one constitutive model to represent different scenarios of aortic micro-structures allowing future studies of collagen reorientation during disease progression.     

Digital image correlation and finite element modelling as a method to determine mechanical properties of human soft tissue in vivo

The mechanical properties of human soft tissue are crucial for impact biomechanics, rehabilitation engineering, and surgical simulation. Validation of these constitutive models using human data remains challenging and often requires the use of non-invasive imaging and inverse finite element (FE) analysis. Post-processing data from imaging methods such as tagged magnetic resonance imaging (MRI) can be challenging. Digital image correlation (DIC), however, is a relatively straightforward imaging method. DIC has been used in the past to study the planar and superficial properties of soft tissue and excised soft tissue layers. However, DIC has not been used to non-invasive study of the bulk properties of human soft tissue in vivo. Thus, the goal of this study was to assess the use of DIC in combination with FE modelling to determine the bulk material properties of human soft tissue. Indentation experiments were performed on a silicone gel soft tissue phantom. A two camera DIC setup was then used to record the 3D surface deformation. The experiment was then simulated using a FE model. The gel was modelled as Neo-Hookean hyperelastic, and the material parameters were determined by minimising the error between the experimental and FE data. The iterative FE analysis determined material parameters (μ=1.80 kPa, K=2999 kPa) that were in close agreement with parameters derived independently from regression to uniaxial compression tests (μ=1.71 kPa, K=2857 kPa). Furthermore the FE model was capable of reproducing the experimental indentor force as well as the surface deformation found (R²=0.81). It was therefore concluded that a two camera DIC configuration combined with FE modelling can be used to determine the bulk mechanical properties of materials that can be represented using hyperelastic Neo-Hookean constitutive laws.     

Orbital and maxillofacial computer aided surgery: patient-specific finite element models to predict surgical outcomes

This paper addresses an important issue raised for the clinical relevance of Computer-Assisted Surgical applications, namely the methodology used to automatically build patient-specific finite element (FE) models of anatomical structures. From this perspective, a method is proposed, based on a technique called the mesh-matching method, followed by a process that corrects mesh irregularities. The mesh-matching algorithm generates patient-specific volume meshes from an existing generic model. The mesh regularization process is based on the Jacobian matrix transform related to the FE reference element and the current element.This method for generating patient-specific FE models is first applied to computer-assisted maxillofacial surgery, and more precisely, to the FE elastic modelling of patient facial soft tissues. For each patient, the planned bone osteotomies (mandible, maxilla, chin) are used as boundary conditions to deform the FE face model, in order to predict the aesthetic outcome of the surgery. Seven FE patient-specific models were successfully generated by our method. For one patient, the prediction of the FE model is qualitatively compared with the patient's post-operative appearance, measured from a computer tomography scan. Then, our methodology is applied to computer-assisted orbital surgery. It is, therefore, evaluated for the generation of 11 patient-specific FE poroelastic models of the orbital soft tissues. These models are used to predict the consequences of the surgical decompression of the orbit. More precisely, an average law is extrapolated from the simulations carried out for each patient model. This law links the size of the osteotomy (i.e. the surgical gesture) and the backward displacement of the eyeball (the consequence of the surgical gesture).     

Biaxial mechanical response of bioprosthetic heart valve biomaterials to high in-plane shear

Utilization of novel biologically-derived biomaterials in bioprosthetic heart valves (BHV) requires robust constitutive models to predict the mechanical behavior under generalized loading states. Thus, it is necessary to perform rigorous experimentation involving all functional deformations to obtain both the form and material constants of a strain-energy density function. In this study, we generated a comprehensive experimental biaxial mechanical dataset that included high in-plane shear stresses using glutaraldehyde treated bovine pericardium (GLBP) as the representative BHV biomaterial. Compared to our previous study (Sacks, JBME, v.121, pp. 551-555, 1999), GLBP demonstrated a substantially different response under high shear strains. This finding was underscored by the inability of the standard Fung model, applied successfully in our previous GLBP study, to fit the high-shear data. To develop an appropriate constitutive model, we utilized an interpolation technique for the pseudo-elastic response to guide modification of the final model form. An eight parameter modified Fung model utilizing additional quartic terms was developed, which fitted the complete dataset well. Model parameters were also constrained to satisfy physical plausibility of the strain energy function. The results of this study underscore the limited predictive ability of current soft tissue models, and the need to collect experimental data for soft tissue simulations over the complete functional range.     

A nonlocal constitutive model for trabecular bone softening in compression

Using the three-dimensional morphological data provided by computed tomography, finite element (FE) models can be generated and used to compute the stiffness and strength of whole bones. Three-dimensional constitutive laws capturing the main features of bone mechanical behavior can be developed and implemented into FE software to enable simulations on complex bone structures. For this purpose, a constitutive law is proposed, which captures the compressive behavior of trabecular bone as a porous material with accumulation of irreversible strain and loss of stiffness beyond its yield point and softening beyond its ultimate point. To account for these features, a constitutive law based on damage coupled with hardening anisotropic elastoplasticity is formulated using density and fabric-based tensors. To prevent mesh dependence of the solution, a nonlocal averaging technique is adopted. The law has been implemented into a FE software and some simple simulations are first presented to illustrate its behavior. Finally, examples dealing with compression of vertebral bodies clearly show the impact of softening on the localization of the inelastic process.     

Strategies towards rapid generation of forefoot model incorporating realistic geometry of metatarsals encapsulated into lumped soft tissues for personalized finite element analysis

Use of finite element (FE) foot model as a clinical diagnostics tool is likely to improve the specificity of foot injury predictions in the diabetic population. Here we proposed a novel workflow for rapid construction of foot FE model incorporating realistic geometry of metatarsals encapsulated into lumped forefoot’s soft tissues. Custom algorithms were implemented to perform unsupervised segmentation and mesh generation to directly convert CT data into a usable FE model. The automatically generated model provided higher efficiency and comparable numerical accuracy when compared to the model constructed using a traditional solid-based mesh process. The entire procedure uses MATLAB as the main platform, and makes the present approach attractive for creating personalized foot models to be used in clinical studies.     

Modeling Soft Tissue Damage and Failure Using a Combined Particle/Continuum Approach

Biological soft tissues experience damage and failure as a result of injury, disease, or simply age; examples include torn ligaments and arterial dissections. Given the complexity of tissue geometry and material behavior, computational models are often essential for studying both damage and failure. Yet, because of the need to account for discontinuous phenomena such as crazing, tearing, and rupturing, continuum methods are limited. Therefore, we model soft tissue damage and failure using a particle/continuum approach. Specifically, we combine continuum damage theory with Smoothed Particle Hydrodynamics (SPH). Because SPH is a meshless particle method, and particle connectivity is determined solely through a neighbor list, discontinuities can be readily modeled by modifying this list. We show, for the first time, that an anisotropic hyperelastic constitutive model commonly employed for modeling soft tissue can be conveniently implemented within a SPH framework and that SPH results show excellent agreement with analytical solutions for uniaxial and biaxial extension as well as finite element solutions for clamped uniaxial extension in 2D and 3D. We further develop a simple algorithm that automatically detects damaged particles and disconnects the spatial domain along rupture lines in 2D and rupture surfaces in 3D. We demonstrate the utility of this approach by simulating damage and failure under clamped uniaxial extension and in a peeling experiment of virtual soft tissue samples. In conclusion, SPH in combination with continuum damage theory may provide an accurate and efficient framework for modeling damage and failure in soft tissues.