Porohyperelastic finite element modeling of abdominal aortic aneurysms

Abdominal aortic aneurysm (AAA) is the gradual weakening and dilation of the infrarenal aorta. This disease is progressive, asymptomatic, and can eventually lead to rupture--a catastrophic event leading to massive internal bleeding and possibly death. The mechanical environment present in AAA is currently thought to be important in disease initiation, progression, and diagnosis. In this study, we utilize porohyperelastic (PHE) finite element models (FEMs) to investigate how such modeling can be used to better understand the local biomechanical environment in AAA. A 3D hypothetical AAA was constructed with a preferential anterior bulge assuming both the intraluminal thrombus (ILT) and the AAA wall act as porous materials. A parametric study was performed to investigate how physiologically meaningful variations in AAA wall and ILT hydraulic permeabilities affect luminal interstitial fluid velocities and wall stresses within an AAA. A corresponding hyperelastic (HE) simulation was also run in order to be able to compare stress values between PHE and HE simulations. The effect of AAA size on local interstitial fluid velocity was also investigated by simulating maximum diameters (5.5 cm, 4.5 cm, and 3.5 cm) at the baseline values of ILT and AAA wall permeability. Finally, a cyclic PHE simulation was utilized to study the variation in local fluid velocities as a result of a physiologic pulsatile blood pressure. While the ILT hydraulic permeability was found to have minimal affect on interstitial velocities, our simulations demonstrated a 28% increase and a 20% decrease in luminal interstitial fluid velocity as a result of a 1 standard deviation increase and decrease in AAA wall hydraulic permeability, respectively. Peak interstitial velocities in all simulations occurred on the luminal surface adjacent to the region of maximum diameter. These values increased with increasing AAA size. PHE simulations resulted in 19.4%, 40.1%, and 81.0% increases in peak maximum principal wall stresses in comparison to HE simulations for maximum diameters of 35 mm, 45 mm, and 55 mm, respectively. The pulsatile AAA PHE FEM demonstrated a complex interstitial fluid velocity field the direction of which alternated in to and out of the luminal layer of the ILT. The biomechanical environment within both the aneurysmal wall and the ILT is involved in AAA pathogenesis and rupture. Assuming these tissues to be porohyperelastic materials may provide additional insight into the complex solid and fluid forces acting on the cells responsible for aneurysmal remodeling and weakening.     

Identification and determination of material properties for porohyperelastic analysis of large arteries.

A "porohyperelastic" (PHE) material model is described and the theoretical framework presented that allows identification of the necessary material properties functions for soft arterial tissues. A generalized Fung form is proposed for the PHE constitutive law in which the two fundamental Lagrangian material properties are the effective strain energy density function, W(e), and the hydraulic permeability, kij. The PHE model is based on isotropic forms using W(e) = Ue (phi) = 1/2C0(e phi - 1) and the radial component of permeability, kRR = kRR(phi), with phi = C1'(I1 - 3) + C2'(I2 - 3) + K'(J - 1)2. The methods for determination of these material properties are illustrated using experimental data from in situ rabbit aortas. Three experiments are described to determine parameters in Ue and kRR for the intima and media of the aortas, i.e., (1) undrained tests to determine C0, C1', and C2'; (2) drained tests to determine K'; and (3) steady-state pressurization tests of intact and de-endothelialized vessels to determine intimal and medial permeability (adventitia removed in these models). Data-reduction procedures are presented that allow determination of kRR for the intima and media and Ue for the media using experimental data. The effectiveness and accuracy of these procedures are studied using input "data" from finite element models generated with the ABAQUS program. The isotropic theory and data-reduction methods give good approximations for the PHE properties of in situ aortas. These methods can be extended to include arterial tissue remodeling and anisotropic behavior when appropriate experimental data are available.     

Structural models for human spinal motion segments based on a poroelastic view of the intervertebral disk

Analytical and finite element models (FEMs) were used to quantify poroelastic material properties for a human intervertebral disk. An axisymmetric FEM based on a poroelastic view of disk constituents was developed for a representative human spinal motion segment (SMS). Creep and steady-state response predicted by FEMs agreed with experimental observations, i.e., long-time creep occurs with flow in the SMS, whereas for rapid steady-state loading an "undrained," nearly incompressible response is evident. A relatively low value was determined for discal permeability. Transient and long-term creep FE analyses included the study of deformation, pore fluid flow, stress, and pore fluid pressure. Relative fluid motion associated with transient creep is related to nuclear nutrition and the overall mechanical response in the normal disk. Degeneration of the disk may be associated with an increase in permeability.     

A poroviscohyperelastic model for numerical analysis of mechanical behavior of single chondrocyte

The aim of this paper is to use a poroviscohyperelastic (PVHE) model, which is developed based on the porohyperelastic (PHE) model to explore the mechanical deformation properties of single chondrocytes. Both creep and relaxation responses are investigated by using finite element analysis models of micropipette aspiration and atomic force microscopy experiments, respectively. The newly developed PVHE model is compared thoroughly with the standard neo-Hookean solid and PHE models. It has been found that the PVHE can accurately capture both creep and stress relaxation behaviors of chondrocytes better than other two models. Hence, the PVHE is a promising model to investigate mechanical properties of single chondrocytes.     

Analysis of indentation: implications for measuring mechanical properties with atomic force microscopy.

Indentation using the atomic force microscope (AFM) has potential to measure detailed micromechanical properties of soft biological samples. However, interpretation of the results is complicated by the tapered shape of the AFM probe tip, and its small size relative to the depth of indentation. Finite element models (FEMs) were used to examine effects of indentation depth, tip geometry, and material nonlinearity and heterogeneity on the finite indentation response. Widely applied infinitesimal strain models agreed with FEM results for linear elastic materials, but yielded substantial errors in the estimated properties for nonlinear elastic materials. By accounting for the indenter geometry to compute an apparent elastic modulus as a function of indentation depth, nonlinearity and heterogeneity of material properties may be identified. Furthermore, combined finite indentation and biaxial stretch may reveal the specific functional form of the constitutive law--a requirement for quantitative estimates of material constants to be extracted from AFM indentation data.     

Drug transport in artery walls: a sequential porohyperelastic-transport approach

A simulation framework for drug-eluting stents (DES) is presented that simulates the two distinct operational phases of a DES: stent deployment is simulated first, a mechanical porohyperelastic/elasto-plastic/contact analysis. This analysis calculates the interstitial fluid velocity as the result of interstitial fluid pressure gradients and mechanical deformations of the vessel wall. The deformed geometry, interstitial fluid velocity field and porosity field are extracted and used as input for the drug release simulation: a reaction-advection-diffusion (RAD) transport analysis calculating the spatial and temporal drug distribution. The advantage of this approach is that the deformed geometry and interstitial fluid velocity field are not assumed a priori, but are actually calculated using a stent deployment simulation. The framework is demonstrated simulating a DES in an idealised, 3D vessel. Varying mechanical and transport properties based on literature data are assigned to each of the three layers in the wall. The results of the drug release simulation for a period of one week show that the drug distributes longitudinally but will remain in the proximity of the stented area.     

3D finite element model of the diabetic neuropathic foot: A gait analysis driven approach

Diabetic foot is an invalidating complication of diabetes that can lead to foot ulcers. Three-dimensional (3D) finite element analysis (FEA) allows characterizing the loads developed in the different anatomical structures of the foot in dynamic conditions. The aim of this study was to develop a subject specific 3D foot FE model (FEM) of a diabetic neuropathic (DNS) and a healthy (HS) subject, whose subject specificity can be found in term of foot geometry and boundary conditions. Kinematics, kinetics and plantar pressure (PP) data were extracted from the gait analysis trials of the two subjects with this purpose. The FEM were developed segmenting bones, cartilage and skin from MRI and drawing a horizontal plate as ground support. Materials properties were adopted from previous literature. FE simulations were run with the kinematics and kinetics data of four different phases of the stance phase of gait (heel strike, loading response, midstance and push off). FEMs were then driven by group gait data of 10 neuropathic and 10 healthy subjects. Model validation focused on agreement between FEM-simulated and experimental PP.     

Characterizing the interaction among bullet, body armor, and human and surrogate targets

This study used a combined experimental and modeling approach to characterize and quantify the interaction among bullet, body armor, and human surrogate targets during the 10-1000 μs range that is crucial to evaluating the protective effectiveness of body armor against blunt injuries. Ballistic tests incorporating high-speed flash X-ray measurements were performed to acquire the deformations of bullets and body armor samples placed against ballistic clay and gelatin targets with images taken between 10 μs and 1 ms of the initial impact. Finite element models (FEMs) of bullet, armor, and gelatin and clay targets were developed with material parameters selected to best fit model calculations to the test measurements. FEMs of bullet and armor interactions were then assembled with a FEM of a human torso and FEMs of clay and gelatin blocks in the shape of a human torso to examine the effects of target material and geometry on the interaction. Test and simulation results revealed three distinct loading phases during the interaction. In the first phase, the bullet was significantly slowed in about 60 μs as it transferred a major portion of its energy into the body armor. In the second phase, fibers inside the armor were pulled toward the point of impact and kept on absorbing energy until about 100 μs after the initial impact when energy absorption reached its peak. In the third phase, the deformation on the armor's back face continued to grow and energies inside both armor and targets redistributed through wave propagation. The results indicated that armor deformation and energy absorption in the second and third phases were significantly affected by the material properties (density and stiffness) and geometrical characteristics (curvature and gap at the armor-target interface) of the targets. Valid surrogate targets for testing the ballistic resistance of the armor need to account for these factors and produce the same armor deformation and energy absorption as on a human torso until at least about 100 μs (maximum armor energy absorption) or more preferably 300 μs (maximum armor deformation).     

Hydro-mechanical coupling in the periodontal ligament: A porohyperelastic finite element model

Harmonic tension–compression tests at 0.1, 0.5 and 1 Hz on hydrated bovine periodontal ligament (PDL) were numerically simulated. The process was modeled by finite elements (FE) within the framework of poromechanics, with the objective of isolating the contributions of the solid- and fluid phases. The solid matrix was modeled as a porous hyperelastic material (hyperfoam) through which the incompressible fluid filling the pores flowed in accordance with the Darcy’s law. The hydro-mechanical coupling between the porous solid matrix and the fluid phase circulating through it provided an apparent time-dependent response to the PDL, whose rate of deformation depended on the permeability of the porous solid with respect to the interstitial fluid. Since the PDL was subjected to significant deformations, finite strains were taken into account and an exponential dependence of PDL permeability on void ratio – and therefore on the deformation state – was assumed. PDL constitutive parameters were identified by fitting the simulated response to the experimental data for the tests at 1 Hz. The values thus obtained were then used to simulate the tests at 0.1 and 0.5 Hz. The results of the present simulation demonstrate that a porohyperelastic model with variable permeability is able to describe the two main aspects of the PDL’s response: (1) the dependency on strain-rate—the saturated material can develop volumetric strains by only exchanging fluid and (2) the asymmetry between tension and compression, which is due to the effect of both the permeability and the elastic properties on deformation.     

Effect of collagen fibre orientation on intervertebral disc torsion mechanics

The intervertebral disc is a complex fibro-cartilaginous material, consisting of a pressurized nucleus pulposus surrounded by the annulus fibrosus, which has an angle-ply structure. Disc injury and degeneration are noted by significant changes in tissue structure and function, which significantly alters stress distribution and disc joint stiffness. Differences in fibre orientation are thought to contribute to changes in disc torsion mechanics. Therefore, the objective of this study was to evaluate the effect of collagen fibre orientation on internal disc mechanics under compression combined with axial rotation. We developed and validated a finite element model (FEM) to delineate changes in disc mechanics due to fibre orientation from differences in material properties. FEM simulations were performed with fibres oriented at \(\pm 30^{\circ }\) throughout the disc (uniform by region and fibre layer). The initial model was validated by published experimental results for two load conditions, including \(0.48\,\hbox {MPa}\) axial compression and \(10\,\hbox {Nm}\) axial rotation. Once validated, fibre orientation was rotated by \(4^{\circ }\) or \(8^{\circ }\) towards the horizontal plane, resulting in a decrease in disc joint torsional stiffness. Furthermore, we observed that axial rotation caused a sinusoidal change in disc height and radial bulge, which may be beneficial for nutrient transport. In conclusion, including anatomically relevant fibre angles in disc joint FEMs is important for understanding stress distribution throughout the disc and will be important for understanding potential causes for disc injury. Future models will include regional differences in fibre orientation to better represent the fibre architecture of the native disc.     

Numerical simulation of fibrous biomaterials with randomly distributed fiber network structure

This paper presents a computational framework to simulate the mechanical behavior of fibrous biomaterials with randomly distributed fiber networks. A random walk algorithm is implemented to generate the synthetic fiber network in 2D used in simulations. The embedded fiber approach is then adopted to model the fibers as embedded truss elements in the ground matrix, which is essentially equivalent to the affine fiber kinematics. The fiber–matrix interaction is partially considered in the sense that the two material components deform together, but no relative movement is considered. A variational approach is carried out to derive the element residual and stiffness matrices for finite element method (FEM), in which material and geometric nonlinearities are both included. Using a data structure proposed to record the network geometric information, the fiber network is directly incorporated into the FEM simulation without significantly increasing the computational cost. A mesh sensitivity analysis is conducted to show the influence of mesh size on various simulation results. The proposed method can be easily combined with Monte Carlo (MC) simulations to include the influence of the stochastic nature of the network and capture the material behavior in an average sense. The computational framework proposed in this work goes midway between homogenizing the fiber network into the surrounding matrix and accounting for the fully coupled fiber–matrix interaction at the segment length scale, and can be used to study the connection between the microscopic structure and the macro-mechanical behavior of fibrous biomaterials with a reasonable computational cost.     

Solute transport across a contact interface in deformable porous media

A finite element formulation of neutral solute transport across a contact interface between deformable porous media is implemented and validated against analytical solutions. By reducing the integral statements of external virtual work on the two contacting surfaces into a single contact integral, the algorithm automatically enforces continuity of solute molar flux across the contact interface, whereas continuity of the effective solute concentration (a measure of the solute mechano-chemical potential) is achieved using a penalty method. This novel formulation facilitates the analysis of problems in biomechanics where the transport of metabolites across contact interfaces of deformable tissues may be of interest. This contact algorithm is the first to address solute transport across deformable interfaces, and is made available in the public domain, open-source finite element code FEBio (http://www.febio.org).     

Dependence of mechanical behavior of the murine tail disc on regional material properties: a parametric finite element study

In vivo rodent tail models are becoming more widely used for exploring the role of mechanical loading on the initiation and progression of intervertebral disc degeneration. Historically, finite element models (FEMs) have been useful for predicting disc mechanics in humans. However, differences in geometry and tissue properties may limit the predictive utility of these models for rodent discs. Clearly, models that are specific for rodent tail discs and accurately simulate the disc's transient mechanical behavior would serve as important tools for clarifying disc mechanics in these animal models. An FEM was developed based on the structure, geometry, and scale of the mouse tail disc. Importantly, two sources of time-dependent mechanical behavior were incorporated: viscoelasticity of the matrix, and fluid permeation. In addition, a novel strain-dependent swelling pressure was implemented through the introduction of a dilatational stress in nuclear elements. The model was then validated against data from quasi-static tension-compression and compressive creep experiments performed previously using mouse tail discs. Finally, sensitivity analyses were performed in which material parameters of each disc subregion were individually varied. During disc compression, matrix consolidation was observed to occur preferentially at the periphery of the nucleus pulposus. Sensitivity analyses revealed that disc mechanics was greatly influenced by changes in nucleus pulposus material properties, but rather insensitive to variations in any of the endplate properties. Moreover, three key features of the model-nuclear swelling pressure, lamellar collagen viscoelasticity, and interstitial fluid permeation-were found to be critical for accurate simulation of disc mechanics. In particular, collagen viscoelasticity dominated the transient behavior of the disc during the initial 2200 s of creep loading, while fluid permeation governed disc deformation thereafter. The FEM developed in this study exhibited excellent agreement with transient creep behavior of intact mouse tail motion segments. Notably, the model was able to produce spatial variations in nucleus pulposus matrix consolidation that are consistent with previous observations in nuclear cell morphology made in mouse discs using confocal microscopy. Results of this study emphasize the need for including nucleus swelling pressure, collagen viscoelasticity, and fluid permeation when simulating transient changes in matrix and fluid stress/strain. Sensitivity analyses suggest that further characterization of nucleus pulposus material properties should be pursued, due to its significance in steady-state and transient disc mechanical response.     

Load sharing between solid and fluid phases in articular cartilage: II--Comparison of experimental results and u-p finite element predictions.

Experimental measurements in conjunction with theoretical predictions were used to determine the extent of load supported by the fluid phase of cartilage at the articular surface. The u-p finite element model was used to simulate the loading of six separate porcine knee joints and to predict surface deformations of the cartilage layer on the lateral femoral condyle. Representative geometry for the condyle, contact pressures, and intrinsic material properties of the cartilage layer were supplied from experimental measures (see Part I). The u-p finite element predictions for surface deformations of the cartilage layer were obtained for several load partitioning states between the solid and fluid phases of cartilage at the articular surface. These were then compared to actual surface deformations obtained experimentally. It appeared from the comparison that approximately 75 percent of the applied load was borne by the fluid phase at the articular surface under this loading regime. This was qualitatively in agreement with the hypothesis that an applied load to articular joints is partitioned at the surface to the two phases according to the surface area ratios of the solid and fluid phases. It appeared that the solid phase was shielded from the total applied stress on the articular surface by the fluid and could be a reason for the excellent durability of the tissue under the demanding conditions in a diarthrodial joint.     

Finite element algorithm for frictionless contact of porous permeable media under finite deformation and sliding

This study formulates and implements a finite element contact algorithm for solid-fluid (biphasic) mixtures, accommodating both finite deformation and sliding. The finite element source code is made available to the general public. The algorithm uses a penalty method regularized with an augmented Lagrangian method to enforce the continuity of contact traction and normal component of fluid flux across the contact interface. The formulation addresses the need to automatically enforce free-draining conditions outside of the contact interface. The accuracy of the implementation is verified using contact problems, for which exact solutions are obtained by alternative analyses. Illustrations are also provided that demonstrate large deformations and sliding under configurations relevant to biomechanical applications such as articular contact. This study addresses an important computational need in the biomechanics of porous-permeable soft tissues. Placing the source code in the public domain provides a useful resource to the biomechanics community.     

Multiphysics simulation of blood flow and LDL transport in a porohyperelastic arterial wall model

Atherosclerosis localizes at a bend andor bifurcation of an artery, and low density lipoproteins (LDL) accumulate in the intima. Hemodynamic factors are known to affect this localization and LDL accumulation, but the details of the process remain unknown. It is thought that the LDL concentration will be affected by the filtration flow, and that the velocity of this flow will be affected by deformation of the arterial wall. Thus, a coupled model of a blood flow and a deformable arterial wall with filtration flow would be invaluable for simulation of the flow field and concentration field in sequence. However, this type of highly coupled interaction analysis has not yet been attempted. Therefore, we performed a coupled analysis of an artery with multiple bends in sequence. First, based on the theory of porous media, we modeled a deformable arterial wall using a porohyperelastic model (PHEM) that was able to express both the filtration flow and the viscoelastic behavior of the living tissue, and simulated a blood flow field in the arterial lumen, a filtration flow field and a displacement field in the arterial wall using a fluid-structure interaction (FSI) program code by the finite element method (FEM). Next, based on the obtained results, we further simulated LDL transport using a mass transfer analysis code by the FEM. We analyzed the PHEM in comparison with a rigid model. For the blood flow, stagnation was observed downward of the bends. The direction of the filtration flow was only from the lumen to the wall for the rigid model, while filtration flows from both the wall to the lumen and the lumen to the wall were observed for the PHEM. The LDL concentration was high at the lumenwall interface for both the PHEM and rigid model, and reached its maximum value at the stagnation area. For the PHEM, the maximum LDL concentration in the wall in the radial direction was observed at the position of 3% wall thickness from the lumenwall interface, while for the rigid model, it was observed just at the lumenwall interface. In addition, the peak LDL accumulation area of the PHEM moved about according to the pulsatile flow. These results demonstrate that the blood flow, arterial wall deformation, and filtration flow all affect the LDL concentration, and that LDL accumulation is due to stagnation and the presence of filtration flow. Thus, FSI analysis is indispensable.     

Gating mechanisms of mechanosensitive channels of large conductance, I: a continuum mechanics-based hierarchical framework

A hierarchical simulation framework that integrates information from molecular dynamics (MD) simulations into a continuum model is established to study the mechanical response of mechanosensitive channel of large-conductance (MscL) using the finite element method (FEM). The proposed MD-decorated FEM (MDeFEM) approach is used to explore the detailed gating mechanisms of the MscL in Escherichia coli embedded in a palmitoyloleoylphosphatidylethanolamine lipid bilayer. In Part I of this study, the framework of MDeFEM is established. The transmembrane and cytoplasmic helices are taken to be elastic rods, the loops are modeled as springs, and the lipid bilayer is approximated by a three-layer sheet. The mechanical properties of the continuum components, as well as their interactions, are derived from molecular simulations based on atomic force fields. In addition, analytical closed-form continuum model and elastic network model are established to complement the MDeFEM approach and to capture the most essential features of gating. In Part II of this study, the detailed gating mechanisms of E. coli-MscL under various types of loading are presented and compared with experiments, structural model, and all-atom simulations, as well as the analytical models established in Part I. It is envisioned that such a hierarchical multiscale framework will find great value in the study of a variety of biological processes involving complex mechanical deformations such as muscle contraction and mechanotransduction.     

Thermoelastic modelling of the skin at finite deformations

The modelling and computation of the coupled thermal and mechanical response of human skin at finite deformations is considered. The model extends current thermal models to account for thermally- and mechanically-induced deformations. Details of the solution of the highly nonlinear system of governing equations using the finite element method are presented. A representative numerical example illustrates the importance of considering the coupled response for the problem of a rigid, hot indenter in contact with the skin.     

An extensive numerical simulation of the cephalic furrow formation in Drosophila embryo

Here, we propose a new finite element model of the Drosophila embryo to simulate the cephalic furrow (CF) formation. Two different 3D geometries are introduced to reproduce the embryo. The two of them are parametrised through a novel system of curvilinear coordinates, well adapted for biological structures, which is obtained by three fictive electrostatic potentials. A gradient decomposition method is used to take into account both the active and the passive deformations occurring to the cells. Four simulations of the CF are proposed, which in turn considers singular and coupled aspects of the problem.