PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Biomech Eng. Author manuscript; available in PMC 2010 October 1.
Published in final edited form as:
PMCID: PMC2860886
NIHMSID: NIHMS166600

CONTINUUM MODELING OF BIOLOGICAL TISSUE GROWTH BY CELL DIVISION, AND ALTERATION OF INTRACELLULAR OSMOLYTES AND EXTRACELLULAR FIXED CHARGE DENSITY

Abstract

A framework is formulated within the theory of mixtures for continuum modeling of biological tissue growth that explicitly addresses cell division, using a homogenized representation of cells and their extracellular matrix (ECM). The model relies on the description of the cell as containing a solution of water and osmolytes, and having a porous solid matrix. The division of a cell into two nearly identical daughter cells is modeled as the doubling of the cell solid matrix and osmolyte content, producing an increase in water uptake via osmotic effects. This framework is also generalized to account for the growth of ECM-bound molecular species that impart a fixed charge density (FCD) to the tissue, such as proteoglycans. This FCD similarly induces osmotic effects, resulting in extracellular water uptake and osmotic pressurization of the ECM interstitial fluid, with concomitant swelling of its solid matrix. Applications of this growth model are illustrated in several examples.

1. INTRODUCTION

Biological tissue growth can occur through a variety of mechanisms, including cell division and deposition of extracellular matrix (ECM) products synthesized by cells. Biological tissue growth may also be regulated by alterations in the intracellular concentration of osmolytes, which regulate cell volume. To date, continuum modeling of biological tissue growth appears to have focused mainly on the challenge of modeling deposition of solid matrix products [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], although some analyses have also addressed the growth of charged species, such as proteoglycans, that are fixed to the solid matrix [11, 12]. Modeling of chemical reactions among solutes is pervasive in chemistry and chemical engineering, and the incorporation of such reactions in the analysis of tissue growth has been demonstrated [12, 13, 14].

Other than in plant growth [15], less attention has been given to the porous nature of biological tissues and the role of osmotic effects in the control of water uptake and release. Since the growth of biological tissues involves significant uptake of intracellular and extracellular water, it is important to also consider osmotic mechanisms in the study of tissue growth.

In the current study, we formulate a framework for continuum modeling of biological tissue growth that explicitly addresses cell division, using a homogenized representation of cells and their ECM. The essential elements of this model rely on the description of the cell as containing a solution of water and osmolytes, and having osmotically inactive solid constituents that may be generically described as a porous solid matrix. The division of a cell into two nearly identical daughter cells normally starts with the duplication of cell contents during the synthesis phase, followed by cell division during the mitosis phase. Thus, ultimately, cell division is equivalent to doubling of the cell solid matrix and osmolyte content, and a resulting increase in water uptake via osmotic effects. In a homogenized representation of the tissue, the geometry of individual cells is not modeled explicitly, but their solid matrix and intracellular osmolyte content can be suitably incorporated into the analysis of the tissue response, thereby accounting for their osmotic effects. Thus, cell division can be described by the growth of these cell constituents, including the accumulation of osmotically active content, and the resultant uptake of water.

This framework is also generalized to account for the growth of ECM-bound molecular species that impart a fixed charge density (FCD) to the tissue, such as proteoglycans. This FCD similarly induces osmotic effects, resulting in extracellular water uptake and osmotic pressurization of the ECM interstitial fluid, with concomitant swelling of its solid matrix.

The formulation of this biological tissue growth framework relies on a homogenization procedure of the cell-ECM system, to overcome the challenge of modeling each cell explicitly. The essential elements of this homogenization rely on the existence of consistent frameworks for modeling the response of cells and porous-hydrated charged ECM to osmotic and mechanical loading.

The classical mathematical treatment for the passive response of isolated cells to osmotic loading with non-electrolytes was presented by Kedem and Katchalsky [16], who modeled the cell as a fluid-filled semi-permeable membrane. In a recent theoretical study using mixture theory, we extended their approach by modeling the protoplasm as a gel that may limit the intracellular solubility of a membrane-permeant osmolyte [17]; experimental support of this concept was subsequently provided in a biomimetic study of osmotic loading of alginate beads [18].

To model swelling of a cell within its surrounding matrix it is also necessary to take into account the swelling response of the matrix to osmotic changes in its external bathing environment, as shown recently by Haider et al. [19]. The triphasic theory of Lai et al. [20], and related theoretical frameworks [21, 22], may be used to analyze a charged matrix osmotically loaded with monovalent counter-ions. Similarly, the related mixture approach of Mauck et al. [23] may be used to analyze osmotic loading of a matrix by neutral solutes, accounting for their solubility in the porous matrix.

The current study combines many of these treatments and proposes a homogenization method for describing the swelling response of a tissue whose cellularity may vary from 0 to 100 percent of its volume. Furthermore, the framework is extended to account for biological tissue growth by cell division, via alteration of the intracellular concentration of osmolytes and solid matrix content, and by growth of ECM-bound charged molecular species.

2. BASIC DEFINITIONS

2.1. Apparent Density

In mixture theory [24, 25], multiple constituents α can co-exist within an elemental region, including solid and fluid constituents. The content of each constituent may be given by the apparent density ρα, where

ρα=dmαdV
(1)

Here, dmα represents the mass of constituent α in the elemental mixture volume dV. If the mixture includes solid constituents constrained to move together (collectively denoted by α = s), a natural choice is to define dV as a material region on this solid matrix; there is no flow of solid across the boundary surface of the elemental region dV. Therefore the only change that may occur in the elemental solid mass dms would be from chemical reactions (anabolic or catabolic reactions leading to resorption or growth, respectively). However, fluid constituents of the mixture (α ≠ s) may flow in or out of the region dV , on the understanding that the solid constituents of a solid-fluid mixture form a porous-permeable matrix.1

In the mixture theory framework employed here, each constituent α is considered to be intrinsically incompressible [26]. This means that the true density ρTα of the constituent is invariant, where

ρTα=dmαdVα
(2)

and dVα is the volume of constituent α in the elemental region dV. In particular, due to intrinsic incompressibility, dVs can only change as a result of chemical reactions. However, the elemental volume dV can change as a result of fluid flow in or out of the porous solid material region. Thus, even though the individual constituents are intrinsincally incompressible, the mixture is compressible as long as its pores remain open (also see Section 2.6 below).

2.2. Volume Fraction

There are no voids in a saturated mixture, thus

dV=αdVα
(3)

It is also possible to define the volume fraction of constituent α as

φα=dVαdV
(4)

such that

0φα1
(5)

The mixture saturation condition may also be written as

αφα=1
(6)

With this definition of [var phi]α, it follows that the apparent and true densities are related according to

ρα=φαρTα
(7)

implying that the apparent density is also bounded, 0ραρTα.

In biological mixtures, the fluid constituents generally consist of one solvent, typically water (α = w), and multiple solutes, generically denoted by α ≠ s, w. In the mixture framework employed here, it is assumed that the volume fraction of solutes is negligible, [var phi]α [double less-than sign] 1 (α ≠ s, w), consistent with the fluid mixture being a dilute solution, so that the mixture saturation condition may be approximated by

φs+φw1
(8)

The solid and solvent volume fractions are natural variables for describing solid and solvent content in a mixture of intrinsically incompressible constituents.

2.3. Concentration

For solutes, a natural variable describing their content is the concentration. Concentration may be defined on a mixture volume basis,

c˜α=dnαdV
(9)

or on a solution volume basis,

cα=dnα(1φs)dVdnαφwdVαs,w
(10)

where dnα is the number of moles of solute α in the elemental region dV. The latter definition is consistent with the common usage in chemistry; as shown in Section 2.5 below, the former definition is more suitable for modeling growth of solute content. Based on these definitions, concentration is related to apparent density via

ρα=Mαc˜α=Mα(1φs)cαMαφwcα
(11)

where the molecular weight of constituent α,

Mα=dmαdnα
(12)

is invariant.

2.4. Kinematics

In the current configuration the spatial position of an elemental mixture region is denoted by x. Since different mixture constituents may have different motions, each constituent α which is currently at x originated from possibly distinct reference locations Xα. The motion of each constituent is thus described by its own function

x=χα(Xα,t)
(13)

The solid s may itself consist of a mixture of solid constituents denoted generically by α = σ. In the current treatment, the solid constituents are constrained to move together at all times [7], thus they all share the same velocity vσ = vs. This constraint makes it meaningful to define the net solid apparent density as ρs = ∑σ ρσ and the net solid volume fraction as [var phi]s = ∑σ [var phi]σ.

Whether the reference configurations of the solid constituents are the same represents a critical consideration here. In the more general case, for distinct reference configurations Xσ, the motions χσ (Xσ, t) would be distinct, albeit constrained to all satisfy

χσ(Xσ,t)t=vs
(14)

In this case, the deformation gradient of each solid constituent would also be distinct and given by

Fσ=χσXσ
(15)

Clearly, the stress Teσ resulting from the strain in each solid constituent σ (the effective stress), each being a function of its own Fσ, would also be distinct and not reduce to zero in a common reference configuration. Furthermore, as the proportion [var phi]σ of each solid constituent evolves due to growth or resorption, the net effective stress in the solid, Tes=σTeσ, may similarly evolve due to evolving contributions from each constituent.2 This makes it impossible to define a time-invariant stress-free reference configuration Xs for the mixture of solid constituents, in a framework that accounts for growth. Therefore, in a general framework, it is necessary to keep track of multiple solid reference configurations, with motions satisfying Eq.(14) and distinct deformation gradients given by Eq.(15). We will address this general framework in greater detail in a future study, with the understanding that the specialized presentation in the remainder of this treatment represents an essential building block for this more general framework.

In the current study, we consider the special framework where all the solid constituents share the same reference configuration Xσ = Xs, such that all their motions are the same,

χσ(Xσ,t)=χs(Xs,t)
(16)

and all solid constituents share the same deformation gradient

Fs=χsXs
(17)

Thus, each of the stresses Teσ reduces to zero in the common reference configuration, and we need to keep track of only one deformation gradient, Fs.

2.5. Interstitial Growth

Interstitial growth is fundamentally described by the equation of conservation of mass,

ραt+div(ραvα)=ρ^α
(18)

where vα is the constituent’s velocity, and [rho with circumflex]α is the mass density supply to constituent α from chemical reactions. Supply terms are rate variables; in the notation of this study, they are denoted with a circumflex accent (^). The conservation of mass for the entire mixture requires that

αρ^α=0
(19)

In many studies of biological tissue growth, this constraint is not explicitly enforced because not all tissue constituents are modeled explicitly.

It is possible to dissociate changes in ρα resulting from growth and those resulting from solid deformation, by defining the apparent density, and the density supply, of α constituents relative to the stress-free reference configuration as

ρrα=Jsραρ^rα=Jsρ^α
(20)

where

Js=detFs=dVdVr
(21)

and dVr is the elemental mixture volume in its stress-free reference configuration, which remains invariant [10, 27, 28]. Thus,

ρrα=dmαdVr
(22)

In the case of solid constituents or molecular species fixed to the solid matrix, denoted by α = σ, dmσ may vary over time due to growth, thus ρrσ can describe this change in mass using an intensive variable, normalized to the invariant volume of the mixture in the reference configuration. Substituting Eq.(20) into Eq.(18) produces an evolution equation for ρrσ,

DsρrσDtρrσt+gradρrσ·vs=ρ^rσ
(23)

where we used vσ = vs. This relation is expressed in the spatial frame, where ρrσ=ρrσ(x,t). In a material frame, where ρrσ=ρrσ(Xs,t), this relation can be rewritten as

ρrσ(Xs,t)t=ρ^rσ(Xs,t)
(24)

and integrated to yield

ρrσ(Xs,t)=ρrσ(Xs,0)+0tρ^rσ(Xs,τ)dτ
(25)

It is evident from this relation that in the absence of growth of constituent σ(ρ^rσ=0), its apparent density relative to the stress-free reference configuration does not change over time.

In this relation, ρ^rσ must be provided via a constitutive relation that may depend on the density of any or all constituents present in the mixture, as well as the solid matrix strain, the ambient temperature, and the respective gradients of all these parameters; the rate of deformation of any constituent, its diffusion velocity relative to the mean mixture velocity; and the solid matrix texture tensors in the reference configuration [28].

For constituents α whose content is described more naturally by their volume fraction, we may define

φrα=dVαdVr
(26)

to represent the volume fraction of α relative to the reference configuration, where dV α is the elemental volume of constituent α in the current configuration; using the definitions of Eqs.(1), (2), and (4), it follows that

ρrα=φrαρTα
(27)

and

φrα=Jsφα
(28)

Since [var phi]α is bounded according to Eq.(5), φrα also satisfies

0φrαJs
(29)

The mixture saturation condition, Eq.(6), now yields

αφrα=Js
(30)

In the case of solid constituents α = σ, substituting Eq.(27) into Eq.(25) and recognizing the invariance of ρTσ yields the growth evolution equation for φrσ,

φrσ(Xs,t)=φrσ(Xs,0)+0tφ^rσ(Xs,τ)dτ
(31)

where φ^rσρ^rσ/ρTσ. Recognize that, unlike [var phi]σ, the upper bound on φrσ is not unity but Js, as given in Eq.(29). Thus, if Js > 1, φrσ can exceed unity. This is consistent with the expectation that a swollen tissue has more pore space available for interstitial solid growth.

The net solid volume fraction relative to the reference configuration is given by

φrs=σφrσ=Jsφs
(32)

Using the mixture saturation condition under the assumption of negligible solute volume fraction implies that

φrs+φrwJs
(33)

For molecules bound to the solid matrix, whose content is described more naturally by their concentration, we may define the concentration relative to the reference configuration as

c˜rσdnσdVr=Jsc˜σ
(34)

such that

ρrσ=Mσc˜rσ=Mσ(Jsφrs)cσMσφrwcσ
(35)

Thus, c˜rσ provides a measure of the change in the number of moles of constituent σ as a result of growth, in the form of an intensive variable normalized by the invariant mixture volume in the reference configuration. Substituting Eq.(35) into Eq.(25) yields the growth evolution equation for c˜rσ,

c˜rσ(Xs,t)=c˜rσ(Xs,0)+0tc^rσ(Xs,τ)dτ
(36)

where c^rσρ^rσ/Mσ. Since c˜rσ=(Jsφrs)cσ, this equation shows that, even in the absence of growth in constituent σ(c^rσ=0), changes can occur in the solution volume-based concentration cσ, either as a result of volume-altering deformation (changing Js) or the evolution of φrs from the growth of other solid constituents.

2.6. Pore Closure

Closure of the pores in a porous permeable solid matrix represents a limiting condition in a mixture of intrinsically incompressible constituents, whose implications are clarified here. Pore closure can occur as a result of deformation or growth.

Pore closure occurs when [var phi]s = 1, implying that Js=φrs according to Eq.(32). Once closure has occurred, there can be no further decrease in Js since there can be no further exudation of fluid, and the remaining solid is intrinsically incompressible. Thus, to accommodate pore closure, an analysis must recognize this transition between a compressible and incompressible mixture, with all the resulting implications (such as the enforcement of Jsφrs and the transition of the stress constitutive relation from compressible to incompressible formulations when Js=φrs). In transient analyses, pore closure can be avoided by formulating the constitutive relation for the permeability of the porous-permeable solid matrix such that it decreases asymptotically to zero as [var phi]s → 1. In this case, it will take an infinite time for the fluid to completely exude from the pores as they close, indefinitely postponing the need for transitioning from a compressible to an incompressible mixture.

An important implication of the current analysis is that interstitial growth of solid constituents occurs by progressively filling the pores of the tissue, thereby displacing the interstitial fluid. If the tissue swells as a result of growth, as may occur via osmotic mechanisms as shown below, the resulting pore swelling will compete with solid matrix deposition such that the pores do not necessarily close with increasing solid matrix deposition.

3. BIOLOGICAL TISSUE MODEL

The mixture framework can describe cells and the ECM, albeit using different mixture constituents. The total stress in a mixture is given by

T=pI+Tes
(37)

where p is the fluid pressure, I is the identity tensor, and Tes is the effective or extra stress, resulting from the deformation of the solid matrix. Under quasi-static conditions, and in the absence of external body forces, the conservation of linear momentum for the mixture requires that

divT=0
(38)

The current treatment is limited to the equilibrium response to osmotic swelling and mechanical loading, when all transient processes resulting from the flow of interstitial solvent and solutes have subsided. In this case, p represents only the osmotic pressure in the fluid, as the transient contribution from mechanical loading would have subsided. This simplification enables a clearer presentation of the salient physico-chemical concepts, and can be extended to transient processes subsequently. It also implies that the conservation of linear momentum for the fluid constituents (solvent and solutes), in the absence of external body forces, requires that their mechano-electrochemical potentials be uniform [20, 29].

The constitutive relations adopted for the mechano-chemical potential of solutes and solvent are those of ideal solutions [30], to avoid a burdensome notation; non-ideal behavior may be easily incorporated subsequently. Thus, for water, the mechano-(electro)chemical potential is given by

μ˜w=μ0w(θ)+1ρTw(pp0Rθαs,wcα)
(39)

where R is the universal gas constant, θ is the absolute ambient temperature (assumed uniform), and μ0w(θ) is the chemical potential in the reference physico-chemical state when p = p0 and cα = 0 [20, 22]. The summation term is taken over all solute constituents of the mixture.

Similarly, the mechano-electrochemical potential of the solutes is given by

μ˜α=μ0α(θ)+1Mα(zαFc(ψψ0)+Rθlncακαc0α)αs,w
(40)

where zα is the charge number and κα is the solubility of solute α in the mixture (0 < κα ≤ 1) [17, 20, 23]; Fc is Faraday’s constant and ψ is the electrical potential in the mixture; μ0α(θ) is the chemical potential of the solute in the reference physico-chemical state when ψ = ψ0 and cα=καc0α.3 For a tissue whose pore sizes vary over some distribution, κα may be thought of as the fraction of the pore volume that is accessible to the solute [23]. This constitutive relation is valid for dilute solutions where the volume fraction of solutes is negligible compared to the solvent volume fraction.

In general the solid matrix of the mixture may be electrically charged due to ionization of matrix-bound molecular species. The associated net FCD is denoted by cF (evaluated on a solution volume basis), and its charge sign by zF (zF = ±1).4 It is assumed that the mixture satisfies the electroneutrality condition,

zFcF+αs,wzαcα=0
(41)

The analysis of the tissue relies on a subdivision of the system into three compartments: the cells, the ECM, and the external bathing environment of the tissue. A parameter f is denoted by fc, fm and f* when associated with each of these respective compartments. The jump condition for the mechano-electrochemical potential of interface-permeant solutes and solvent across the interface between two compartments is given by

[[μ˜α]]=0αs
(42)

This jump condition,5 derived from the conservation of linear momentum and energy, and the entropy inequality across the interface, is valid in the absence of chemical reactions that produce or remove that particular solute or solvent on that interface, and in the absence of net interfacial free energy production or dissipation [28]. No jump condition applies for interface-impermeant species. In addition, a jump condition is provided by the conservation of linear momentum for the mixture,

[[pI+Tes]]·n=0
(43)

where n is the unit normal to the interface. This jump condition is valid in the absence of linear momentum supply on the interface; thus it neglects tension in the cell membrane [17].

The cell is treated as a gel filled semi-permeable membrane where the gel can limit the solubility of the solute in the protoplasm. In some cases, as indicated below, the stresses generated in the solid matrix of the protoplasmic gel may be neglected, Tce0, consistent with prior treatments [16, 17]. Also consistent with these prior treatments, the electrical charge of the protoplasmic molecules is not modeled explicitly here. Therefore, in the case of membrane-permeant solutes, the analysis presented here is limited to neutral osmolytes only; membrane-permeant neutral solutes are denoted by α = p, with zp = 0.

For membrane-impermeant osmolytes, their electrical charge does not influence the cell’s osmotic response, thus they are included in the analysis. In this treatment, the definition of “membrane-impermeant solutes” includes solutes that may slowly leak across the cell membrane via passive diffusion, but whose leakage is counter-balanced by active transport via membrane pumps (such as Na+, Cl, and K+ ions). If net transmembrane transport of these solutes occurs due to an imbalance between leaking and pumping, the resulting change in their concentration may be treated mathematically as a growth process. For simplicity, we limit the ECM solutes that are cell membrane-impermeant to monovalent counter-ions (such as Na+ and Cl), denoted by α = + and α = − , respectively, and a single neutral solute denoted by α = n (Table 2), with z+ = +1, z = −1, and zn = 0.6 Intracellular membrane-impermeant osmolytes, which play an important role in the regulation of the cell volume, are denoted by α = i.

TABLE 2
Notation for the constituents of a mixture.

It is further assumed that the monovalent counter-ions are so small relative to the pore space of the ECM that they are not excluded from any of the pores, κm+=κm=1. Note that in the external bath, devoid of a porous solid matrix, κ*α=1 for all solutes (α ≠ s, w); furthermore, since there can be no bath FCD, Eq.(41) requires that c*+=c*c*.

3.1. Single Cell and ECM

Consider a single cell embedded within the ECM. Since the ECM, cell membrane and protoplasm are all permeable to the water solvent, equilibrium conditions and the jump conditions on the water mechano-electrochemical potential, Eq.(42), yield

μ˜*w=μ˜mw=μ˜cw
(44)

Similarly, for the membrane-permeant neutral solute,

μ˜*p=μ˜mp=μ˜cp
(45)

For the membrane-impermeant ECM solutes, the conditions are limited to

μ˜*+=μ˜m+,μ˜*=μ˜m,μ˜*n=μ˜mn
(46)

Combining these relations with the constitutive equations (39)(40) and using the electroneutrality condition, Eq.(41), in the ECM yields

cmp=κmpc*p,ccp=κcpc*p
(47)

cmn=κmnc*n
(48)

cm+=zFcmF+(cmF)2+(2c*)22,  cm=zFcmF+(cmF)2+(2c*)22
(49)

ψm=ψ*+Rθ2FclnzFcmF+(cmF)2+(2c*)2zFcmF+(cmF)2+(2c*)2
(50)

pm=p*+Rθ[(cmF)2+(2c*)2+(κmn1)c*n+(κmp1)c*p2c*]
(51)

pc=p*+Rθ[(κcp1)c*p+cci(c*n+2c*)]
(52)

These expressions are arranged such that the right-hand-side depends on prescribed quantities in the external bath (c*,c*p,c*n,p*,ψ*),7 solubilities in the ECM and cell (κmp,κmn,κcp), the FCD in the ECM (cmF), and the intracellular concentration of membrane-impermeant solutes (cci).

Since the charges fixed to the solid matrix of the ECM must move with it, Eq.(36) constrains the variation in cmF according to

c˜mrF(t)=c˜mrF(0)+0tc^mrF(τ)dτ
(53)

where c˜mrF(t) is the FCD on a mixture volume basis in the current configuration, evaluated relative to the mixture volume in the stress-free reference configuration; the dependence of each term on Xs is made implicit in this expression, for notational simplicity. According to Eq.(35), the expression for cmF can be obtained from this relation via

cmF=c˜mrF/(Jsφmrs)
(54)

evaluated in the current configuration.

Similarly, since the intracellular membrane-impermeant solutes cannot leave the cell, they can also be considered solid matrix-bound so that

c˜cri(t)=c˜cri(0)+0tc^cri(τ)dτ
(55)

and

cci=c˜cri/(Jsφcrs)
(56)

Example 1

Consider the case of an isolated cell (no ECM), subjected to a change in external bath osmolarity (osmotic loading) using a membrane-impermeant neutral solute. In this case, c*p=c*=0 and Eq.(52) reduces to pc=p*+Rθ(ccic*n). Furthermore, when neglecting effective stresses in the protoplasm (Te,cs0), the jump condition of Eq.(43) reduces to pc = p*. Combining these two results implies that cci=c*n which, when substituted into Eq.(55), yields

Js=[1φcrs(0)]c*n(0)c*n(t)+φcrs(t)+1c*n(t)0tc^cri(τ)dτ
(57)

where c*n(0) represents the osmolarity of the external bath solution in the reference state (Js = 1, assumed to occur at t = 0). In the absence of intracellular solute and solid growth, c^cri=0andφ^crs=0, the resulting linear relationship between the relative cell volume Js and the ratio c*n(0)/c*n(t),

Js=[1φcrs(0)]c*n(0)c*n(t)+φcrs(0)
(58)

is known as the Boyle-van’t Hoff relation [31]. Most cells exhibit this linear response under osmotic loading and are called perfect osmometers [32]. Experimental measurements of Js for various values of c*n(0)/c*n(t) can be used to extract the value of the cell solid volume fraction φcrs(0); note that 1φcrs(0)=φcrw(0) is also known as the fraction of osmotically active water.

Alternatively, if the external bath environment remains unchanged, c*n(t)=c*n(0), but in the presence of growth of intracellular solute or solid, the cell volume ratio is given by

Js=1φcrs(0)+φcrs(t)+1c*n(0)0tc^cri(τ)dτ=1+0t[φ^crs(τ)+c^cri(τ)c*n(0)]dτ
(59)

where we used the relation of Eq.(31). This relation shows that an increase in intracellular solute concentration (c^cri>0), or solid content (φ^crs>0), increases the cell volume (Js > 1). This is a mathematical representation of the cell’s well-known volume regulation mechanism [33], illustrating how cells can regulate their volume by altering intracellular membrane-impermeant solute and solid content. This example also illustrates how growth of constituents produces volume changes via osmotic mechanisms; note in particular that Js increases linearly with increasing φcrs(t) such that there is no fixed upper bound on φcrs(t) according to Eq.(29).

3.2. Homogenization

The goal of homogenization is to model the cells collectively, together with the ECM, using a mixture representation. To avoid an excess of diacritics in our notation, the symbol f for a parameter may represent either the generic form of that parameter, or the homogenized value in the tissue; the meaning should be evident from the context.

Consider that the volume fraction of cells in an elemental region dV of the tissue is χ, so that 1 − χ is the volume fraction of the ECM (0 ≤ χ ≤ 1). The homogenized apparent density ρα of constituent α is given by

ρα=χρcα+(1χ)ρmα
(60)

Substituting Eq.(7) into this relation and recalling that ρTα is invariant produces an expression for the homogenized volume fraction of solvent in the cells and ECM,

φw=χφcw+(1χ)φmw
(61)

Similarly, substituting Eq.(11) into Eq.(60) and recalling that Mα is invariant, the homogenized concentration cα is found to be related to the concentrations in each compartment via

φwcα=χφcwccα+(1χ)φmwcmααs,w
(62)

For solutes that are not present in a particular compartment, the corresponding concentration is simply set to zero (ccα=0orcmα=0). The form of Eq.(62) can also be specialized to evaluate the homogenized FCD in the tissue, letting ccF=0 and yielding

φwcF=(1χ)φmwcmF
(63)

The homogenized solubility is similarly obtained,

φwκα=χφcwκcα+(1χ)φmwκmααs,w
(64)

In particular, substituting the solutions of Eqs.(47)(49) for concentrations in the single cell-ECM-bath system into these equations, it is found that the corresponding solutions for solute concentrations in the homogenized system reduce to

cp=κpc*pcn=κnc*nc+=12[zFcF+(cF)2+[(1χφcwφw)2c*]2]c=12[zFcF+(cF)2+[(1χφcwφw)2c*]2]
(65)

To find the homogenized interstitial fluid pressure p and electric potential ψ, note that the mechano-electrochemical potential is defined here in terms of energy per mass, so that the homogenized mechano-electrochemical potential of any fluid constituent is given by

ραμ˜α=χρcαμ˜cα+(1χ)ρmαμ˜mααs
(66)

Thus, using Eqs.(7) and (11) to relate ρw to [var phi]w, and ρα to cα, it is found that

φwμ˜w=χφcwμ˜cw+(1χ)φmwμ˜mw
(67)

φwcαμ˜α=χφcwccαμ˜cα+(1χ)φmwcmαμ˜mααs,w
(68)

Given that p and ψ appear as linear terms with constant coefficients in the constitutive relations (39)(40) for the mechano-electrochemical potential, the homogenized expressions for these quantities are given by

φwp=χφcwpc+(1χ)φmwpm
(69)

φwcαψ=χφcwccαψc+(1χ)φmwcmαψm
(70)

Using the cell-ECM-bath equilibrium solutions of Eqs.(50)(52), and recognizing that ccα=0 for α = +,− in the current treatment, the solutions for the homogenized system are thus given by

p=p*+Rθ[(cF)2+[(1χφcwφw)2c*]2+(κp1)c*p+χφcwφwcci+(1χφcwφw)κmnc*nc*n2c*]
(71)

ψ=ψ*+Rθ2FclnzFcF+(cF)2+[(1χφcwφw)2c*]2zFcF+(cF)2+[(1χφcwφw)2c*]2
(72)

A key aspect of this homogenization procedure is to recognize that solid deformations of the ECM and cells cannot remain distinct within an elemental region dV, as ECM and cells become mixture constituents coexisting at every point in the homogenized region. In this homogenized mixture, as seen above, the solid matrix constituents are described by φms,φcs, and the solid-bound molecular species by cmFandcci, where φmF1andφci1. Since all solid constituents in an elemental region are constrained to move together, they must share the same current configuration; as explained in Section (2.4), we also assume in the current treatment that they share the same reference configuration, so that the deformation gradient Fs, and relative change in the volume, Js = det Fs, also coincide for all solid constituents within dV. An important implication of this constraint is that the cell volume fraction χ does not vary with deformation in this homogenization scheme.

The homogenized tissue stress may be given by

T=χTc+(1χ)Tm=pI+χTe,cs+(1χ)Te,ms
(73)

where p is the homogenized interstitial fluid pressure, given in Eq.(71). Since the cell volume fraction χ also represents the area fraction on any arbitrary plane through the tissue, according to Delesse’s law [34], this homogenization scheme is consistent with the jump condition of Eq.(43). Since we opted in Section 2.4 to simplify the analysis by letting the cells and matrix share the same reference configuration, it follows that Te,cs=0andTe,ms=0 simultaneously in that configuration.

The last remaining requirement is to identify a suitable set of ambient conditions for the reference (stress-free) configuration. In a pure elasticity analysis, in the absence of residual stresses, the reference configuration is typically defined to be traction-free.8 Since the current analysis also includes osmotic pressure, which can produce residual stresses, we may additionally define the osmotic environment that can produce a stress-free state under traction-free conditions. In previous studies of charged tissues bathing in a monovalent salt environment [20], in the absence of explicit modeling of cells, an infinitely hypertonic state (c* → ∞) could be defined in which the fluid pressure p reduced to the ambient pressure p*, as can be observed from Eq.(71) under these suitable limiting conditions (χ=0,c*p=c*n=0). In that case the tissue would not be subjected to osmotic swelling, so that Tes could reduce to zero under traction-free conditions (Tn=pn+Tesn=p*n according to Eq.(43)). A hypertonic state can be reproduced experimentally by using a sufficiently large salt concentration c* to overwhelm the effect of the FCD in the tissue, making it plausible to achieve this reference configuration under experimental conditions.9

In the current analysis, where cells are explicitly modeled, this hypertonic state is mathematically plausible though it leads to a reference configuration for cells that is not physiologically desirable. This can be illustrated with the special case of a homogenized tissue with an uncharged ECM (cF = 0), subjected to the same bathing conditions (c*p=c*n=0). In this case, Eq.(71) predicts that

p=p*+Rθχφcwφw(cci2c*)
(74)

when using Eqs.(61)(62). Looking at this expression, it is not possible to expect pp* as c* → ∞, unless cci=2c*. In other words, the intracellular concentration of membrane-impermeant solute must be infinite in the reference configuration, which implies that the cell has been emptied of its water content. Not only is this condition physiologically undesirable, but it would also be too restrictive to assume that Te,msandTe,cs reduce to zero only in this limiting case. Furthermore, as explained in Section (2.6), emptying a cell of its water content implies that φcs=1andJs=φcrs; however, in a homogenized model of the tissue, the lower bound on J is given by Js|minmax(φcrs,φmrs), thus it may not even be possible to empty the cell of its water content if the ECM gets depleted first (φmrs>φcrs).

An alternative is to propose that the reference configuration corresponds to a charged tissue bathing in a NaCl solution of finite tonicity (c*p=c*n=0andc*=c*r), in the limiting case when its FCD is zero (cF = 0). The resulting osmotic pressure is also given by Eq.(74) and we similarly find that p = p* when cci=2c*r in the reference state, except that c*r is now finite.10 Thus the cell need not be devoid of its water under these conditions. In this reference state, both Te,msandTe,cs reduce to zero.

4. INTERSTITIAL GROWTH OF BIOLOGICAL TISSUE

In the current model, there are four solid or solid-bound constituents whose content may change as a result of growth: The solid matrix of the ECM, whose content is described by φmrs; the molecular species imparting the FCD to the ECM, whose content is described by c˜mrF; the solid matrix of the cells, described by φcrs; and the membrane-impermeant intracellular solute, described by c˜cri. In this homogenized representation of the tissue, the intracellular solute is considered solid-bound since it remains in the cell, and the cell is embedded in the solid matrix of the ECM. It is implicit that the molecular building blocks, cytokines, morphogens, and other factors required for growth, are available even though they are not modeled explicitly.

For the ECM and protoplasm solid matrix, the growth evolution relation of Eq.(31) can be used to track changes in φmrsandφcrs, given constitutive relations for φ^mrsandφ^crs. For the FCD and intracellular membrane-impermeant solute, Eq.(36) can similarly be used to track changes in c˜mrFandc˜cri, given constitutive relations for c^mrFandc^cri. (Recall that cmFandcci may change as a result of growth or volume deformation, according to Eqs.(54) and (56), but changes in c˜mrFandc˜cri can only occur due to growth.) Growth of any one of these four constituents may alter the osmotic environment of the tissue, leading to changes in water content and tissue volume, as illustrated in Example 1. In the sections below, we address specific biological processes in the context of changes in these constituents.

In some of the examples presented below, a finite element analysis is performed to illustrate specific phenomena. Finite element modeling for the framework presented in this study is conceptually similar to the approach required for finite elasticity [35], since the equation to be solved is that of Eq.(38), where the unknown is the solid matrix deformation, subject to the natural boundary condition of Eq.(43). The standard finite element approach of elasticity theory is extended to incorporate the osmotic pressure p into the stress tensor, as per Eq.(37), where p is given by Eq.(71). The dependence of p on Js produces a contribution to the elasticity tensor of the mixture, called the osmotic modulus, as described in our recent study [36]. The custom-written finite element code used in that study was extended here, to include the osmotic pressure of a homogenized tissue of cells and ECM, and to allow material properties to evolve arbitrarily over time in response to growth.

4.1. Extracellular Matrix Growth

Extracellular matrix growth can be described by an increase in φmrs, though the current framework restricts this mechanism to the assumption that the reference configuration of the newly deposited material is the same as that of the pre-existing matrix (Section 2.4). A common occurrence in tissue growth processes, where this assumption may not be so restrictive, is the case of ECM degradation, since the loss of solid matrix material points that all share a common reference configuration does not require the identification of new reference configurations.

The loss of solid constituent in the ECM is described by a decrease in φmrs(thusφ^mrs0). Since ECM material properties may depend on the solid content φmrs(orρmrs), typically in a nonlinear manner [37, 38], changes in the amount of matrix can be equivalently described by alterations in intrinsic material properties. Thus, if μs(ρmrs) represents a generic material property of the solid matrix, a constitutive relation validated from experiments may describe its dependence on ρmrs. (More generally, the material properties may be functions of the content of any of the solid or fluid constituents of the mixture, and the ambient temperature.) An example is used to illustrate this ECM degradation mechanism in the case of articular cartilage.

Example 2

The solid matrix of articular cartilage consists primarily of type II collagen and proteoglycans. During the early onset of osteoarthritis it has been observed that the collagen fibrillates and the cartilage matrix swells. This swelling is generally attributed to the Donnan osmotic pressure of the proteoglycans, whose content may not have decreased substantially in the early stage of the disease, encountering less resistance against swelling from the weakened collagen matrix.

To simulate this mechanism with the current framework, we assume for simplicity that the cell volume fraction, normally ~5% in adult cartilage, is negligible for the purposes of this analysis (χ = 0). Thus, the only properties that need to be specified are those of the ECM. Let φmrs=0.2andc˜mrF=120mEq/L based on a representative composition of cartilage [39], and assume that the external bathing environment of the tissue is isotonic saline (2c* = 300 mOsM, c*p=c*n=0) at body temperature (θ = 310 K). For simplicity, let the porous solid matrix be described by an isotropic neo-Hookean model, with a Poisson’s ratio of 0 and Young’s modulus initially set to 10 MPa, to represent the tensile properties of healthy cartilage.11 Under traction-free conditions, the tissue swells isotropically and homogeneously by 2.7% (Js = 1.027) relative to the reference configuration, under an osmotic pressure difference pp* = 86 kPa. This solution is obtained by solving the jump condition for the total traction, (pI+Tes)n=p*, where p is given by Eq.(71), for any arbitrary direction n.

It is now assumed that φmrs decreases to 0.18 as a result of collagen degradation and, as a further consequence of concomitant fibrillation, Young’s modulus decreases to 2 MPa, while c˜mrF remains the same. Then Js is observed to increases to 1.115, implying that the tissue swelled by 9% relative to its healthy initial state. The osmotic pressure difference decreases slightly as a result of this volume increase, pp* = 68 kPa; though c˜mrF was held constant, the FCD normalized to the interstitial fluid volume in the current configuration, cmF, reduced from 145 mEq/L to 128 mEq/L. As expected, this simple analysis shows that the initial swelling observed in the early stage of osteoarthritis can be explained by the influence of collagen degradation on the tensile properties of the cartilage ECM.

4.2. Synthesis of Solid-Bound Charged Species

An increase in the FCD, c˜mrF, may occur by cell synthesis of charged molecules which then bind to the ECM. All else remaining the same, the resulting increase in cmF will increase the osmotic pressure in the ECM, according to Eq.(71). This type of growth is illustrated with an example related to alterations in residual stresses in the aorta.

Example 3

The aortic wall is divided into the intima, media and adventitia. Unlike the collagenous adventitia, the intima and media have a significant proteoglycan content, synthesized by endothelial and vascular smooth muscle cells [40, 41]. In a recent experimental and theoretical study [36] we demonstrated that the associated FCD can play an important role in regulating transmural residual stresses, and the “opening angle” measured by cutting an arterial ring along a radial direction [42]. Our study demonstrated that decreasing the NaCl concentration in the external bath produced an increasingly larger opening angle, for a prescribed FCD in the intima and media.

Since proteoglycan content is known to increase during atherosclerosis [43, 44, 45], and since the opening angle has also been shown to increase with this advancing disease [46, 47, 48], it is reasonable to hypothesize that these two phenomena are directly related. In this example we simulate the growth of proteoglycan in the rat aorta, using some of the geometric and material parameters of our recent study [36]. The cell volume fraction in the intima and media is taken to be χ = 0.3, while that of the adventitia is assumed to be negligible (χ = 0). The intracellular and extracellular solid volume fractions relative to the reference state are constant and taken to be φmrs=φcrs=0.3 throughout the wall thickness. The intracellular concentration of membrane-impermeant solutes is set to c˜cri=210mOsM (a value that yields cci=c˜cri/(Jsφcrs)=300mOsM when Js = 1). The FCD is non-zero only in the intima and media; to model FCD growth, c˜mrF is increased from 40 mEq/L to 120 mEq/L over time. The external bath is physiological saline, with 2c* = 300 mOsM, c*p=0andc*n=0. Room temperature is assumed, θ = 298 K. The ECM solid matrix is modeled with a neo-Hookean constitutive relation, with a shear modulus of 140 kPa and Poisson’s ratio of 0 throughout the entire wall thickness. The finite element model represents an aortic ring cut along its radius; the objective is to analyze the evolution of the opening angle with growth of the FCD. Due to symmetry considerations, only one quarter of the ring is modeled.

Results of the finite element analysis are presented in Figure 1, showing that the opening angle increases significantly with increasing FCD, from 13° to 95° in this case. This effect can be attributed to the Donnan osmotic swelling induced by the FCD in the intima and media; since no swelling occurs in the adventitia, the intimal and medial swelling produce an inhomogeneous transmural stress that manifests itself as an opening angle upon sectioning of the aortic ring. This analysis provides theoretical support for the hypothesis that the increase in intimal and medial proteoglycan content (a growth process) accompanying atherosclerosis can drive the increase in opening angle.

FIGURE 1
Effect of growth of fixed charge density on opening angle of rat aorta. Due to symmetry, only one-quarter of a cut ring is displayed. The opening angle increases from 13° to 95° as c˜mrF increases from 40 mEq/L to 120 mEq/L in ...

4.3. Cell Volume Regulation

Intracellular concentration of membrane-impermeant solutes may change as part of the well-known mechanism for regulating cell volume [33]; this mechanism is primarily driven by alterations in c˜cri, either via active membrane transport processes, or the binding and release of intracellular solutes from substrates. There are many problems in growth and remodeling of biological tissues where such a process may play an important role. Two illustrative examples are provided below, relevant to conditions encountered in biological tissues.

Example 4

Increased intracranial pressure may be caused by osmotic swelling of cells or their ECM. An increase by 10 mmHg is clinically significant and potentially life threatening [49]. It has been shown that brain cells can increase their intracellular osmolarity via active uptake of sodium, potassium and chloride ions, or production of organic osmolytes generically described as “idiogenic osmoles” [33].

In the current framework, an elementary analysis can be performed, using a simplified representation of the brain, that can predict this increase in intracranial pressure. Consider that the volume of the brain does not change in situ (Js =constant) due to the physical constraints of the cranium. The blood supply represents the “external bath” in this simplification, with pressure p* and osmolarity 2c*, and c*p=c*n=0. Though brain tissue contains proteoglycans, it is assumed for simplicity that cmrF=0 in this example, so that we can let Js = 1. Then, Eq.(71) simplifies to

pp*=Rθχφcwφw(cci2c*)=Rθχφcwφw(c˜cri1φcrs2c*)
(75)

Clearly, all else being equal, an increase in c˜cri will lead to an increase in the intracranial pressure relative to the blood pressure, pp*. For example, if we let χ = 0.85, φcrs=φmrs=0.3 (so that φcw/φw=1), θ = 310 K, 2c* = 300 mOsM, and let the initial value c˜cri yield zero pressure difference (c˜cri(0)=(1φcrs)2c*=210mOsM), then an increase in pressure by 10 mmHg (pp* = 1.3 kPa) can be achieved if c˜cri increases by as little as 0.42 mOsM.

This analysis is simplified of course, since it does not take into account the compliance introduced by vasoconstriction, dura mater deformation, and drainage of cerebrospinal fluid. Thus, in reality, it would take larger increases in c˜cri to produce this change in intracranial pressure. However, this analysis conveys the salient aspect of the role of cell volume regulatory mechanisms on intracranial pressure.

Example 5

In Example 3, it was shown that the opening angle of the aorta increases with increasing FCD. In that analysis, in the limit of zero FCD, it can also be shown that the opening angle reduces to zero, since all residual stresses subside in this limiting condition. However, experimental measurements have shown that the opening angle may vary significantly along the length of the aorta, and in some cases it is observed to be negative. This outcome, which cannot be predicted from Donnan swelling, suggests that residual stresses in the aortic wall cannot be explained exclusively by the osmotic swelling resulting from the FCD, as acknowledged in our earlier study [36].

In the cardiovascular biomechanics literature, the aortic wall residual stresses manifested by the opening angle have been attributed to other factors, including “incompatible growth” of the elastin and collagen constituents [50], and osmotic swelling of the endothelial and vascular smooth muscle cells of the intima and media [51]. Both of these hypotheses are plausible explanations for these residual stresses; in the current tissue modeling framework we are able to address the latter mechanism (Section 2.4).

In this example we test the hypothesis that the negative opening angle observed in the aorta could result from a decreasing concentration of the intracellular membrane-impermeant solutes, via the normal mechanisms of cell volume regulation. A model of the rat aorta is created, similar to that used in Example 3, except that c˜mrF is constant and set to 40 mEq/L in the intima and media. In this analysis, the parameter allowed to change as a result of growth is the intracellular concentration c˜cri of membrane-impermeant solute in the intima and media, which is decreased from 210 mOsM to 175 mOsM.

Results show that the initial opening angle is positive and equal to 13° when c˜cri=210mOsM. However, as c˜cri decreases, the cell volume also decreases as a result of water exudation, and the opening angle becomes negative, reducing to −37° when c˜cri=175mOsM (Figure 2). This outcome supports the hypothesis that the opening angle of arteries may be regulated by the intracellular concentration of membrane-impermeant solutes, complementing other mechanisms such as Donnan swelling due to the FCD.

FIGURE 2
Effect of growth of intracellular membrane-impermeant solute concentration on opening angle of rat aorta. Due to symmetry, only one-quarter of a cut ring is displayed. The opening angle decreases from 13° to −37° as c˜ ...

Though this example illustrates the reduction in opening angle with decreasing c˜cri, it can also be shown that the opening angle can increase with increasing c˜cri. This potential mechanism is of particular interest in light of the observation that the opening angle of the aorta can change over a period as little as two days following a ligation procedure [46], a time duration consistent with the relatively rapid process of cell volume regulation.

4.4. Cell Division

The net effect of cell division can be modeled in the current framework by allowing three parameters to change simultaneously, namely, the cell solid volume fraction φcrs, the intracellular membrane-impermeant solute concentration, c˜cri, and the cell volume fraction, χ. The evolution in φcrsandc˜cri occurs because the division of a cell into two daughter cells follows an initial synthesis phase during which the parent cell doubles its content of osmotically inactive and osmotically active molecules. Consequently, in the case of cell division, it is expected that the growth rates φ^crsandc^cri should be proportional to each other.

Everything else being equal, the evolution in χ occurs because the total volume of cells increases relative to that of the extracellular matrix. The cell volume fraction may be expressed as

χ=dVcdVm+dVcξ1+ξ
(76)

where dVm is the elemental mixture volume of the ECM and dVc is the elemental mixture volume of cells in the elemental homogenized mixture volume dV, and ξ = dVc/dVm; since dV = JsdVr, dVm = JsdVmr and dVc = JsdVcr in this homogenization framework, the expressions for χ and ξ can be written in the same form when using the elemental volumes in the reference configuration. Two limiting cases can be noted here. In the limit when the matrix volume is negligible (confluent cells with negligible ECM), ξ tends to infinity and χ approaches unity; in this case, though cell division further increases ξ, it does not cause measurable changes in χ. In the opposite limit, when there are no cells, ξ = 0 and χ = 0.

The rate at which cells divide may be represented by the function [Xi w/ hat] = [partial differential]ξ (X, t)/[partial differential]t, such that

ξ(X,t)=ξ(X,0)+0tξ^(X,τ)dτ
(77)

On the assumption that the volume of daughter cells is nearly identical to the parent cell, it is also reasonable to assume that [Xi w/ hat] is proportional to φ^crsandc^cri,

c^criφ^crsξ^
(78)

This constraint is just a guideline that facilitates the analysis of cell division in this homogenization framework. Deviations from this rule can occur and may be addressed in the appropriate context.

Example 6

In long bone morphogenesis, the early stage of bone formation starts with the development and growth of a hyaline cartilage model, followed by the development of primary and secondary ossification centers, then the formation of the articular layers and epiphyseal plate. In this example, the growth of the hyaline cartilage model by chondrocyte division is described using a cell division analysis. For simplicity it is assumed that χ = 1; thus, the only parameters that need to be specified explicitly are those related to the cells. The external bath is assumed to be normal saline.

The intracellular solid content φcrs and concentration of membrane-impermeant solute c˜cri are both assumed to increase by an arbitrary factor of six at a constant rate, φcrs=0.2to1.2,andc˜cri=240mOsM to1440mOsM. From symmetry considerations, only one octant of the geometry is modeled. Results demonstrate that the cartilage model grows homogeneously and isotropically (Figure 3), increasing in volume by a factor of six (from Js = 1 to 6). This outcome represents the simplest analysis of growth by cell division, where the ECM (non-existent in this example) has no influence on the outcome. Since there are no extra stresses in the ECM (Te,ms=0) or cell (Te,cs=0), the osmotic pressure in the tissue is equal to the ambient bath pressure, pp* = 0. Importantly, the intracellular solute concentration in the current configuration, cci=c˜cri/(Jsφcrs), remains constant and equal to 300 mOsM during this growth process. This example also illustrates a problem where φcrs exceeds unity without violating the constraint of Eq.(29).

FIGURE 3
Growth of hyaline cartilage model in long bone morphogenesis. Due to symmetry, only one octant of the model is shown. In this analysis growth occurs by cell division. Sixfold increases in φcrsandc˜cri lead to a sixfold increase in tissue ...

Example 7

In this example we reprise the analysis of growth of a hyaline cartilage model, this time including an ECM. Let φmrs=0.2, and let the ECM be described by a neo-Hookean model for Te,ms, with shear modulus equal to 0.5 MPa and Poisson’s ratio equal to 0. The intracellular contents grow at a constant rate over time as described in Example 6. The cell volume fraction is initially χ = 0.8, equivalent to ξ = 4. ξ is assumed to increase by a factor of six at a constant rate, to ξ = 24, yielding a final cell volume fraction of χ = 0.96. The increase in χ over time has a non-constant rate dictated by Eq.(76).

Results show that the growth is substantially similar to the ECM-free case described in Figure 3, with a final value of Js = 5.95, homogeneous over the entire hyaline cartilage model. The small difference in Js, as compared to the ECM-free case, occurs because the tensile stress in the ECM offers some resistance to the osmotic swelling mechanism accompanying cell division. This tensile stress is balanced by a small non-zero osmotic swelling pressure difference, pp* = 8 kPa in this example. This outcome suggests that osmotic swelling easily overcomes the tensile stresses in the ECM, for the material parameters used here.

5. DISCUSSION

The objective of this study was to formulate a mathematical framework for modeling growth of biological tissues, with a particular emphasis on cell division and the growth of constituents that alter the osmotic environment of the interstitial fluid. This framework was developed in the context of the theory of mixtures, which can accommodate any number of solid and fluid constituents. Though the general framework of growth is applicable to any materials undergoing chemical reactions, the specific application to biological tissues stems primarily from the modeling of cells within their (possibly charged) ECM, with a mixture whose constituents are outlined in Table 2.

Incorporating cells within their ECM was achieved by focusing on mechano-electrochemical balances between intracellular, extracellular, and external bathing environments. A tissue homogenization procedure was introduced to combine the contributions of intracellular and extracellular compartments to all the dependent variables in the analysis. The influence of cells on the elasticity of the tissue’s solid matrix was not examined explicitly in this homogenization framework, where a simple rule of mixtures was adopted for Tes and the dependence of solid matrix elasticity on tissue content, including cell volume fraction, was relegated to suitable constitutive modeling. For example, other studies have examined the topic of modeling cells as inclusions in an infinite medium [52], and such frameworks could be adopted here.

A basic principle of this study is that the growth of solid constituents is described by changes in their apparent density, based on the equation of conservation of mass, Eq.(18). The effects of growth and deformation on ρσ are dissociated by introducing the apparent density relative to the reference configuration, ρrσ, in Eq.(20). Thus, changes in ρrσ only result from growth, according to the growth evolution equation (25).

Changes in ρrσ do not necessarily imply changes in the volume of the tissue. The examples presented in Section 4 illustrate various possibilities. In the case of cartilage collagen degradation, tissue volume increased due to swelling caused by proteoglycans, which encountered less resistance against the weakened collagen matrix; had the cartilage been devoid of proteoglycans, the collagen degradation in this example would not have caused changes in tissue volume, but only an increase in tissue porosity. In the example of cell volume regulation in the brain, the boundary conditions were set to prevent changes in brain volume, thus growth only led to changes in interstitial fluid pressure; however, if the brain had been allowed to swell, as occurs in emergency treatments when decompressive craniectomy is performed thereby removing a portion of the skull to relieve intracranial pressure [53], the analysis could predict this volume increase by releasing the displacement constraints on the brain surface. In the example of long bone growth, an increase in volume (mostly by water uptake) was indeed predicted in direct proportion to changes in ρrσ for the intracellular solid matrix content and membrane-impermeant solutes; however, as further illustrated, this swelling growth could be constrained by increasing the stiffness of the ECM (as also shown in the cartilage degradation example). Thus, interstitial growth does not necessarily imply volume changes; it can influence tissue volume, or interstitial osmotic fluid pressure, or both.

The growth mechanisms described in this study may result in the evolution of residual stresses in tissues in a predictable manner. When the interstitial osmotic pressure increases inhomogeneously as a result of selective growth of constituents, the stress distribution in the tissue solid matrix may also become inhomogeneous under a traction-free configuration, leading to residual stresses, as illustrated in the case of proteoglycan growth or cell volume regulation in the intima and media of the aorta. Similarly, inhomogeneous cell division may also produce such residual stresses, though this case was not specifically illustrated here.

The framework presented here was specialized to the case where all solid matrix constituents share a common reference configuration. The topic of differential or incompatible growth of solid constituents and evolving reference configurations remains an active area of investigation, but this study did not place a focus on it, addressing instead the role of osmotic effects on growth. Nevertheless, the presentation of Section 2.4 provides the foundation of an approach that we believe can address the topic of growth of multiple generations of solid constituents, each having its own reference configuration, though all constrained to share the same current configuration, see Eqs.(14)(15), a foundation shared with the approach proposed by Humphrey and Rajagopal [7]. In this context, a “generation” consists of any number of solid constituents that share the same reference configuration, though their content may vary over time due to growth or degradation; thus, the specialization adopted in our subsequent treatment may be viewed as that applicable to a single generation of solid matrix constituents.

When a framework becomes available for describing multiple generations of solid matrix constituents, it will become possible to address such problems as the initial growth of a tissue by cell division, followed by selective cell apoptosis and their replacement with newly synthesized extracellular matrix. Other applications may include the remodeling of actin cytoskeletal structure in response to osmotically driven volume changes in cells. The other mode of tissue growth not addressed in this study is that of appositional, or surface growth [3, 54]. A framework for describing surface growth in the context of mixture theory has been proposed in our earlier study [28].

In summary, the current framework emphasizes the role of osmotic effects on the uptake of water into the tissue. A significant benefit of giving due consideration to these effects is the ability to model cell division as shown above. Since all biological soft tissues are hydrated, the premise of our approach is that the inclusion of osmotic effects and water uptake is essential to capture salient growth mechanisms in such tissues. This study provides a formulation for this framework and several illustrations that demonstrate its application.

TABLE 1
Nomenclature.

ACKNOWLEDGMENTS

This study was supported with funds from the National Institute of Arthritis and Musculoskeletal and Skin Diseases of the U.S. National Institutes of Health (AR46532, AR52871).

Footnotes

1In the special case when a fluid cannot flow, due to being entirely trapped within pores or otherwise bound to the solid, it may be modeled as part of the solid matrix.

2In mixture theory, the effective stress Teσ in a constituent is defined such that its associated traction vector on a surface represents the force acting on that constituent, normalized by the elemental area of the mixture. It is thus an apparent stress; apparent stresses can be added together to provide a net apparent stress, and its associated traction vector on an elemental mixture area.

3The values of p0, ψ0, c0 and the corresponding reference chemical potentials μ0α(θ) at a given θ are prescribed standard states that remain invariant [30].

4If the FCD results from a homogeneous molecular species σ of concentration cσ and charge number zσ, then zσcσ = zF cF.

5Let the interface divide the domain into two regions, denoted by ’+’ and ’−’, then [[f]] n [equivalent] f+n+ + fn = (f+f) n, where n is the unit outward normal to the ’+’ side.

6Using monovalent counter-ions yields a closed-form solution for the Donnan osmotic pressure in the ECM, as shown below; multivalent ions can be analyzed similarly but require numerical solutions. Analyzing multiple neutral solutes (whether membrane-permeant or impermeant) requires a straightforward superposition of the analysis of a single solute.

7It is common to let p* = 0 such that pm and pc represent gauge pressures; it is similarly convenient to let ψ*= 0.

8In this treatment where we explicitly model the ambient bath pressure p*, traction-free should be understood to mean that the applied traction is only that of the ambient pressure.

9However, there may be other physico-chemical factors, such as osmotic pressure resulting from configurational entropy of the matrix-bound charged molecular species, that may not vanish under hypertonic conditions. Thus, as often encountered in continuum mechanics, a true stress-free reference configuration may not be achievable experimentally, but only mathematically.

10For example, c*r = 0.15 M to reproduce common physiological conditions, though any other value may be chosen.

11Recall that the porous solid matrix is compressible because its pores may change their volume during deformation, even though the matrix skeleton and the interstitial fluid are modeled as intrinsically incompressible. Thus, it is acceptable to have a Poisson’s ratio equal to zero in this framework.

REFERENCES

1. Hsu FH. The influences of mechanical loads on the form of a growing elastic body. J Biomech. 1968;1(4):303–311. [PubMed]
2. Cowin SC, Hegedus DH. Bone remodeling - 1. theory of adaptive elasticity. J Elasticity. 1976;6(3):313–326.
3. Skalak R, Dasgupta G, Moss M, Otten E, Dullumeijer P, Vilmann H. Analytical description of growth. J Theor Biol. 1982;94(3):555–577. [PubMed]
4. Rodriguez EK, Hoger A, McCulloch AD. Stress-dependent finite growth in soft elastic tissues. J Biomech. 1994;27(4):455–467. [PubMed]
5. Epstein M, Maugin GA. Thermomechanics of volumetric growth in uniform bodies. Int J Plasticity. 2000;16(7–8):951–978.
6. Klisch S, Van Dyke T, Hoger A. A theory of volumetric growth for compressible elastic biological materials. Math Mech Solids (USA) 2001;6(6):551–575.
7. Humphrey JD, Rajagopal KR. A constrained mixture model for growth and remodeling of soft tissues. Math Mod Meth Appl S. 2002;12(3):407–430.
8. Garikipati K, Arruda E, Grosh K, Narayanan H, Calve S. A continuum treatment of growth in biological tissue: The coupling of mass transport and mechanics. J Mech Phys Solids. 2004;52(7):1595–1625.
9. Volokh KY. Mathematical framework for modeling tissue growth. Biorheology. 2004;41(3–4):263–269. [PubMed]
10. Guillou A, Ogden RW. Growth in soft biological tissue and residual stress development. In: Holzapfel GA, Ogden RW, editors. Mechanics of Biological Tissue. Springer; 2006. pp. 47–62.
11. Klisch SM, Chen SS, Sah RL, Hoger A. A growth mixture theory for cartilage with application to growth-related experiments on cartilage explants. J Biomech Eng. 2003;125(2):169–179. [PubMed]
12. Dimicco MA, Sah RL. Dependence of cartilage matrix composition on biosynthesis, diffusion, and reaction. Transport Porous Med. 2003;50(1–2):57–73.
13. Radisic M, Deen W, Langer R, Vunjak-Novakovic G. Mathematical model of oxygen distribution in engineered cardiac tissue with parallel channel array perfused with culture medium containing oxygen carriers. Am J Physiol Heart Circ Physiol. 2005;288(3):H1278–H1289. [PubMed]
14. Lemon G, King JR, Byrne HM, Jensen OE, Shakesheff KM. Mathematical modelling of engineered tissue growth using a multiphase porous flow mixture theory. J Math Biol. 2006;52(5):571–594. [PubMed]
15. Cosgrove D. Biophysical control of plant cell growth. Annu Rev Plant Physiol. 1986;37:377–405. [PubMed]
16. Kedem O, Katchalsky A. Thermodynamic analysis of the permeability of biological membranes to non-electrolytes. Biochim Biophys Acta. 1958;27(2):229–246. [PubMed]
17. Ateshian GA, Likhitpanichkul M, Hung CT. A mixture theory analysis for passive transport in osmotic loading of cells. J Biomech. 2006;39(3):464–475. [PMC free article] [PubMed]
18. Albro MB, Chahine NO, Caligaris M, Wei VI, Likhitpanichkul M, Ng KW, Hung CT, Ateshian GA. Osmotic loading of spherical gels: a biomimetic study of hindered transport in the cell protoplasm. J Biomech Eng. 2007;129(4):503–510. [PMC free article] [PubMed]
19. Haider MA, Schugart RC, Setton LA, Guilak F. A mechano-chemical model for the passive swelling response of an isolated chondron under osmotic loading. Biomech Model Mechanobiol. 2006;5(2–3):160–171. [PubMed]
20. 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–258. [PubMed]
21. Huyghe JM, Janssen JD. Quadriphasic mechanics of swelling incompressible porous media. Int J Eng Sci. 1997;35(8):793–802.
22. Gu WY, Lai WM, Mow VC. A mixture theory for charged-hydrated soft tissues containing multi-electrolytes: passive transport and swelling behaviors. J Biomech Eng. 1998;120(2):169–180. [PubMed]
23. Mauck RL, Hung CT, Ateshian GA. Modeling of neutral solute transport in a dynamically loaded porous permeable gel: implications for articular cartilage biosynthesis and tissue engineering. J Biomech Eng. 2003;125(5):602–614. [PMC free article] [PubMed]
24. Truesdell C, Toupin R. The classical field theories. In: Flugge S, editor. Handbuch der Physik. Vol. III/1. Berlin: Springer-Verlag; 1960.
25. Bowen RM. Theory of mixtures. In: Eringen AE, editor. Continuum physics. Vol. 3. New York: Academic Press; 1976. pp. 1–127.
26. Bowen RM. Incompressible porous media models by use of the theory of mixtures. Int J Eng Sci. 1980;18(9):1129–1148.
27. Bowen RM. The thermochemistry of a reacting mixture of elastic materials with diffusion. Arch Ration Mech An. 1969;34(2):97–127.
28. Ateshian GA. On the theory of reactive mixtures for modeling biological growth. Biomech Model Mechanobiol. 2007;6(6):423–445. [PMC free article] [PubMed]
29. Katzir-Katchalsky A, Curran PF. Nonequilibrium thermodynamics in biophysics. Harvard books in biophysics ; no. 1. Cambridge: Harvard University Press; 1965.
30. Tinoco I, Jr, Sauer K, Wang JC. Physical chemistry : principles and applications in biological sciences. Prentice Hall: 1995.
31. Nobel PS. The boyle-van’t hoff relation. J Theor Biol. 1969;23(3):375–379. [PubMed]
32. Weiss TF. Cellular biophysics. Cambridge, Mass: MIT Press; 1996.
33. McManus ML, Churchwell KB, Strange K. Regulation of cell volume in health and disease. N Engl J Med. 1995;333(19):1260–1266. [PubMed]
34. Delesse A. Pour déterminer la composition des roches. Annales de mines. 1848;4(13):379.
35. Bonet J, Wood RD. Nonlinear continuum mechanics for finite element analysis. Cambridge ; New York, NY: Cambridge University Press; 1997.
36. Azeloglu EU, Albro MB, Thimmappa VA, Ateshian GA, Costa KD. Heterogeneous transmural proteoglycan distribution provides a mechanism for regulating residual stresses in the aorta. Am J Physiol Heart Circ Physiol. 2008;294(3):H1197–H1205. [PubMed]
37. Carter DR, Hayes WC. Bone compressive strength: the influence of density and strain rate. Science. 1976;194(4270):1174–1176. [PubMed]
38. Gibson LJ, Ashby MF. Cellular solids : structure and properties. 2nd ed. Cambridge ; New York: Cambridge solid state science series. Cambridge University Press; 1997.
39. Mow VC, Huiskes R. Basic orthopaedic biomechanics & mechano-biology. 3rd ed. Philadelphia, Pa. ; London: Lippincott Williams & Wilkins; 2005.
40. Wight TN, Ross R. Proteoglycans in primate arteries. i. ultrastructural localization and distribution in the intima. J Cell Biol. 1975;67(3):660–674. [PMC free article] [PubMed]
41. Yao LY, Moody C, Schonherr E, Wight TN, Sandell LJ. Identification of the proteoglycan versican in aorta and smooth muscle cells by dna sequence analysis, in situ hybridization and immunohistochemistry. Matrix Biol. 1994;14(3):213–225. [PubMed]
42. Chuong CJ, Fung YC. On residual stresses in arteries. J Biomech Eng. 1986;108(2):189–192. [PubMed]
43. Porterfield SP, Calhoon TB, Weiss HS. Changes in connective tissue colloidal charge density with atherosclerosis and age. Am J Physiol. 1968;215(2):324–329. [PubMed]
44. Volker W, Schmidt A, Oortmann W, Broszey T, Faber V, Buddecke E. Mapping of proteoglycans in atherosclerotic lesions. Eur Heart J. 1990;11 Suppl E:29–40. [PubMed]
45. Evanko SP, Raines EW, Ross R, Gold LI, Wight TN. Proteoglycan distribution in lesions of atherosclerosis depends on lesion severity, structural characteristics, and the proximity of platelet-derived growth factor and transforming growth factor-beta. Am J Pathol. 1998;152(2):533–546. [PubMed]
46. Matsumoto T, Hayashi K, Ide K. Residual strain and local strain distributions in the rabbit atherosclerotic aorta. J Biomech. 1995;28(10):1207–1217. [PubMed]
47. Valenta J, Svoboda J, Valerianova D, Vitek K. Residual strain in human atherosclerotic coronary arteries and age related geometrical changes. Biomed Mater Eng. 1999;9(5–6):311–317. [PubMed]
48. Gregersen H, Zhao J, Lu X, Zhou J, Falk E. Remodelling of the zero-stress state and residual strains in apoe-deficient mouse aorta. Biorheology. 2007;44(2):75–89. [PubMed]
49. Treggiari MM, Schutz N, Yanez ND, Romand JA. Role of intracranial pressure values and patterns in predicting outcome in traumatic brain injury: a systematic review. Neurocrit Care. 2007;6(2):104–112. [PubMed]
50. Humphrey JD, Rajagopal KR. A constrained mixture model for arterial adaptations to a sustained step change in blood flow. Biomech Model Mechanobiol. 2003;2(2):109–126. [PubMed]
51. Guo X, Lanir Y, Kassab GS. Effect of osmolarity on the zero-stress state and mechanical properties of aorta. Am J Physiol Heart Circ Physiol. 2007;293(4):H2328–H2334. [PubMed]
52. Guilak F, Mow VC. The mechanical environment of the chondrocyte: a biphasic finite element model of cell-matrix interactions in articular cartilage. J Biomech. 2000;33(12):1663–1673. [PubMed]
53. Rangel-Castilla L, Gopinath S, Robertson CS. Management of intracranial hypertension. Neurol Clin. 2008;26(2):521–541. x. [PMC free article] [PubMed]
54. Skalak R, Farrow DA, Hoger A. Kinematics of surface growth. J Math Biol. 1997;35(8):869–907. [PubMed]