Modelling Soft Tissue Using Biphasic Theory — A Word of Caution

In recent years the biphasic approach initiated by Mow and coworkers, has been very popular in modelling soft, hydrated, cartilage tissues as well as other soft tissues, such as the brain. This work points out that due to the inherent inability of biphasic models in their present form to account for stress-strain rate dependence resulting from the viscoelasticity of the solid phase, the applicability of these models is limited to the loading conditions producing large relative velocities of phases.     

Strain energy descriptions of biological swelling. I: Single fluid compartment models

Strain energy functions are derived from biphasic soft tissue models in order to describe large-deformation, large-swelling, elastic behavior of nonlinear materials. The resulting analysis leads to calculations of stress-extension relations and tissue fluid pressure. Also explored are the elastic stability of the biphasic tissue models and the manner in which tissue pressure is altered by material deformation.     

An evaluation of three-dimensional diarthrodial joint contact using penetration data and the finite element method

We have developed an approximate method for simulating the three-dimensional contact of soft biphasic tissues in diarthrodial joints under physiological loading. Input to the method includes: (i) kinematic information describing an in vitro joint articulation, measured while the cartilage is deformed under physiological loads, (ii) geometric properties for the relaxed (undeformed) cartilage layers, obtained for the analyses in this study via stereophotogrammetry, and (iii) material parameters for the biphasic constitutive relations used to represent cartilage. Solid models of the relaxed tissue layers are assembled in physiological positions, resulting in a mathematical overlap of the cartilage layers. The overlap distribution is quantified and converted via the biphasic governing equations into applied traction boundary conditions for both the solid and fluid phases for each of the contacting layers. Linear, biphasic, three-dimensional, finite element analysis is performed using the contact boundary conditions derived for each of the contacting layers. The method is found to produce results consistent with the continuity requirements of biphasic contact. Comparison with results from independent, biphasic contact analyses of axisymmetric problems shows that the method slightly underestimates the contact area, leading to an overestimation of the total traction, but yields a good approximation to elastic stress and solid phase displacement.     

A transversely isotropic biphasic model for unconfined compression of growth plate and chondroepiphysis.

Using the biphasic theory for hydrated soft tissues (Mow et al., 1980) and a transversely isotropic elastic model for the solid matrix, an analytical solution is presented for the unconfined compression of cylindrical disks of growth plate tissues compressed between two rigid platens with a frictionless interface. The axisymmetric case where the plane of transverse isotropy is perpendicular to the cylindrical axis is studied, and the stress-relaxation response to imposed step and ramp displacements is solved. This solution is then used to analyze experimental data from unconfined compression stress-relaxation tests performed on specimens from bovine distal ulnar growth plate and chondroepiphysis to determine the biphasic material parameters. The transversely isotropic biphasic model provides an excellent agreement between theory and experimental results, better than was previously achieved with an isotropic model, and can explain the observed experimental behavior in unconfined compression of these tissues.     

Biphasic finite element modeling of hydrated soft tissue contact using an augmented Lagrangian method

A study of biphasic soft tissues contact is fundamental to understanding the biomechanical behavior of human diarthrodial joints. To date, biphasic-biphasic contact has been developed for idealized geometries and not been accessible for more general geometries. In this paper a finite element formulation is developed for contact of biphasic tissues. The augmented Lagrangian method is used to enforce the continuity of contact traction and fluid pressure across the contact interface, and the resulting method is implemented in the commercial software COMSOL Multiphysics. The accuracy of the implementation is verified using 2D axisymmetric problems, including indentation with a flat-ended indenter, indentation with spherical-ended indenter, and contact of glenohumeral cartilage layers. The biphasic finite element contact formulation and its implementation are shown to be robust and able to handle physiologically relevant problems.     

Hydrostatic pressurization and depletion of trapped lubricant pool during creep contact of a rippled indenter against a biphasic articular cartilage layer

This study presents an analysis of the contact of a rippled rigid impermeable indenter against a cartilage layer, which represents a first simulation of the contact of rough cartilage surfaces with lubricant entrapment. Cartilage was modeled with the biphasic theory for hydrated soft tissues, to account for fluid flow into or out of the lubricant pool. The findings of this study demonstrate that under contact creep, the trapped lubricant pool gets depleted within a time period on the order of seconds or minutes as a result of lubricant flow into the articular cartilage. Prior to depletion, hydrostatic fluid load support across the contact interface may be enhanced by the presence of the trapped lubricant pool, depending on the initial geometry of the lubricant pool. According to friction models based on the biphasic nature of the tissue, this enhancement in fluid load support produces a smaller minimum friction coefficient than would otherwise be predicted without a lubricant pool. The results of this study support the hypothesis that trapped lubricant decreases the initial friction coefficient following load application, independently of squeeze-film lubrication effects.     

A penetration-based finite element method for hyperelastic 3D biphasic tissues in contact. Part II: finite element simulations

The penetration method allows for the efficient finite element simulation of contact between soft hydrated biphasic tissues in diarthrodial joints. Efficiency of the method is achieved by separating the intrinsically nonlinear contact problem into a pair of linked biphasic finite element analyses, in which an approximate, spatially and temporally varying contact traction is applied to each of the contacting tissues. In Part I of this study, we extended the penetration method to contact involving nonlinear biphasic tissue layers, and demonstrated how to derive the approximate contact traction boundary conditions. The traction derivation involves time and space dependent natural boundary conditions, and requires special numerical treatment. This paper (Part II) describes how we obtain an efficient nonlinear finite element procedure to solve for the biphasic response of the individual contacting layers. In particular, alternate linearization of the nonlinear weak form, as well as both velocity-pressure, v-p, and displacement-pressure, u-p, mixed formulations are considered. We conclude that the u-p approach, with linearization of both the material law and the deformation gradients, performs best for the problem at hand. The nonlinear biphasic contact solution will be demonstrated for the motion of the glenohumeral joint of the human shoulder joint.     

A finite deformation theory for cartilage and other soft hydrated connective tissues--I. Equilibrium results

The determination of valid stress-strain relations for articular cartilage under finite deformation conditions is a prerequisite for constructing models for synovial joint lubrication. Under physiological conditions of high strain rates and/or high stresses in the joint, large strains occur in cartilage. A finite deformation theory valid for describing cartilage, as well as other soft hydrated connective tissues under large loads, has been developed. This theory is based on the choice of a specific Helmholtz energy function which satisfies the generalized Coleman-Noll (GCN0) condition and the Baker-Ericksen (B-E) inequalities established in finite elasticity theory. In addition, the finite deformation biphasic theory includes the effects of strain-dependent porosity and permeability. These nonlinear effects are essential for properly describing the biomechanical behavior of articular cartilage, even when strain rates are low and strains are infinitesimal. The finite deformation theory describes the large strain behavior of cartilage observed in one-dimensional confined compression experiments at equilibrium, and it reduces to the linear biphasic theory under infinitesimal strain and slow strain rate conditions. Using this theory, we have determined the material coefficients of both human and bovine articular cartilages under large strain conditions at equilibrium. The theory compares very well with experimental results.     

A study of preconditioned Krylov subspace methods with reordering for linear systems from a biphasic v-p finite element formulation

A study was conducted on combinations of preconditioned iterative methods with matrix reordering to solve the linear systems arising from a biphasic velocity-pressure (v-p) finite element formulation used to simulate soft hydrated tissues in the human musculoskeletal system. Krylov subspace methods were tested due to the symmetric indefiniteness of our systems, specifically the generalized minimal residual (GMRES), transpose-free quasi-minimal residual (TFQMR), and biconjugate gradient stabilized (BiCGSTAB) methods. Standard graph reordering techniques were used with incomplete LU (ILU) preconditioning. Performance of the methods was compared on the basis of convergence rate, computing time, and memory requirements. Our results indicate that performance is affected more significantly by the choice of reordering scheme than by the choice of Krylov method. Overall, BiCGSTAB with one-way dissection (OWD) reordering performed best for a test problem representative of a physiological tissue layer. The preferred methods were then used to simulate the contact of the humeral head and glenoid tissue layers in glenohumeral joint of the shoulder, using a penetration-based method to approximate contact. The distribution of pressure and stress fields within the tissues shows significant through-thickness effects and demonstrates the importance of simulating soft hydrated tissues with a biphasic model.     

A finite element analysis of the indentation stress-relaxation response of linear biphasic articular cartilage.

The indentation problem of a thin layer of hydrated soft tissue such as cartilage or meniscus by a circular plane-ended indenter is investigated. The tissue is represented by a biphasic continuum model consisting of a solid phase (collagen and proteoglycan) and a fluid phase (interstitial water). A finite element formulation of the linear biphasic continuum equations is used to solve an axisymmetric approximation of the indentation problem. We consider stress-relaxation problems for which analytic solution is intractable; where the indenter is impermeable (solid) and/or when the interface between the indenter and tissue is perfectly adhesive. Thicknesses corresponding to a thin and thick specimen are considered to examine the effects of tissue thickness. The different flow, pressure, stress and strain fields which are predicted within the tissue, over time periods typically used in the mechanical testing of soft tissues, will be presented. Results are compared with the case of a porous free-draining indenter with a perfectly lubricated tissue-indenter interface, for which an analytic solution is available, to show the effects of friction at the tissue-indenter interface, and the effects of an impermeable indenter. While these effects are present for both thin and thick tissues, they are shown to be more significant for the thin tissue. We also examine the effects of the stiffness of the subchondral bone on the response of the soft tissue and demonstrate that the subchondral bone substrate can be modeled as a rigid, impermeable boundary. The effects of a curved tissue-subchondral bone interface, and the early time response are also studied. For physiologically reasonable levels of curvature, we will show that the curved tissue-subchrondal bone interface has negligible influence on the tissue response away from the interface. In addition, the short-time stress-relaxation responses of the tissue (e.g., at times less than 1s) demonstrate the essential role of the fluid phase in supporting the load applied to the tissue, and by extrapolation to shorter times characteristics of normal joint motion, suggest the essential role of a biphasic model in representing soft tissue behavior in joint response.     

A Study of preconditioned Krylov subspace methods with reordering for linear systems from a biphasic v–p finite element formulation

A study was conducted on combinations of preconditioned iterative methods with matrix reordering to solve the linear systems arising from a biphasic velocity–pressure (v–p) finite element formulation used to simulate soft hydrated tissues in the human musculoskeletal system. Krylov subspace methods were tested due to the symmetric indefiniteness of our systems, specifically the generalized minimal residual (GMRES), transpose-free quasi-minimal residual (TFQMR), and biconjugate gradient stabilized (BiCGSTAB) methods. Standard graph reordering techniques were used with incomplete LU (ILU) preconditioning. Performance of the methods was compared on the basis of convergence rate, computing time, and memory requirements. Our results indicate that performance is affected more significantly by the choice of reordering scheme than by the choice of Krylov method. Overall, BiCGSTAB with one-way dissection (OWD) reordering performed best for a test problem representative of a physiological tissue layer. The preferred methods were then used to simulate the contact of the humeral head and glenoid tissue layers in glenohumeral joint of the shoulder, using a penetration-based method to approximate contact. The distribution of pressure and stress fields within the tissues shows significant through-thickness effects and demonstrates the importance of simulating soft hydrated tissues with a biphasic model.     

Compression testing of very soft biological tissues using semi-confined configuration--a word of caution

We analyse semi-confined (i.e. using no-slip boundary conditions) compression experiment of very soft tissue sample using finite element method. We show that the assumption that the planes perpendicular to the direction of the applied force remain plane during the experiments is not satisfied for compression levels lower than previously stated in Miller [2005. Method for testing very soft biological tissues in compression. Journal of Biomechanics 38, 153-158]. Therefore, we recommend that the parameters for constitutive models of very soft tissues be determined by fitting a solution of the finite element models of the experimental set-up to the measurements obtained using semi-confined compression experiments.     

Updated Lagrangian finite element formulations of various biological soft tissue non-linear material models: a comprehensive procedure and review

Simplified material models are commonly used in computational simulation of biological soft tissue as an approximation of the complicated material response and to minimize computational resources. However, the simulation of complex loadings, such as long-duration tissue swelling, necessitates complex models that are not easy to formulate. This paper strives to offer the updated Lagrangian formulation comprehensive procedure of various non-linear material models for the application of finite element analysis of biological soft tissues including a definition of the Cauchy stress and the spatial tangential stiffness. The relationships between water content, osmotic pressure, ionic concentration and the pore pressure stress of the tissue are discussed with the merits of these models and their applications.     

A comparison between mechano-electrochemical and biphasic swelling theories for soft hydrated tissues

Biological tissues like intervertebral discs and articular cartilage primarily consist of interstitial fluid, collagen fibrils and negatively charged proteoglycans. Due to the fixed charges of the proteoglycans, the total ion concentration inside the tissue is higher than in the surrounding synovial fluid (cation concentration is higher and the anion concentration is lower). This excess of ion particles leads to an osmotic pressure difference, which causes swelling of the tissue. In the last decade several mechano-electrochemical models, which include this mechanism, have been developed. As these models are complex and computationally expensive, it is only possible to analyze geometrically relatively small problems. Furthermore, there is still no commercial finite element tool that includes such a mechano-electrochemical theory. Lanir (Biorheology, 24, pp. 173-187, 1987) hypothesized that electrolyte flux in articular cartilage can be neglected in mechanical studies. Lanir's hypothesis implies that the swelling behavior of cartilage is only determined by deformation of the solid and by fluid flow. Hence, the response could be described by adding a deformation-dependent pressure term to the standard biphasic equations. Based on this theory we developed a biphasic swelling model. The goal of the study was to test Lanir's hypothesis for a range of material properties. We compared the deformation behavior predicted by the biphasic swelling model and a full mechano-electrochemical model for confined compression and 1D swelling. It was shown that, depending on the material properties, the biphasic swelling model behaves largely the same as the mechano-electrochemical model, with regard to stresses and strains in the tissue following either mechanical or chemical perturbations. Hence, the biphasic swelling model could be an alternative for the more complex mechano-electrochemical model, in those cases where the ion flux itself is not the subject of the study. We propose thumbrules to estimate the correlation between the two models for specific problems.