How Compressibility Influences the Mechanical Bidomain Model

Compressibility influences the mechanical bidomain model describing the elastic properties of tissue. Displacements of the intracellular and extracellular spaces are analyzed individually, and differences in these displacements produce membrane forces. Two length constants are associated with the membrane spring constant, one contains the shear moduli and the other contains the bulk moduli. The analytical solutions in these examples indicate that the monodomain part does not contribute to the membrane force. Accounting for compressibility has its largest impact on the intracellular and extracellular pressures. The bidomain contribution to the pressure obeys the Helmholtz equation rather than Laplace’s equation. This model predicts membrane forces that might cause tissue remodeling or mechanotransduction. Keywords-biomechanics; mechanical bidomain model; mechanotransduction; pressure; remodeling.


I. INTRODUCTION
The mechanical bidomain model is a new mathematical description of the biomechanics of tissue, which distinguishes between displacements in the intracellular and extracellular spaces and focuses on forces across the cell membrane [1][2][3][4][5].Such membrane forces may affect transmembrane proteins, such as integrins, that are responsible for tissue remodeling and mechanotransduction [6][7][8].Membrane forces may also play a role in tissue engineering because mechanical stresses influence the growth of replacement tissue [9], in remodeling of blood vessels [10], and in development because mechanical stresses guide the growth of fetal tissue [11].The unique feature of the mechanical bidomain model is that it is macroscopic (representing the tissue rather than individual cells) but because it tracks the intracellular and extracellular spaces individually it can predict membrane forces.The mechanical bidomain model is analogous to the electrical biodmain model [12], which is currently the most widely used model for simulating defibrillation.For a review of the mechanical bidomain model, see [4].
In a previous version of the mechanical bidomain model [4], both the intracellular and extracellular spaces were incompressible, so tissue displacements did not change the tissue volume (dilatation).The pressure distribution in a tissue is defined as the bulk modulus times the dilatation [13,14].The incompressible limit corresponds to the finite product of these quantities as the bulk modulus goes to infinity and the dilatation goes to zero.Tissue is largely water, so the incompressible assumption is often valid.Tissue bulk moduli are much larger than shear moduli.The bulk modulus is on the order of 2 × 10 9 Pa, whereas the shear modulus is about 2 × 10 4 Pa, a difference of a factor of 100,000 [15], so one can typically assume incompressibility.
The mechanical bidomain model is complicated because there are two bulk moduli: one for the intracellular space, χ, and one for the extracellular space, δ.Although one can easily imagine a limiting process to arrive at intracellular and extracellular pressures, p and q, similar to that outlined above, the ratio χ δ may impact the tissue behavior, even if χ and δ are both much greater than the intracellular and extracellular shear moduli, ν and µ.Furthermore, the pressures contribute to the boundary conditions at the tissue surface, and the appropriate boundary conditions remain uncertain.
The goal of this paper is to resolve questions surrounding the mechanical bidomain model by rederiving the model without the assumption of incompressibility.We solve several biomechanics problems, analyze their solutions, and then impose incompressibility to determine the correct behavior in that limit.This procedure introduces new qualitative behavior into the model.For instance, previous studies revealed a bidomain length constant that depended on the membrane spring constant K coupling the intracellular and extracellular spaces [1,4].In this paper, there are two length constants, both involving K; one containing the shear moduli and one the bulk moduli.Our results change the way we calculate and interpret the intracellular and extracellular pressures, impact the predicted displacements, and affect the membrane force distribution.

A. Mechanical Bidomain Model
Consider an isotropic tissue and ignore any active tension.To keep things simple, we consider a two-dimensional Cartesian coordinate system (x, (3) The relationships between intracellular and extracellular stresses and strains are τ eyy = δ( exx + eyy ) + 2µ eyy (7) τ ixy = 2ν ixy (8) The equations of mechanical equilibrium are Plugging (1-9) into (10-13) and rearranging, we obtain the equations governing the mechanical bidomain model Fig. 1 shows a schematic drawing of a twodimensional sheet of tissue representing the mechanical bidomain model.The intracellular and extracellular spaces are each represented by a grid of springs (green and blue), and their behavior is described by the left-hand sides of (14-17).These grids are connected across the cell membrane by an array of springs (red) coupling the two spaces.The force produced by these springs is represented by the right-hand sides of (14)(15)(16)(17).The membrane force depends on the difference between the intracellular and extracellular displacements and the spring constant K.

B. Shear Displacement
We first analyze shear displacements using methods similar to those presented previously [4].Take the y-derivative of ( 14) and subtract the x-derivative of ( 16), and similarly take yderivative of (15) and subtract the x-derivative of (17).The results are simpler when expressed in terms of the z-components of the curl of the displacements, ω i = (∇ × u) z = ∂uy ∂x − ∂ux ∂y and ω e = (∇ × w) z = ∂wy ∂x − ∂wx ∂y , where ν and µ are the shear moduli (terms containing the bulk moduli χ and δ cancel out).If we add (18) and (19), we get If we divide (18) by ν and (19) by µ and subtract, we obtain where σ = νµ K(ν+µ) .
Define α = ω i + µ ν ω e and β = ω i − ω e , so ω i = ν ν+µ (α + µ ν β) and ω e = ν ν+µ (α − β).These functions obey ∇ 2 α = 0 and ∇ 2 β = 1 σ 2 β.The function α is a weighted sum of the intracellular and extracellular spaces, obeys Laplace's equation, and describes the "monodomain" behavior of the system; β is the difference between the intracellular and extracellular spaces, obeys the Helmholtz equation so it will decay with distance by the length constant σ, and specifies the "bidomain" behavior.Because β represents the difference between u and w, it is responsible for the membrane force.

C. Volume Dilatation
We next consider volume changes in the tissue, which is the new feature presented in this paper.Take the x-derivative of ( 14) and add it to the yderivative of (16), and similarly take x-derivative of (15) and add it to the y-derivative of (17), and get where e i = ∇•u and e e = ∇•w are the intracellular and extracellular dilatations.Now assume that the bulk moduli are much greater than the shear moduli (χ, δ >> ν, µ) so δ∇ 2 e e = −K(e i − e e ).
If we add the two equations, we get If we divide (24) by χ and (25) by δ and subtract, we obtain where Define the intracellular and extracellular pressures as p = −χe i and q = −δe e .To examine the incompressible limit, let χ and δ go to infinity and e i and e e go to zero in such a way that p and q remain finite.Furthermore, define two auxiliary pressures: P = p+q is the monodomain pressure, and Q = p − χ δ q is a weighted difference of the pressures.The pressures p and q are expressed in terms of P and Q as The pressures P and Q obey the equations ∇ 2 P = 0, ∇ 2 Q = 1 ξ 2 Q.P obeys Laplaces equation and represents the monodomain behavior.Q obeys the Helmholtz equation, decays with distance by a length constant ξ, and represents the bidomain behavior.

D. Summary of the Methods
This analysis shows that the displacements caused by shearing and the dilatations caused by volume changes obey analogous equations.Both can be separated into a monodomain part and a bidomain part.Both introduce new length constants, σ and ξ, that characterize the behavior.The analysis of shear was presented previously [4] and results in a boundary layer near the tissue edge [5].The analysis of volume dilatation is new: it represents an advance beyond the previous model.
The equations above suggest that the pressure and displacement are independent.However, they are in fact coupled by the boundary conditions at the tissue surface.For instance, if the tissue is perfused by an adjacent bath, the boundary conditions are: the normal component of the extracellular stress is continuous with the stress in the bath, and the normal component of the intracellular stress is zero.If the tissue-bath surface were specified by the plane y = constant, then the boundary conditions would be τ exx = −P bath and τ exy = τ ixx = τ ixy = 0. Here, P bath is the hydrostatic pressure in the bath, which we take as a static fluid.
Our next goal is to examine three examples, from the very simple to the increasingly complex, that illustrate the behavior of the mechanical bidomain model, with our focus primarily on the dilatation and pressure.

A. Example 1: Cylinder in a Bath
Consider a cylinder of isotropic tissue (of radius b) immersed in a fluid bath at pressure P 0 (Fig. 2).This might represent, for example, a cardiac papillary muscle lying in a superfusing bath.We examine this case first because shear forces play no role in the behavior; we need to determine only the pressure distribution in the tissue.In this case, the displacement is in the radial direction r and does not depend on the angle θ, so in cylindrical coordinates the strains are related to the displacement by irr = ∂ur ∂r , iθθ = ur r , irθ = 0, err = ∂wr ∂r , eθθ = wr r , and erθ = 0 [16].The stresses and strains are related by τ irr = χ( irr + iθθ )+2ν irr , τ iθθ = χ( irr + iθθ )+2ν iθθ , τ irθ = 0, τ err = δ( err + eθθ ) + 2µ err , τ eθθ = δ( err + eθθ ) + 2µ eθθ , and τ erθ = 0.
The equations of mechanical equilibrium written in cylindrical coordinates are [16] At the surface of the tissue (r = b), the boundary conditions are τ irr = 0 and τ err = −P 0 .The displacements that obey the equations of mechanical equilibrium and boundary conditions are where I 0 and I 1 are modified Bessel functions of the first kind.Membrane forces arise due to the difference in displacements, uw.In this case, u r − w r = P 0 The membrane force is largest near the boundary r = b.However, the displacements go to zero as χ and δ go to infinity, so the membrane force vanishes for an incompressible tissue.
The pressures p = −χe i and q = −δe e remain finite when χ and δ go to infinity Fig. 3 shows p and q as functions of r.These pressures can be recast in terms of P and Q, P = P 0 and The monodomain pressure P is just the bath pressure, while the bidomain pressure Q decays with length constant ξ as one moves inward from the tissue surface.If ξ << b, then several length constants below the surface Q is approximately zero and the intracellular and extracellular pressures are constants; p = χ χ+δ P 0 and q = δ χ+δ P 0 .Although the pressures remain finite as the bulk moduli go to infinity, they still depend on the ratio of the bulk moduli, χ δ .The tissue behavior is independent of the shear moduli ν and µ.
The length constant ξ depends on the bulk moduli, and as χ and δ go to infinity ξ grows.However, the value of the membrane spring constant K is not known but is expected to be large [4].Therefore, the size of ξ = χδ K(χ+δ) relative to the cylinder radius b is not obvious.If ξ << b, the pressure changes over a thin boundary layer near the tissue surface, like shown in Fig. 3.If ξ >> b, p = 0 and q = P 0 , and are constant.The limit ξ >> b corresponds to the prediction of the mechanical bidomain model presented previously [4].Ultimately, the size of ξ is an experimental question that cannot be resolved until we have the necessary data.One virtue of this first example is that it highlights a key new feature of our model, the new length constant ξ, and indicates what quantities need to be measured to assess the model.

B. Example 2: Blood Vessel
The next example is of a cylinder of fluid, of radius a, surrounded by tissue, which is a model for a blood vessel (Fig. 4).This example is slightly more complicated than the first example, because the shear moduli now play a role in the solution.However, the approach is similar to Example 1.The displacements that obey the equations of equi- librium and the boundary conditions are where K 0 and K 1 are modified Bessel functions of the second kind.These expressions each have two terms: a first term, P0 ν+µ a 2 2r , that governs the monodomain behavior and does not go to zero as the bulk moduli go to infinity, and a second bidomain term containing Bessel functions.Because the monodomain terms are the same in the intracellular and extracellular spaces, they contribute nothing to the membrane force.The monodomain part (Fig. 5a) implies that a high pressure inside the vessel causes it to expand.The difference between the intracellular and extracellular displacements is δ , which gives rise to membrane forces near the tissue boundary r = a that fall off with a length constant ξ (Fig. 5b).The difference u r − w r vanishes when χ and δ go to infinity.Therefore in the incompressible limit, there is no membrane force.
The dilatations e i = P0  and c) the pressures p and q, all as functions of r, for ν = 2µ , χ = 2δ, and ξ = a 10 .The displacements and pressures have been normalized to their values at r = a.q = P 0 ν ν+µ K0(r/ξ) K0(a/ξ) are independent of χ and δ except though their dependence on ξ (Fig. 5c).The monodomain pressure is zero, P = 0, so there is only a bidomain contribution.In the limit when ξ >> a (as presented by [4]), the pressures become constant: p = −P 0 ν ν+µ and q = P 0 ν ν+µ .This example emphasizes two points.First, like Example 1, it illustrates how the pressures fall off from the boundary with length constant ξ, the new parameter introduced in this paper.Second, it shows how the displacement is divided into two parts, one of which gives rise to membrane forces (although in this case the membrane force goes to zero as χ and δ go to infinity).

C. Example 3: Sheet of Active Tissue
The last example is of a two-dimensional sheet of cardiac tissue undergoing an active contraction (Fig. 6).This example is more complicated than the first two, but we include here it for several reasons.First, it shows how the two length constants, σ and ξ, can both contribute to the solution.Second, it contains a non-zero membrane force even when the bulk moduli go to infinity.Third, it was analyzed in detail using the incompressible model [5], so redoing the calculation with compressibility highlights those changes compressibility introduces.Finally, the example shows how the active tension enters the equations, which is a crucial element when modeling cardiac tissue.
We represent the active contraction as a uniform tension T added to the intracellular stress tensor [13,14].The tension acts along the myocardial fibers, which we assume are straight, uniform and oriented in the x direction.Other than this active tension, we assume that the tissue is isotropic.In Cartesian coordinates, the tension would be represented by a constant T added to τ ixx .In cylindrical coordinates, the tension enters the stress tensor in a more complicated way, shown below.
In cylindrical coordinates, the strains are expressed in terms of the displacements by [16] irr = ∂u r ∂r err = ∂w r ∂r (39 The stress and strain are related by [5] τ irr = χ( irr + iθθ )+2ν irr + T 2 (1+cos 2θ) (43) The equations of equilibrium are [16] ∂τ irr ∂r + 1 r At the surface of the tissue r = R, the boundary is stress free: τ irr = τ err = τ irθ = τ erθ = 0.The displacements obeying these equations of equilibrium and boundary conditions are where Both u and w contain a leading monodomain term that is in general larger than the subsequent bidomain terms (Fig. 7a).The tissue contracts along the fiber direction, and expands perpendicular to the fiber direction.(Note: Fig. 7a appears different than Fig. 1 in [5] because Fig. 1 in [5] is incorrect.It should look exactly like our Fig.7a).As χ and δ go to infinity, the difference in the displacements is In general, σ << R. The membrane force is largest within a few length constants σ of the edge r = R, and the θ component is larger than the r component (Fig. 7b).(Note: our Fig.7b is not identical to Fig. 2 in [5] because we use a different value of σ R ; [5] used σ = R 10 , whereas we use σ = R 100 .) Taking the limit as χ and δ go to infinity and e i and e e go to zero, the pressures become The pressures can be expressed as P = T 2 and Fig. 8 shows p and q as functions of position.These pressure distributions are very different than those shown in Figs. 3 and 4 of [5].Away from the tissue edge the pressures in our Fig. 8 are both constant, with deviations from this constant value only within a few length constants of the boundary.

IV. DISCUSSION
The analysis in the Methods and Results illustrates how compressibility affects the mechanical bidomain model.It is useful to compare the previous derivation of the model with the version developed in this paper.If we assume an isotopic tissue T = 0 and analyze the model developed previously [4] using the methods derived here, the mechanical bidomain equations would be ∇ 2 α = 0, ∇ 2 β = 1 σ 2 β, ∇ 2 P = 0, and ∇ 2 Q = 0.In our revised model, we obtain the same equations, except for the equation governing Q: ∇ 2 α = 0, ∇ 2 β = 1 σ 2 β, ∇ 2 P = 0, and The monodomain part of the model is the same in both cases (∇ 2 α = 0 and ∇ 2 P = 0) but the bidomain part is different.In the original model, the bidomain behavior was characterized by one length constant σ = νµ K(ν+µ) .In the revised model, the bidomain behavior is characterized by two length constants, σ and ξ = χδ K(χ+δ) .As discussed previously [4], the value of the membrane spring constant K is not known but is expected to be large, implying that σ is small.The size of ξ is not as clear.K is large, but so are χ and δ.In the truly incompressible limit, χ and δ go to infinity, so ξ becomes very large and the equation for Q becomes ∇ 2 Q = 0, as it was in the previous model.Note that in many of the expressions for the tissue displacement (e.g., (33)), a factor of ξ χ+δ is present, which goes to zero as χ and δ become very large, even if ξ is large.The length constant governing shearing, σ, is in general much smaller than the length constant for dilatation, ξ, because the bulk moduli are much greater than the shear moduli.The ratio ξ σ is about 300, independent of the value of K. We often speak of boundary layers of both pressure and displacement arising from the mechanical bidomain model.This view of a thin boundary layer is only useful in the limit when both σ and ξ are small compared to other length scales, such as the radius of the papillary muscle or blood vessel being modeled.
The solutions for the displacement in all three examples have a similar form: a monodomain part that depends on some power of the radius and which is the same in the intracellular and extracellular spaces, and a bidomain part that falls off with length constant σ or ξ.The monodomain part is typically larger than the bidomain part, but because it is the same in both spaces it contributes nothing to the membrane force, which is proportional to the difference between the intracellular and extracellular spaces.The bidomain part is different in the two spaces, so it alone contributes to the membrane force.In these three examples, the bidomain part contains a modified Bessel function, which solves the Helmholtz equation in two dimensions.Although these special functions are somewhat unfamiliar, at large values of their argument they behave qualitatively like exponentials (I as an increasing exponential, and K as a decaying exponential).So, if the bidomain length constants are small the bidomain part of the displacement falls off approximately exponentially with distance from the tissue boundary.Example 3 is particularly useful because that problem was solved completely using the previous model [5].The monodomain parts of the calculation are the same in both cases.The solution for the membrane force is similar but not identical.A term containing a modified Bessel function is present in both, and implies that both calculations result in a boundary layer at the edge of the tissue r = R of thickness σ.However, in the previous calculation [5] the membrane force also contained a term proportional to r, implying that the membrane force had a small contribution far from the boundary.This term is not present in our calculation.Thus, the effect of properly accounting for compressibility can make a difference in the model predictions, even in the incompressible limit.
The largest difference in Example 3 between the previous calculation and ours is the pressure.In the previous calculation [5], both the intracellular and extracellular pressures contained a term that varied as r 2 cos 2θ.The intracellular pressure also contained a constant term T  2 , but the extracellular pressure did not.Therefore, there were large differences between p and q throughout the tissue.In our calculation, p and q are both dominated by constant terms, and the only spatial dependence arises from terms containing a Bessel function.If ξ << R these Bessel function terms decay away from the surface, so p = T χ 2(χ+δ) and q = T δ 2(χ+δ) .However, if ξ >> R the Bessel function I 2 is approximately proportional to r 2 , so the pressure varies as r 2 cos 2θ like in the previous calculation.
The physical meaning of the intracellular and extracellular pressures has always been a confusing issue with the mechanical bidomain model.Reference [4] discussed the connection between macroscopic and microscopic properties of the model.According to that discussion, p and q are macroscopic pressures, and are related to the microscopic pressures p micro and q micro by p = Θ i p micro and q = Θ e q micro , where Θ i and Θ e are the intracellular and extracellular volume fractions.The shear moduli ν and µ are also macroscopic parameters.If we similarly take χ and δ to be macroscopic parameters, we find by analogy that χ = Θ i χ micro and δ = Θ e δ micro .If χ micro and δ micro are the same (both the bulk modulus of water) then χ δ = ν µ and the microscopic pressure difference p micro − q micro is proportional to Q and is zero except within a few length constants of the tissue edge.In other words, the bidomain pressure Q represents the difference between the microscopic intracellular and extracellular pressures, and any fluid flow between the two spaces would be driven by Q.
If one wants to use the previous model to do calculations, the primary change is the boundary conditions on the pressures.Instead of requiring that p = 0 at the boundary, one specifies that Q = 0, so p = χ χ+δ P and q = δ χ+δ P .The resulting calculations using the previous model will be correct, except in a boundary layer of thickness ξ.

V. CONCLUSION
In conclusion, consideration of tissue compressibility clarifies the behavior of the mechanical bidomain model.The monodomain contribution to the displacement and pressure are unchanged.The bidomain contribution to the displacement is similar to previous calculations, and it determines the membrane force.This paper shows that the bidomain contribution to the pressure is very different than was thought previously, and it determines the difference in microscopic pressure between the intracellular and extracellular spaces.The bidomain model predicts two new length constants: σ determines the width of the boundary layer for the membrane forces, and ξ determines the width of the boundary layer for the pressure.

Fig. 1 .
Fig.1.A schematic diagram of the mechanical bidomain model for a two-dimensional sheet of tissue.The elastic properties of the intracellular space are depicted by the lower grid of springs (green), and the properties of the extracellular space by the upper grid (blue).The two spaces are coupled by the membrane, shown as an array of springs (red).Because this is a two-dimensional model, stretching of the membrane springs is not caused by displacements in the z direction perpendicular to the sheet.Rather, if the intracellular and extracellular spaces are displaced by different amounts in the x-y plane, the membrane springs will stretch causing forces on the membrane.

Fig. 2 .
Fig. 2. A schematic diagram of Example 1: a cylinder of tissue of radius b immersed in a bath at pressure P0.

Fig. 3 .
Fig. 3.The pressures p and q as functions of r, for χ = 2δ and ξ = b 10 .The blue curve represents the intracellular pressure p and the red curve represents the extracellular pressure q.

Fig. 4 .
Fig. 4. A schematic diagram of Example 2: a cylindrical vessel of radius a containing fluid at pressure P0, surrounded by tissue.
Fig. 5. a) The monodomain part of the displacement, b) the bidomain part of the displacement ur − wr, and c) the pressures p and q, all as functions of r, for ν = 2µ , χ = 2δ, and ξ = a 10 .The displacements and pressures have been normalized to their values at r = a.

Fig. 6 .
Fig. 6.A schematic diagram of Example 3: a sheet of isotropic cardiac tissue of radius R with straight and uniform myocardial fibers oriented in the x-direction.

Fig. 7 .
Fig. 7. a) The monodomain part of the displacement, and b) the bidomain part of the displacement, as functions of position, for χ = 2δ, ν = 2µ, ξ = R 10 and σ = R 100 .The solid circle shows the tissue boundary with zero displacement, and the dashed oval shows how the tissue deforms when an active tension is present.The fibers are horizontal.The arrows in panels a) and b) are scaled differently; without this scaling the arrows in (b) would be much smaller.