- Research
- Open Access

# Finite element analysis of biological soft tissue surrounded by a deformable membrane that controls transmembrane flow

- Satoko Hirabayashi
^{1}Email authorView ORCID ID profile and - Masami Iwamoto
^{1}

**15**:21

https://doi.org/10.1186/s12976-018-0094-9

© The Author(s) 2018

**Received:**6 May 2018**Accepted:**2 October 2018**Published:**10 December 2018

## Abstract

### Background

Many biological soft tissues are hydrated porous hyperelastic materials, which consist of a complex solid skeleton with fine voids and fluid filling these voids. Mechanical interactions between the solid and the fluid in hydrated porous tissues have been analyzed by finite element methods (FEMs) in which the mixture theory was introduced in various ways. Although most of the tissues are surrounded by deformable membranes that control transmembrane flows, the boundaries of the tissues have been treated as rigid and/or freely permeable in these studies. The purpose of this study was to develop a method for the analysis of hydrated porous hyperelastic tissues surrounded by deformable membranes that control transmembrane flows.

### Results

For this, we developed a new nonlinear finite element formulation of the mixture theory, where the nodal unknowns were the pore water pressure and solid displacement. This method allows the control of the fluid flow rate across the membrane using Neumann boundary condition. Using the method, we conducted a compression test of the hydrated porous hyperelastic tissue, which was surrounded by a flaccid impermeable membrane, and a part of the top surface of this tissue was pushed by a platen. The simulation results showed a stress relaxation phenomenon, resulting from the interaction between the elastic deformation of the tissue, pore water pressure gradient, and the movement of fluid. The results also showed that the fluid trapped by the impermeable membrane led to the swelling of the tissue around the platen.

### Conclusions

These facts suggest that our new method can be effectively used for the analysis of a large deformation of hydrated porous hyperelastic material surrounded by a deformable membrane that controls transmembrane flow, and further investigations may allow more realistic analyses of the biological soft tissues, such as brain edema, brain trauma, the flow of blood and lymph in capillaries and pitting edema.

## Keywords

- Finite element method
- Large deformation
- Pore water pressure
- Mixture theory
- Membrane
- Hydrated poroelastic material

## Background

Finite element method (FEM) is an effective tool for the investigation and prediction of these phenomena. However, the complex micro-level solid-fluid boundaries in the hydrated soft tissues require microscopically fine meshes, making FEM analysis unrealistic. Therefore, FEMs based on the classical consolidation theory (or more rigorous and versatile mixture theory) have been developed. According to the biphasic theory, which is a mixture theory considering a two-phase mixture of solid and fluid phases, the movement of each microscopic component is treated separately. By considering the averaged interactions between them, the mixture is then macroscopically analyzed. Since hydrated biological tissues generally experience large deformations, Oomens et al. [7] developed a nonlinear mixed finite element formulation of mixture theory by using pressure and solid displacement as the nodal unknowns. To assume the solid phase to be hyperelastic, Suh et al. [8] developed a new nonlinear penalty finite element formulation by using solid displacement and fluid velocity as the nodal unknowns. Levenston et al. [9] described two newer three-field (solid displacement, fluid velocity, and pressure) mixed finite element formulations and demonstrated the improved performances of the formulations over an analogous two-field (solid displacement and fluid velocity) penalty formulation. These FEMs have been developed to treat more complicated phenomena [10–13] involving multi-physics phenomena such as charged hydrated tissues [14–19]. In these studies, most of the boundaries of the tissues were treated as rigid and/or freely permeable. Despite the importance of the transmembrane flow control by deformable membrane we mentioned above, these studies didn’t analyze the tissues with the boundaries of deformable membranes that control the transmembrane flows.

## Methods

### Preliminaries of continuum mixture theory

_{x}represents the gradient operator with respect to the current configuration and σ is the Cauchy stress tensor for the mixture (total stress). For a fully saturated mixture with incompressible constituents, the Cauchy stress is

where *p* is the fluid (pore) pressure, I is the rank-two identity tensor, and σ^{E} is the stress induced by solid deformation [7–10]. The stress σ^{E} corresponds to the so-called effective stress described in the classical consolidation theory [7, 10, 20].

*α*∈{

*s*,

*f*} denotes a quantity related to a particular constituent

*α*(

*s*: solid and

*f*: fluid). The original micro-volume, \(\mathrm {d}\Omega _{0}^{\alpha }\), is related to the current volume d

*Ω*

_{t}as:

*J*

^{α}is the determinant of the deformation gradient tensor, F

^{α}(Fig. 4). They are defined as

^{α}deforms to dx

^{α}in the current configuration. Since

*α*is intrinsically incompressible, the true volume of

*α*is constant during the deformation. Then, the current and original volume fractions

*ϕ*

^{α}and \(\phi ^{\alpha }_{0}\) can be defined as:

### Incompressible condition

*α*can be represented as:

*α*. The

*α*-velocity field v

^{α}is defined as:

where x^{α} is the position vector of *α*.

*α*in the volume space

*Ω*

_{t}, denoted by

*m*

^{α}, can be written as:

*α*occupying

*Ω*

_{t}at current time. Therefore, using Eq. (8) and the well-known relationship \(\dot {J}^{\alpha } = J^{\alpha } \boldsymbol {\nabla }_{x} \cdot \boldsymbol {v}^{\alpha }\), the principle of the conservation of mass can be rewritten as:

### Finite element formulation

^{s}, \(\mathrm {d}\Omega _{0}^{s}\),

*J*

^{s}, and F

^{s}are abbreviated to x, d

*Ω*

_{0},

*J*, and F, respectively. The solid displacement u can be defined as

^{s}can be rewritten as

*T*indicates transposition. The first and second Piola-Kirchhoff stresses π and S are related to the total Cauchy stress σ by:

^{E}can be defined for the stress induced by solid deformation using the same relation:

*Γ*

_{0}for the mixture as a whole. The boundary conditions for the solid displacement and the surface traction are prescribed on \(\Gamma _{0}^{u}\) and \(\Gamma _{0}^{t}\), respectively. Here, \(\Gamma _{0}^{u}\) and \(\Gamma _{0}^{t}\) represent the complementary portions of

*Γ*

_{0}:

*Γ*

_{t}:

Generally, the solid and fluid boundary partitions do not need to coincide.

_{χ}represents the gradient operator with respect to the reference configuration. An integral formulation can be obtained by multiplying Eq. (30) by an admissible displacement field \(\boldsymbol {\check {u}}\) and integrating the result over a reference volume

*Ω*

_{0}of the mixture:

where *δ*E is defined consistently with *δ*u, and N is the outward unit normal vector on the reference surface.

*Ω*

_{t}:

where \(\tilde q\) is the given outward flow.

### Linearizations

*n*-th element obtained by subdividing the observed continuum into

*n*

_{el}elements, the displacement of the solid phase, u, can be interpolated as follows:

*n*

_{u}is the total number of nodes for u in the element, \(N_{u}^{(i)}\) is the shape function of

*i*-th node (see Appendix 1 for details), and u

^{(i)}is the

*i*-th nodal displacement. In the same way, the fluid pressure,

*p*can be interpolated as

*n*

_{p}is the total number of nodes for

*p*in the element, \(N_{p}^{(i)}\) is the shape function of

*i*-th node (see Appendix 1 for details), and

*p*

^{(i)}is the

*i*-th nodal pressure. Their virtual values,

*δ*u and

*δ*

*p*, can be interpolated using \(N_{u}^{(i)}\) and \(N_{p}^{(i)}\), respectively:

^{(i)}, E, S, and \(\tilde {\boldsymbol {t}} (\equiv \boldsymbol {\varPi }^{T} \cdot \boldsymbol {N})\) at matrix notation, the following vectors can be defined at an arbitrary point:

*δ*

*u*} and {

*δ*

*E*} in the same way, matrix [

*B*] can be obtained as a matrix that satisfies the following relationship:

*N*

_{u}] as

where \(\Omega _{0}^{e}\) and \(\Gamma _{0}^{e}\) denote the parts of *Ω*_{0} and \(\Gamma _{0}^{t}\) included in the *n*-th element.

*δ*

*p*} is defined in the same way. The 1×

*n*

_{p}matrix [

*N*

_{p}] can be defined as:

*z*] is defined as:

*κ*] as the coefficient matrix of

*κ*, Eq. (36) leads to the following equation:

where \(\Omega _{t}^{e}\) and \(\Gamma _{t}^{e}\) denote the parts of *Ω*_{t} and \(\Gamma _{t}^{w}\) included in the *n*-th element.

*F*

_{ext}}, {

*F*

_{int}}, {Λ

_{ext}}, and {Λ

_{int}} represent external force, internal force, external outward flow, and internal outward flow, respectively. Here, the external outward flow {Λ

_{ext}} is the value provided as the Neumann boundary condition representing the characteristics of the membrane which surrounded the mixture. For spatial integration, the Gauss integration can be used. As \(\{\delta \bar u\}\) and \(\{\delta \bar p\}\) are arbitrary vectors,

*Δ*

*t*. By using the backwards Euler algorithm, \(\dot {\boldsymbol {u}}\) at

*t*+

*Δ*

*t*can be replaced:

*t*are known and represented as u

_{t}and

*p*

_{t}, they are adopted as the initial predictors u

^{(0)}and

*p*

^{(0)}for the iterative scheme to obtain the solution at time

*t*+

*Δ*

*t*:

*k*-th predictors u

^{(k)}and

*p*

^{(k)}are

*Δ*u

^{(k)}and

*Δ*

*p*

^{(k)}are the

*k*-th incremental values. By arranging

*Δ*u

^{(k)}and

*Δ*

*p*

^{(k)}at all nodes in a row as the vectors \(\{\Delta \bar {u}\}^{(k)}\) and \(\{\Delta \bar {p}\}^{(k)}\), respectively Eq. (54) can be linearized as:

*F*

_{int}}

^{(k−1)}and {Λ

_{int}}

^{(k−1)}are {

*F*

_{int}} and {Λ

_{int}}, respectively, calculated using the predictors obtained at the previous iteration, u

^{(k−1)}and

*p*

^{(k−1)}. For the

*k*-th iteration, the system of equations that is obtained by coupling these two systems of equations is solved for \(\{\Delta \bar {u}\}^{(k)}\) and \(\{\Delta \bar {p}\}^{(k)}\):

*t*+

*Δ*

*t*. For this present analysis, the following convergence criterion of the norm of the incremental vector was used:

where *TOL* is a user-specified tolerance with the value of 10^{−5} in this study.

## Results

*W*:

^{E}and E at matrix notation were used. We employed the following hyperelastic energy function:

*λ*and

*μ*represent Lamé’s constants. Here, S

^{E}is linear to E:

The hyperelastic properties dominate the solid displacements, fluid velocities, and the pore water pressures throughout the deformation process. Although the dense solid constituting the skeleton is incompressible, the skeleton is compressible because it is porous and contains compressible voids. Thus, we assumed that Poisson’s ratio of the skeleton is 0.0. We also assumed that the skeleton is homogeneous with Young’s modulus = 0.3 kPa (i.e., *λ*=0.0 and *μ*=0.15 kPa), an initial solid volume fraction \(\phi ^{s}_{0}=\)0.2, and a permeability κ=*κ*I (*κ*=8.0 m^{4}/N-s). As indicated in Eqs. (35) and (36), the permeability κ affects the fluid and solid velocities, w and \(\dot {\boldsymbol {u}}\), resulting in a viscous behavior. With smaller κ, the fluid velocity inside the mixture becomes lower. Except specific conditions, the lower fluid velocity makes solid velocity lower. In general, this kind of viscous behavior may absorb a shock of impact by reducing accelerations of objects. The viscous behavior can protect the brain from the hit of skull in the SAS when considering the trabecular network filled with CSF as the hydrated soft tissue.

_{ext}}=0, which has not ever been represented using the FEMs formulated by Suh et al. [8] and Levenston et al. [9] by the reason mentioned previously (Fig. 6). The permeability of the membrane, which is given as the boundary condition in this study, changes the total volume of the mixture by determining the transmembrane flow. It may result in an edema. Although the membrane is assumed to be flaccid, it cannot deform when the top surface is uniformly loaded because the tissue cannot deform due to its incompressibility. Thus, we conducted a three-dimensional test, where the platen covered only the central part of the top surface. The platen was compressed at the velocity of 6 m/s during 1 ms of the compression period, and maintained at that compression level, 6 mm, during the following hold period. The finite element meshes used in the analyses represent 1/4 of the tissue (Fig. 7). The platen axially pushed the areas surrounded by bold lines. At the other area of the top surface, the load perpendicular to the surface was kept at 0, i.e., π

^{T}·N=0. The displacements perpendicular to the surface were confined at the sides and base. We prepared two types of meshes, Mesh A and Mesh B, which consisted of 1200 and 4800 5/4c displacement-pressure elements, respectively. Three different sets of time steps and mesh types used for the numerical examples are presented in Table 1.

Time stepts and mesh types used for the numerical examples

Parameter | Case 1 | Case 2 | Case 3 |
---|---|---|---|

| 0.1 ms | 0.2 ms | 0.2 ms |

Mesh type | Mesh A | Mesh A | Mesh B |

*ϕ*

^{s}increased in the area under the platen due to the outflow of fluid, and decreased under the surface around the platen because of the fluid block and accumulation by the impermeable membrane. The accumulated fluid led to the swelling of the tissue around the platen at the end of the compression period (1 ms). During the hold period, the movement of the fluid relaxed the pore water pressure gradient (Fig. 9), which resulted in the reduction of the reaction force (Fig. 8) and the fluid velocity (Fig. 10). Afterwards, the reaction force converged to the value derived from the elastic deformation of the solid phase, showing a typical stress relaxation (Fig. 8) until the fluid stopped. The similar behaviors were observed in other two analyses.

The total volume decreased to 93.3% at the end of the compression period (Fig. 8), against the expectation that the total volume is maintained due to the impermeable membrane covering the entire tissue. Additionally, the surface around the platen sank during the hold period, while the total volume remained stable and no counter flow occurred, i.e., the fluid flow to compressed area from the surrounding areas.

## Discussion

Here, based on the mixture theory, we performed the finite element formulation of a hyperelastic porous solid filled with fluid, which flowed according to the pore water pressure gradient and interacted with the deformation of the solid. In the method developed and presented here, only the solid displacements and the pore water pressures are treated as the nodal unknowns and the fluid velocities are not regarded as such. Therefore, the fluid flow volume across the membrane can be directly given as the Neumann boundary condition. Our new method enables us to analyze a large deformation of hydrated tissue surrounded by a deformable membrane that controls transmembrane flow. While Oomens et al. [7], who also used the set of the solid displacements and the pore water pressures for the nodal unknowns, employed a rate-type constitutive equation for the solid phase, we assumed that the drained solid skeleton is a hyperelastic material that has elastic potential function.

Using our new method, we analyzed the hydrated soft tissues surrounded by flaccid impermeable membranes. Our new method has the following three limitations based on the simulation results. Firstly, we could neither validate our proposed method using a simple problem with an analytical solution, such as one-dimensional compression, nor compare our results with the results obtained by previously presented numerical methods since there are no studies on numerical methods of the hydrated biological soft tissues surrounded by deformable membranes that control the transmembrane flows. Further studies are needed for validation of our proposed study using some sort of analytical solutions or experimental data. Secondly, although the whole tissue was surrounded by the fully impermeable membrane, the total volume decreased to 93.3% by the end of the compression period, showing that the results contain numerical errors. The errors should be investigated and corrected in the future. Thirdly, the peak values of the reaction force were shown to depend on the mesh. This is probably because too coarse mesh failed to capture the nonlinear pressure gradient around the membrane, generating the inconsistency of the reaction force. Further investigation is necessary for the problem of mesh dependency.

Despite above mentioned limitations of our new method, the fluid extruded by platen was trapped by the impermeable membrane and it led to the swelling of the tissue around the platen during the compression period, while the total volume was well maintained during the hold period, suggesting the ability of our method to reproduce a deformable membrane that controls transmembrane flow.

Fluid flowed according to the pore water pressure gradient caused by the compression, resulting in the relaxation of the gradient. Afterwards, the reaction force and the fluid velocity were reduced and converged to the constant values. As a result, the reaction force showed a typical stress relaxation although the constitutive equation of the solid skeleton did not contain the viscoelastic component and the permeability was constant. This viscoelastic behavior can play important roles in the body. For example, the viscoelasticity of SAS can protect the brain by absorbing shocks. As the behavior involves the volume changes of elements with the total volume preserved by the interactions between elements via the movement of the fluid, it is impossible to reproduce this viscoelastic behavior by substituting an equivalent viscoelastic body for the mixture material used in this study.

## Conclusions

To analyze the hydrated soft tissues surrounded by deformable membranes that control transmembrane flows, which are observed in most of the tissues, we developed a new nonlinear FEM based on mixture theory using the pore water pressure and solid displacement as the nodal unknowns, and controlled the transmembrane flows across a deformable membrane with the Neumann boundary condition. The results obtained by our new method demonstrated that our new method can effectively analyze large deformation of hydrated tissues surrounded by a flaccid impermeable membrane. Although the proposed method has some limitations on the validation, numerical errors, and mesh dependency, it is confirmed that the proposed method can reproduce the viscoelastic behavior characterized in hydrated biological soft tissues that the fluid trapped by the impermeable membrane led to the swelling of the tissue around the platen. Further developments of this method may allow the analysis of the biological phenomena such as brain edema, brain trauma, the flow of blood and lymph in capillaries and pitting edema.

## Appendix 1

### Piola transformation

*Γ*

_{t}on the current configuration, the Nanson’s formula can be denoted as:

*Γ*

_{0}on the reference configuration. Using these equations, the following relationship is obtained:

## Appendix 2

### Shape functions in a 5/4c displacement-pressure element

As the continuity Eq. (36) involved the second derivative of the pressure field, *C*^{0} continuity was required for the pressure field. In addition, the mixed formulation requires a suitable combination of the displacement and pressure interpolations in order to obtain a robust element. In the present analysis, 5/4c (*P*1 + /*P*1) displacement-pressure element, which is in the family of MINI elements, is used. This is a tetrahedral first-order element with a linear continuous interpolation. The pressure can be interpolated by 4 nodes located at vertices. For the interpolation of the displacement, a bubble node located at the centroid was added as the 5th node.

*r*

_{i},(

*i*=1 to 4), the volume coordinate of the P can be defined as (

*r*

_{1},

*r*

_{2},

*r*

_{3}). By using the volume coordinate, the shape functions of the displacement and the pressure are defined as

## Declarations

### Acknowledgements

We thank Prof. Toshiaki Hisada, of the University of Tokyo, for the years of providing advice. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. We also thank Editage (www.editage.jp) for English language editing.

### Availability of data and materials

The datasets generated during the current study are available from the corresponding author on reasonable request.

### Authors’ contributions

SH made substantial contributions to conception and design, analysis and interpretation of data. MI was involved in drafting the manuscript and agreed to be accountable for all aspects of the work. Both authors read and approved the final manuscript.

### Ethics approval and consent to participate

Not applicable.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interest.

### Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## Authors’ Affiliations

## References

- Yamashima T, Friede RL. Light and electron microscopic studies on the subdural space, the subarachnoid space and the arachnoid membrane. Neurol Med Chir. 1984; 24(10):737–46.View ArticleGoogle Scholar
- Killer H, Laeng H, Flammer J, Groscurth P. Architecture of arachnoid trabeculae, pillars, and septa in the subarachnoid space of the human optic nerve: anatomy and clinical considerations. Br J Ophthalmol. 2003; 87(6):777–81.View ArticleGoogle Scholar
- Gupta S, Soellinger M, Boesiger P, Poulikakos D, Kurtcuoglu V. Three-dimensional computational modeling of subject-specific cerebrospinal fluid flow in the subarachnoid space. J Biomech Eng. 2008; 131(2):021010.View ArticleGoogle Scholar
- Saboori P, Sadegh A. Material modeling of the head’s subarachnoid space. Sci Iran. 2011; 18(6):1492–9.View ArticleGoogle Scholar
- Saboori P, Sadegh A. Histology and morphology of the brain subarachnoid trabeculae. Anat Res Int. 2015; 2015:9. Article ID 279814, https://doi.org/10.1155/2015/279814.Google Scholar
- Raslan A, Bhardwaj A. Medical management of cerebral edema. Neurosurg Focus. 2007; 22(5):1–12.View ArticleGoogle Scholar
- Oomens CWJ, van Campen DH, Grootenboer HJ. A mixture approach to the mechanics of skin. J Biomech. 1987; 20(9):877–85.View ArticleGoogle Scholar
- Suh JK, Spilker R, Holmes M. A penalty finite element analysis for nonlinear mechanics of biphasic hydrated soft tissue under large deformation. Int J Numer Methods Eng. 1991; 32(7):1411–39.View ArticleGoogle Scholar
- Levenston ME, Frank EH, Grodzinsky AJ. Variationally derived 3-field finite element formulations for quasistatic poroelastic analysis of hydrated biological tissues. Comput Methods Appl Mech Eng. 1998; 156(1):231–46.View ArticleGoogle Scholar
- Stastna M, Tenti G, Sivaloganathan S, Drake JM. Brain biomechanics: consolidation theory of hydrocephalus. Variable permeability and transient effects. Can Appl Math Q. 1999; 7:111–24.Google Scholar
- Taylor Z, Miller K. Reassessment of brain elasticity for analysis of biomechanisms of hydrocephalus. J Biomech. 2004; 37(8):1263–9.View ArticleGoogle Scholar
- Franceschini G, Bigoni D, Regitnig P, Holzapfel GA. Brain tissue deforms similarly to filled elastomers and follows consolidation theory. J Mech Phys Solids. 2006; 54(12):2592–620.View ArticleGoogle Scholar
- Cheng S, Bilston LE. Unconfined compression of white matter. J Biomech. 2007; 40(1):117–24.View ArticleGoogle Scholar
- Chen Y, Chen X, Hisada T. Non-linear finite element analysis of mechanical electrochemical phenomena in hydrated soft tissues based on triphasic theory. Int J Numer Methods Eng. 2006; 65(2):147–73.View ArticleGoogle Scholar
- Lu XL, Wan LQ, Guo XE, Mow VC. A linearized formulation of triphasic mixture theory for articular cartilage, and its application to indentation analysis. J Biomech. 2010; 43(4):673–9.View ArticleGoogle Scholar
- Lai WM, Hou JS, Mow VC. A triphasic theory for the swelling and deformation behaviors of articular cartilage. J Biomech Eng. 1991; 113(3):245–58.View ArticleGoogle Scholar
- Sun DN, Gu WY, Guo XE, Lai WM, Mow VC. A mixed finite element formulation of triphasic mechano-electrochemical theory for charged, hydrated biological soft tissues. Int J Numer Methods Eng. 1999; 45(10):1375–402.View ArticleGoogle Scholar
- Huyghe JM, Wilson W, Malakpoor K. On the thermodynamical admissibility of the triphasic theory of charged hydrated tissues. J Biomech Eng. 2009; 131(4):044504.View ArticleGoogle Scholar
- Ateshian GA, Chahine NO, Basalo IM, Hung CT. The correspondence between equilibrium biphasic and triphasic material properties in mixture models of articular cartilage. J Biomech. 2004; 37(3):391–400.View ArticleGoogle Scholar
- Biot MA. Mechanics of Deformation and Acoustic Propagation in Porous Media. J Appl Phys. 1962; 33(4):1482–98.View ArticleGoogle Scholar
- Hughes TJR, Marsden JE. Some applications of geometry is continuum mechanics. Rep Math Phys. 1977; 12(1):35–44.View ArticleGoogle Scholar