|Home | About | Journals | Submit | Contact Us | Français|
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.
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 , 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 , 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 ; experimental support of this concept was subsequently provided in a biomimetic study of osmotic loading of alginate beads .
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. . The triphasic theory of Lai et al. , 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.  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.
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
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 . This means that the true density of the constituent is invariant, where
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).
There are no voids in a saturated mixture, thus
It is also possible to define the volume fraction of constituent α as
The mixture saturation condition may also be written as
With this definition of α, it follows that the apparent and true densities are related according to
implying that the apparent density is also bounded, .
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, α 1 (α ≠ s, w), consistent with the fluid mixture being a dilute solution, so that the mixture saturation condition may be approximated by
The solid and solvent volume fractions are natural variables for describing solid and solvent content in a mixture of intrinsically incompressible constituents.
For solutes, a natural variable describing their content is the concentration. Concentration may be defined on a mixture volume basis,
or on a solution volume basis,
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
where the molecular weight of constituent α,
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
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 , 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 s = ∑σ σ.
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
In this case, the deformation gradient of each solid constituent would also be distinct and given by
Clearly, the stress 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 σ of each solid constituent evolves due to growth or resorption, the net effective stress in the solid, , 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,
and all solid constituents share the same deformation gradient
Thus, each of the stresses reduces to zero in the common reference configuration, and we need to keep track of only one deformation gradient, Fs.
Interstitial growth is fundamentally described by the equation of conservation of mass,
where vα is the constituent’s velocity, and α 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
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
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 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 ,
where we used vσ = vs. This relation is expressed in the spatial frame, where . In a material frame, where , this relation can be rewritten as
and integrated to yield
It is evident from this relation that in the absence of growth of constituent , its apparent density relative to the stress-free reference configuration does not change over time.
In this relation, 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 .
For constituents α whose content is described more naturally by their volume fraction, we may define
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
Since α is bounded according to Eq.(5), also satisfies
The mixture saturation condition, Eq.(6), now yields
where . Recognize that, unlike σ, the upper bound on is not unity but Js, as given in Eq.(29). Thus, if Js > 1, 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
Using the mixture saturation condition under the assumption of negligible solute volume fraction implies that
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
Thus, 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 ,
where . Since , this equation shows that, even in the absence of growth in constituent , changes can occur in the solution volume-based concentration cσ, either as a result of volume-altering deformation (changing Js) or the evolution of from the growth of other solid constituents.
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 s = 1, implying that 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 and the transition of the stress constitutive relation from compressible to incompressible formulations when ). 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 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.
The mixture framework can describe cells and the ECM, albeit using different mixture constituents. The total stress in a mixture is given by
where p is the fluid pressure, I is the identity tensor, and 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
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 , to avoid a burdensome notation; non-ideal behavior may be easily incorporated subsequently. Thus, for water, the mechano-(electro)chemical potential is given by
where R is the universal gas constant, θ is the absolute ambient temperature (assumed uniform), and 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
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; is the chemical potential of the solute in the reference physico-chemical state when ψ = ψ0 and .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 . 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,
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
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 . No jump condition applies for interface-impermeant species. In addition, a jump condition is provided by the conservation of linear momentum for the mixture,
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 .
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, , 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.
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, . Note that in the external bath, devoid of a porous solid matrix, for all solutes (α ≠ s, w); furthermore, since there can be no bath FCD, Eq.(41) requires that .
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
Similarly, for the membrane-permeant neutral solute,
For the membrane-impermeant ECM solutes, the conditions are limited to
These expressions are arranged such that the right-hand-side depends on prescribed quantities in the external bath ,7 solubilities in the ECM and cell , the FCD in the ECM , and the intracellular concentration of membrane-impermeant solutes .
Since the charges fixed to the solid matrix of the ECM must move with it, Eq.(36) constrains the variation in according to
where 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 can be obtained from this relation via
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
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, and Eq.(52) reduces to . Furthermore, when neglecting effective stresses in the protoplasm , the jump condition of Eq.(43) reduces to pc = p*. Combining these two results implies that which, when substituted into Eq.(55), yields
where 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, , the resulting linear relationship between the relative cell volume Js and the ratio ,
is known as the Boyle-van’t Hoff relation . Most cells exhibit this linear response under osmotic loading and are called perfect osmometers . Experimental measurements of Js for various values of can be used to extract the value of the cell solid volume fraction ; note that is also known as the fraction of osmotically active water.
Alternatively, if the external bath environment remains unchanged, , but in the presence of growth of intracellular solute or solid, the cell volume ratio is given by
where we used the relation of Eq.(31). This relation shows that an increase in intracellular solute concentration , or solid content , increases the cell volume (Js > 1). This is a mathematical representation of the cell’s well-known volume regulation mechanism , 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 such that there is no fixed upper bound on according to Eq.(29).
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
Substituting Eq.(7) into this relation and recalling that is invariant produces an expression for the homogenized volume fraction of solvent in the cells and ECM,
For solutes that are not present in a particular compartment, the corresponding concentration is simply set to zero . The form of Eq.(62) can also be specialized to evaluate the homogenized FCD in the tissue, letting and yielding
The homogenized solubility is similarly obtained,
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
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
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
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 , and the solid-bound molecular species by , where . 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
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 , 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 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 , 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 . In that case the tissue would not be subjected to osmotic swelling, so that could reduce to zero under traction-free conditions ( 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 . In this case, Eq.(71) predicts that
when using Eqs.(61)–(62). Looking at this expression, it is not possible to expect p → p* as c* → ∞, unless . 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 reduce to zero only in this limiting case. Furthermore, as explained in Section (2.6), emptying a cell of its water content implies that ; however, in a homogenized model of the tissue, the lower bound on J is given by , thus it may not even be possible to empty the cell of its water content if the ECM gets depleted first .
An alternative is to propose that the reference configuration corresponds to a charged tissue bathing in a NaCl solution of finite tonicity , 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 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 reduce to zero.
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 ; the molecular species imparting the FCD to the ECM, whose content is described by ; the solid matrix of the cells, described by ; and the membrane-impermeant intracellular solute, described by . 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 , given constitutive relations for . For the FCD and intracellular membrane-impermeant solute, Eq.(36) can similarly be used to track changes in , given constitutive relations for . (Recall that may change as a result of growth or volume deformation, according to Eqs.(54) and (56), but changes in 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 , 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 . 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.
Extracellular matrix growth can be described by an increase in , 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 . Since ECM material properties may depend on the solid content , 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 represents a generic material property of the solid matrix, a constitutive relation validated from experiments may describe its dependence on . (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.
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 based on a representative composition of cartilage , and assume that the external bathing environment of the tissue is isotonic saline (2c* = 300 mOsM, ) 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 p − p* = 86 kPa. This solution is obtained by solving the jump condition for the total traction, , where p is given by Eq.(71), for any arbitrary direction n.
It is now assumed that 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 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, p − p* = 68 kPa; though was held constant, the FCD normalized to the interstitial fluid volume in the current configuration, , 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.
An increase in the FCD, , may occur by cell synthesis of charged molecules which then bind to the ECM. All else remaining the same, the resulting increase in 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.
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  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 . 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 . 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 throughout the wall thickness. The intracellular concentration of membrane-impermeant solutes is set to (a value that yields when Js = 1). The FCD is non-zero only in the intima and media; to model FCD growth, is increased from 40 mEq/L to 120 mEq/L over time. The external bath is physiological saline, with 2c* = 300 mOsM, . 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.
Intracellular concentration of membrane-impermeant solutes may change as part of the well-known mechanism for regulating cell volume ; this mechanism is primarily driven by alterations in , 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.
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 . 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” .
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 . Though brain tissue contains proteoglycans, it is assumed for simplicity that in this example, so that we can let Js = 1. Then, Eq.(71) simplifies to
Clearly, all else being equal, an increase in will lead to an increase in the intracranial pressure relative to the blood pressure, p − p*. For example, if we let χ = 0.85, (so that ), θ = 310 K, 2c* = 300 mOsM, and let the initial value yield zero pressure difference , then an increase in pressure by 10 mmHg (p − p* = 1.3 kPa) can be achieved if 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 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.
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 .
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 , and osmotic swelling of the endothelial and vascular smooth muscle cells of the intima and media . 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 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 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 . However, as decreases, the cell volume also decreases as a result of water exudation, and the opening angle becomes negative, reducing to −37° when (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.
Though this example illustrates the reduction in opening angle with decreasing , it can also be shown that the opening angle can increase with increasing . 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 , a time duration consistent with the relatively rapid process of cell volume regulation.
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 , the intracellular membrane-impermeant solute concentration, , and the cell volume fraction, χ. The evolution in 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 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
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 = ξ (X, t)/t, such that
On the assumption that the volume of daughter cells is nearly identical to the parent cell, it is also reasonable to assume that is proportional to ,
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.
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 and concentration of membrane-impermeant solute are both assumed to increase by an arbitrary factor of six at a constant rate, . 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 or cell , the osmotic pressure in the tissue is equal to the ambient bath pressure, p − p* = 0. Importantly, the intracellular solute concentration in the current configuration, , remains constant and equal to 300 mOsM during this growth process. This example also illustrates a problem where exceeds unity without violating the constraint of Eq.(29).
In this example we reprise the analysis of growth of a hyaline cartilage model, this time including an ECM. Let , and let the ECM be described by a neo-Hookean model for , 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, p − p* = 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.
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 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 , 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, , in Eq.(20). Thus, changes in only result from growth, according to the growth evolution equation (25).
Changes in 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 , 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 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 . 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 .
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.
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).
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 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 at a given θ are prescribed standard states that remain invariant .
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 f+n+ + f−n− = (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.