Areas of Interest

My academic interests fall under the following relatively broad categories:

Biological growth: Reaction, transport and mechanics

H. Narayanan, E. M. Arruda, K. Grosh, K. Garikipati

Received: date / Accepted: date

Abstract

In this paper, we address some modelling issues related to biological growth. Our treatment is based on a recently-proposed, general formulation for growth (Journal of the Mechanics and Physics of Solids, 52, 2004, 1595-1625) within the context of open system continuum thermodynamics. We aim to enhance this treatment by making it more appropriate for the biophysics of growth in soft tissue, specifically tendon. This involves several modifications to the mathematical formulation to represent the reactions, transport and mechanics, and their interactions. We also reformulate the governing differential equations for reaction-transport to represent the incompressibility constraint on the fluid phase of the tissue. This revision enables a straightforward implementation of numerical stabilisation for the hyperbolic, or advection-dominated, limit. A finite element implementation employing a staggered scheme is utilised to solve the coupled nonlinear partial differential equations that arise from the theory. Motivated by our experimental model, an in vitro scaffold-free engineered tendon formed by self-assembly of tendon fibroblasts [Calve et al., 2004], we solve several numerical examples demonstrating biophysical aspects of growth, and the improved numerical performance of the models.

1  Introduction

Growth involves the addition or depletion of mass in biological tissue. In biological systems, growth occurs in combination with remodelling, which is a change in microstructure, and possibly with morphogenesis, which is a change in form in the embryonic state. The physics of these processes are quite distinct, and for modelling purposes can, and must, be separated. Our previous work [Garikipati et al., 2004], upon which we now seek to build, drew in some measure from Cowin and Hegedus [1976],Epstein and Maugin [2000],Taber and Humphrey [2001] and Kuhl and Steinmann [2003], and was focused upon a comprehensive account of the coupling between transport and mechanics. The origins of this coupling were traced to the balance equations, kinematics and constitutive relations. A major contribution of that work was the identification and discussion of several driving forces for transport that are thermodynamically-consistent in the sense that specification of these relations does not violate the Clausius-Duhem dissipation inequality. Now, we seek to restrict the range of physically-admissible possibilities in order to gain greater physiological relevance for modelling growth in soft biological tissue. The advection-diffusion equations for mass transport require numerical stabilisation in the advection-dominated regime (the hyperbolic limit). We draw upon the enforcement of the incompressibility limit for the fluid phase to facilitate this process. Below, we briefly introduce each aspect that we have considered, but postpone details until relevant sections in the paper.
These issues are treated in detail in relevant sections of the paper, which is laid out as follows: Balance equations and kinematics are discussed in Section 2, constitutive relations for reactions, transport and mechanics in Section 3, and numerical examples are presented in Section 4. Conclusions are drawn in Section 5.

2  Balance equations and kinematics of growth

In this section, the coupled, continuum balance equations governing the behaviour of growing tissue are summarised and specialised as outlined in Section 1. For detailed continuum mechanical arguments underlying the equations, the interested reader is directed to Garikipati et al. [2004].
The tissue of interest is an open subset of R3 with a piecewise smooth boundary. At a reference placement of the tissue, Ω0 , points in the tissue are identified by their reference positions, X Ω0 . The motion of the tissue is a sufficiently smooth bijective map ϕ : Ω 0 ×[0,T] R3 , where Ω 0 := Ω0 Ω0 . At a typical time t[0,T], ϕ ( X ,t) maps a point X to its current position, x . In its current configuration, the tissue occupies a region Ωt = ϕ t ( Ω0 ). These details are depicted in Figure 1. The deformation gradient F := ϕ / X is the tangent map of ϕ .
Figure
Figure 1: The continuum tissue with growing and diffusing species.
The tissue consists of numerous species, of which the following groupings are of importance for the models: A solid species, consisting of solid collagen fibrils and cells, 1 denoted by c, an extra-cellular fluid species denoted by f and consisting primarily of water, and solute species, consisting of precursors to reactions, byproducts, nutrients, and other regulatory chemicals. A generic solute will be denoted by s. In what follows, an arbitrary species will be denoted by ι, where ι=c,f,s.
The fundamental quantities of interest are mass concentrations, ρ0 ι ( X ,t). These are the masses of each species per unit system volume in Ω0 . Formally, these quantities can also be thought of in terms of the maps ρ0 ι : Ω 0 ×[0,T]R, upon which the formulation imposes some smoothness requirements. By definition, the total material density of the tissue at a point is a sum of these concentrations over all species ι ρ0 ι = ρ0 . Other than the solid species, c, all phases have mass fluxes, M ι . These are mass flow rates per unit cross-sectional area in the reference configuration defined relative to the solid phase. Except for the fluid, f, all species have mass sources/sinks, Πι .

2.1  Balance of mass for an open system

As a result of mass transport (via the flux terms) and inter-conversion of species (via the source/sink terms) introduced previously, the concentrations, ρ0 ι , change with time. In local form, the balance of mass for an arbitrary species in the reference configuration is

ρ0 ι t = Πι - DIV [ M ι ],  ι, (1)
recalling that, in particular, M s =0 and Πf =0. Here, DIV[] is the divergence operator in the reference configuration. The functional forms of Πι are abstractions of the underlying biochemistry, physiologically relevant examples of which are discussed in Section 3.4, and the fluxes, M ι , are determined from the thermodynamically motivated constitutive relations described in Section 3.3.
The behaviour of the entire system can be determined by summing Equation (1) over all species ι. Additionally, sources and sinks satisfy the relation

ι Πι =0, (2)
which is consistent [Garikipati et al., 2004] with the Law of Mass Action for reaction rates and with mixture theory [Truesdell and Noll, 1965].

2.1.1  The role of mass balance in the current configuration

Though it is not mathematically incorrect to solve the initial-boundary-value problem in terms of reference coordinates, it is important to note that as soft tissues deform, the current configuration, Ωt , and its boundary, Ωt , change in time. Even as the pore structure at the boundary deforms with the tissue, the fluid concentration with respect to Ωt remains constant if the boundary is in contact with a fluid bath. Accordingly, this is the appropriate Dirichlet boundary condition to impose under normal physiological conditions. This is shown in an idealised manner in Figure 2.
Figure
Figure 2: If the pore structure at the boundary deforms with the tissue and this boundary is in contact with a fluid bath, the fluid concentration with respect to the current configuration, i.e., ρf , remains constant.
In order to apply boundary conditions (either specification of species flux or concentration) that are physical, it is straightforward to use the local form of the balance of mass in the current configuration,

d ρι dt = πι - div [ m ι ]- ρι div [ v ],  ι, (3)
where ρι ( x ,t), πι ( x ,t), and m ι ( x ,t) are the current mass concentration, source and mass flux of species ι respectively and v is the velocity of the solid phase. The spatial divergence operator is div[] , and the time derivative on the left hand-side in Equation (3) is the material time derivative, that may be written explicitly as t |X , implying that the reference position is held fixed.

2.2  The kinematics of growth (changes in concentration)

Figure
Figure 3: The kinematics of growth.
Local volumetric changes are associated with changes in the concentrations of species. The material of the species swells with an increase in concentration (mass of the species per unit system volume), and shrinks as its concentration decreases. This leads to the notion of a growth deformation gradient. One aspect of the coupling between mass transport and mechanics stems from this phenomenon. In the setting of finite strain kinematics, the total deformation gradient is decomposed into the growth deformation gradient, a geometrically-necessitated elastic deformation accompanying growth, and an additional elastic deformation due to external stress. This split is analogous to the classical decomposition of multiplicative plasticity [Lee, 1969] and is similar to the approach followed in existing literature on biological growth (see, for e.g., Taber and Humphrey [2001],Ambrosi and Mollica [2002]).
The split itself is visualised in Figure 3. Assuming that the volume changes associated with growth described above are isotropic, a simple form for the growth deformation gradient tensor is

F gι = ρ0 ι ρ 0ini ι 1, (4)
where ρ 0ini ι ( X ) can be interpreted as the initial reference state where the species would be stress free in the absence of a deformation, and 1 is the second-order isotropic tensor. This being a local definition, the action of F gι alone can result in incompatibility, which is eliminated by the geometrically-necessary elastic deformation F ˜ eι . Thus, the total deformation gradient F = F e F ˜ eι F gι , (where F e arises from the external stress) and internal stresses in the tissue arise due to the compatibility restoring tensor F ˜ eι .

2.2.1  Saturation and tissue swelling

Figure
Figure 4: Only under conditions of fluid saturation does an increase in fluid concentration cause tissue swelling.
Upon closer examination of the solid phase of the tissue as a porous medium, it is observed that its degree of saturation plays a fundamental role in determining whether the tissue responds by swelling or shrinking to an infusion or expulsion of fluid. In particular, the isotropic swelling law defined by Equation (4) has to be generalised to handle the case in which the solid phase is not saturated by fluid.
Figure 4 schematically depicts two potential scenarios. If the tissue is initially unsaturated (as in A), this corresponds to the fact that, on a microscopic scale, it still contains unfilled voids. It is thus capable of allowing an influx of fluid, which tends to increase its degree of saturation (to reach B), but does not cause the tissue to swell, as there exists space for incoming fluid to occupy. However, if the tissue was initially saturated (as in C), an increase in the amount of fluid will result in swelling (as depicted in D), as there is no free volume for the entering fluid to occupy. It is this second case that is modelled by (4).
The isotropic swelling law can be extended to the unsaturated case by introducing a degree of saturation arising in a natural fashion from the current concentrations, ρι . These quantities can also be thought of as the product of the intrinsic density of the species, ρ ~ ι , and the corresponding volume fraction, v ~ ι , in the current configuration. Upon solution of the mass balance equation (3) for ρι , the species volume fractions, v ~ ι , can be computed since the intrinsic densities are known material properties. The sum of these volume fractions is our required measure of saturation, and clearly cannot exceed unity in the current configuration. We thus proceed to redefine the growth deformation gradient tensor as follows:

F gι ={ 1, ι v ~ ι <1 ρ0 ι ρ 0ini ι 1,otherwise. (5)
However, it is important to note that under normal physiological conditions, soft tissues are fully saturated by the fluid and Equation (4) is appropriate.

2.3  Balance of momenta

In soft tissues, the species production rate and flux that appear on the right hand-side in Equation (1), are strongly dependent on the local state of stress. To correctly model this coupling, the balance of linear momentum should be solved to determine the local state of strain and stress.
The deformation of the tissue is characterised by the map ϕ ( X ,t). Recognising that, in tendons, the solid collagen fibrils and fibroblasts do not undergo mass transport, 2 the material velocity of this species, V = ϕ /t, is used as the primitive variable for mechanics. The motion of each remaining species is split into a deformation along with the solid species, and mass transport relative to it. To this end, it is useful to define the material velocity of a species ι relative to the solid skeleton as: V ι =(1/ ρ0 ι ) F M ι . Thus, the total material velocity of a species ι is V + V ι .
The total first Piola-Kirchhoff stress tensor, P , is the sum of the partial stresses P ι (borne by a species ι) over all the species present. 3 With the introduction of these quantities, the balance of linear momentum in local form for a species ι in Ω0 is
where g is the body force per unit mass, and q ι is an interaction term denoting the force per unit mass exerted upon ι by all other species present. The final term with the (reference) gradient denotes the contribution of the flux to the balance of momentum. In practise, the relative magnitude of the fluid mobility (and hence flux) is small, so the final term on the right hand side of Equation (6) is negligible, resulting in a more classical form of the balance of momentum. Furthermore, in the absence of significant acceleration of the tissue during growth, the left hand-side can also be neglected, reducing (6) to the quasi-static balance of linear momentum.
The balance of momentum of the entire tissue is obtained by summing Equation (6) over all ι. Additionally, recognising that the rate of change of momentum of the entire tissue is affected only by external agents and is independent of internal interactions, the following relation arises.

ι ( ρ0 ι q ι + Πι V ι )=0. (6)
This is also consistent with classical mixture theory [Truesdell and Noll, 1965]. See Garikipati et al. [2004] for further details on balance of linear momentum, and the formulation of balance of angular momentum. We only note here that the latter principle leads to a symmetric partial Cauchy stress, σ ι for each species in contrast with the unsymmetric Cauchy stress of Epstein and Maugin, [2000].

3  Constitutive framework and modelling choices

As is customary in field theories of continuum physics, the Clausius-Duhem inequality is obtained by multiplying the entropy inequality (the Second Law of thermodynamics) by the temperature field, θ, and subtracting it from the balance of energy (the First Law of thermodynamics). We assume the internal energy per unit reference volume of species ι to be of a sufficiently general form: eι = e ^ ι ( F eι , η0 ι , ρ0 ι ), where η0 ι is the entropy per unit system volume. Substituting this in the Clausius-Duhem inequality results in a form of this inequality that the specified constitutive relations must not violate. Only the valid constitutive laws relevant to the examples that follow are listed here. For details, see Garikipati et al., [2004].

3.1  An anisotropic network model based on entropic elasticity

Each constituent of the tissue has a mass-specific Helmholtz free energy density, ψι . Utilising the material response of a hyperelastic material, the partial first Piola-Kirchhoff stress of collagen is P c = ρc ψc / F ec . Here, F ec = F F g c-1 is the elastic deformation gradient, and F gc is the growth deformation gradient, of collagen. Along the lines of Equation (4), if we were considering unidirectional growth of collagen along a unit vector e , we have F gc = ρ0 c ρ 0ini c e e , with ρ 0ini c denoting the initial concentration of collagen at the point.
The mechanical response (function) of tendons in tension is determined by their dominant structural component, highly oriented fibrils of collagen. In our preliminary formulation, the strain energy density for collagen has been obtained from hierarchical multi-scale considerations based upon an entropic elasticity-based worm-like chain (WLC) model [Kratky and Porod, 1949]. The WLC model has been widely used for long chain single molecules, most prominently for DNA [Marko and Siggia, 1995,Rief et al., 1997,Bustamante et al., 2003], and recently for the collagen monomer [Sun et al., 2002]. The central parameters of this model are the chain's contour length, L, and persistence length, A. The latter is a measure of its stiffness and given by A=χ/kθ, where χ is the bending rigidity, k is Boltzmann's constant and θ is the temperature [Landau and Lifshitz, 1951]. We have fitted the WLC response function derived by Marko and Siggia [1995] to the collagen fibril data of Graham et al. [2004] with A=6 nm and L=3480 nm. This is to be compared with A=14.5 nm and L=309 nm, reported by Sun et al. [2002], for a single collagen molecule. Taken together, these results demonstrate that the WLC analysis correctly predicts a collagen fibril to be longer and slightly more compliant than its constituent molecule due to compliant intermolecular cross-links in a fibril. To model the possibility of a collagen network structure, the WLC model has been embedded as a single constituent chain of an eight-chain model [Bischoff et al., 2002a,Bischoff et al., 2002b], depicted in Figure 5. Homogenisation via averaging then leads to a continuum Helmholtz strain energy function, ψF c : 4
Here, ρ0 c is the concentration of collagen, N is the density of chains, and a,b and c are lengths of the unit cell sides aligned with the principal stretch directions. The material model is isotropic only if a=b=c.
Figure
Figure 5: The eight-chain model incorporating worm-like chains.
The elastic stretches along the unit cell axes are respectively denoted by λ1 ec , λ2 ec and λ3 ec , C ec = F e cT F ec is the elastic right Cauchy-Green strain tensor of collagen and 1 is the second-order isotropic tensor. The factors γ and β control the bulk compressibility of the model. The end to end chain length is given by r= 1 2 a2 λ1 2 + b2 λ2 2 + c2 λ3 2 , where λI = N I · C ec N I , and N I ,I=1,2,3 are the unit vectors along the three unit cell axes, respectively. In our numerical simulations that appear below in Section 4, the numerical values used for the parameters introduced in (8) have been based off those in Kuhl et al. [2005].

3.2  A nearly incompressible ideal fluid

In this preliminary work, the fluid phase is treated as a nearly incompressible, ideal, i.e., inviscid, fluid. The partial Cauchy stress in the fluid is
σ f =det( F ef )-1 P f F efT =h( ρf )1, (7)
where a large value of h' ( ρf )= dh( ρf ) d ρf ensures near-incompressibility.

3.2.1  Response of the fluid in tension; cavitation

The response of the ideal fluid, as defined by Equation (9), does not explicitly distinguish between the cases where the fluid is subjected to tension or compression, i.e., whether det( F ef )1. When the fluid phase is subjected to compression, being (nearly) incompressible, it can develop compressive stresses without bound and is modelled accurately. Under tension, the fluid can develop at most a small tensile stress [Brennen, 1995], 5 and the bulk of the tensile stiffness arises from the collagen phase.This is not accurately represented by (9), which predicts a tensile response in the fluid similar to the compressive response.
Here, we preclude all tensile load carrying by the fluid by limiting det( F ef )1. For consistency, we first introduce an additional component to the mixture, a void species, v. Denoting its deformation gradient by F v , we now only require

F v F ef = F , (8)
where F is the pointwise homogeneous deformation gradient. If det( F )<1, 6 we set F v =1, denoting a collapse of voids and a return to saturation. Otherwise, we set

det( F ef )=1 (9)
and determine det( F v ) using the relation

F v = F F e f-1 . (10)
In this manner, we allow the pores to evacuate, (i.e. allow det( F v )>1), giving us an additional measure of the unsaturation in the system.

3.3  Constitutive relations for fluxes

From Garikipati et al. [2004], the constitutive relation for the flux of extra-cellular fluid relative to the collagen takes the following form,

M f = D f ( ρ0 f F T g + F T DIV [ P f ]- GRAD [ ef -θ ηf ]), (11)
where D f is the positive semi-definite mobility of the fluid and isothermal conditions are assumed to approximate the physiological ones. Experimentally determined transport coefficients (e.g. for rat tail skin [Swartz et al., 1999] and rabbit Achilles tendons [Han et al., 2000]) are used for the fluid mobility values. The terms in the parenthesis on the right hand-side of Equation (13) sum up to give the total driving force for transport. The first term is the contribution due to gravitational acceleration. In order to maintain physiological relevance, this term has been neglected in the following treatment. The second term arises from stress divergence; for an ideal fluid, it reduces to a pressure gradient, thereby specifying that the fluid moves down a compressive pressure gradient, which is Darcy's Law. The third term can be thought of as the gradient of a chemical potential. The included entropy gradient in this term results in classical Fickean diffusion if only mixing entropy exists. For a detailed derivation and discussion of Equation (13), the reader is directed to Garikipati et al. [2004].
Since the driving forces in (13) originate from different physics, it proves useful (as seen in the following subsection) to split the mobility, D f , into permeability (stress-gradient driven) and diffusion (chemical potential gradient driven) mobilities:

M f = D P f F T DIV [ P f ]- D μ f GRAD [ ef -θ ηf ], (12)
where D P f is now the permeability of the tissue and D μ f is the mobility of the fluid phase through the porous solid.

3.3.1  Saturation and Fickean diffusion

Figure Figure
Figure 6: Only unsaturated tissues can undergo Fickean diffusion.
As depicted in Figure 6, only when pores are unsaturated are there multiple configurations available to the fluid molecules at a fixed fluid concentration. This leads to a non-zero mixing entropy. In contrast, if saturated, there is a single available configuration (degeneracy), resulting in zero mixing entropy. Consequently, Fickean diffusion, which arises from the gradient of mixing entropy can exist only in the unsaturated case. However, even a saturated pore structure can demonstrate concentration gradient-dependent mass transport phenomenologically: The fluid stress depends on fluid concentration, see Equation (9), and fluid stress gradient-driven flux appears as concentration gradient-driven flux.
The saturation dependence of Fickean diffusion is modelled by using the measure of saturation introduced in Subsection 2.2.1. We first rewrite the chemical potential gradient term in (14) as

D μ f GRAD [ ef -θ ηf ]= D e f GRAD [ ef ]- D η f GRAD [θ ηf ], (13)
where D e f and D η f are positive semi-definite mobilities, and the mobility for the Fickean contribution is

D η f ={ D ˜ η f ,   if    ι v ~ ι <1 0,otherwise. (14)
It is again important to note that under physiological conditions, soft tissues are fully saturated by fluid, and it is appropriate to set D η f =0.

3.3.2  Transport of solute species

The numerous dissolved solute species (proteins, sugars, nutrients, ...), denoted by s, undergo long range transport primarily by being advected by the perfusing fluid. In addition to this, they undergo transport relative to the fluid. This motivates an additional velocity split of the form V s = V s ˜ + V f , where V s ˜ denotes the velocity of the solute relative to the fluid. The constitutive relation for the corresponding flux, denoted by M s ˜ , has the following form, similar to Equation (13) defined for the fluid flux.

M s ˜ = D s (- GRAD [ es -θ ηs ]), (15)
where D s is the positive semi-definite mobility of the solute relative to the fluid, and again, isothermal conditions are assumed to approximate the physiological ones. Being in solution, this phase does not bear appreciable stress, and the stress divergence term is absent from the constitutive relation.

3.3.3  Objectivity and the contribution from acceleration

In our earlier treatment [Garikipati et al., 2004], the constitutive relation for the fluid flux had a driving force contribution arising from the acceleration of the solid phase, - ρ0 f F T V t . This term, being motivated by the reduced dissipation inequality, does not violate the Second Law and supports our intuitive understanding that accelerating the solid skeleton in one direction must result in an inertial driving force on the fluid in the opposite direction. However, as defined, this acceleration is obtained by the time differentiation of kinematic quantities, 7 and does not transform in a frame-indifferent manner. Unlike the superficially similar term arising from the gravity vector, 8 the acceleration term presents an improper dependence on the frame of the observer. Thus, its use in constitutive relations is inappropriate, and the term has been dropped in Equation (13).

3.3.4  Incompressible fluid in a porous solid

Upon incorporation of the additional velocity split, V s = V s ˜ + V f , described in Subsection 3.3.2, the resulting mass transport equation (3) for the solute species is

d ρs dt = πs - div [ m s ˜ + ρs ρf m f ]- ρs div [ v ], (16)
an advection-diffusion equation. In the hyperbolic limit-where advection dominates-spatial oscillations emerge in numerical solutions of this equation [Brooks and Hughes, 1982,Hughes et al., 1987]. However, the form in which the equation is obtained is not amenable for the application of standard stabilisation techniques [Hughes et al., 1987]. In part, this is because although the (near) incompressibility of the fluid phase is imbibed in the balance of linear momentum, it has not yet been explicitly incorporated into the transport equations. This subsection proceeds to impose the fluid incompressibility condition and deduces implications for the solute mass transport equation, including a crucial simplification allowing for its straightforward numerical stabilisation.
From Equation (3), the local form of the balance of mass for the fluid species (recalling that Πf =0) in the current configuration is

d ρf dt =- div [ m f ]- ρf div [ v ]. (17)
In order to impose the incompressibility of the fluid, we first denote by ρ 0ini f the initial value of the fluid reference concentration, and recognise that
which results in a very large pressure gradient driven flux due to near incompressibility. In (20), J denotes det ( F ), J fg denotes det( F gf ) and ρini f is the push-forward of ρ 0ini f to the current configuration.
Restricting the argument to a non-growing solid, i.e. F in the solid and F fe are uniform,

t ( ρ 0ini f ( X ))0 t ( ρf ( x ϕ ,t))=0, (18)
which is the hidden implication of our assumption of a homogeneous deformation, i.e., F is the deformation gradient of all phases. This leads to ρf t =0. We therefore proceed to treat our fluid mass transport at steady state. Rewriting the flux m f from Equation (19) as the product ρf v f and using the result derived above,
Thus, using the incompressibility condition (22), we get the simplified form of the balance of mass for an arbitrary solute species, s,

d ρs dt = πs - div [ m s ˜ ]- m f · grad [ ρs ] ρf + ρs m f · grad [ ρf ] ρ f2 . (19)
Using the pushed-forward form of (17), this is now in standard advection-diffusion form, and is well suited for stabilisation schemes such as the streamline upwind Petrov-Galerkin (SUPG) method (see, for e.g., Hughes et al., [1987]) described briefly below.

3.3.5  Stabilisation of the simplified solute transport equation

SUPG methods are a class of finite element methods, originally developed for the scalar advection-diffusion equation, which have proven efficient in the solution of a variety of flow problems [Hughes, 1987]. The standard form of the scalar advection-diffusion is

ϕ t + a ·grad[ϕ]=div[κgrad[ϕ]]+f, (20)
where ϕ=ϕ( x ,t) is a scalar field, a is the advective velocity, κ is a diffusivity and f is a volumetric source term. Here, if κ>0, we have the parabolic case and if κ=0, we have the hyperbolic case. In terms of numerics, the greatest challenge is posed when the element Peclet number,

α=max | a |h 2κ , (21)
is large; h being the mesh size. Spatial instabilities are observed in the numerical solution of (25) when the element Peclet number is large [Hughes, 1987].
In contrast to the standard Galerkin method, the SUPG offers greater control over the advective-derivative term by adding an upwinded, weighted residual contribution which acts only along the direction of the streamline. This stabilising control therefore maintains consistency [Hughes, 1987]. In weak form, the method for pde (25) is defined as
where Γh is the Neumann boundary, and this equation introduces a numerical stabilisation parameter τ, which we have calculated from the L2 norms of element level matrices, as described in Tezduyar and Sathe, [2003].
Figures 7 and 8 present the solute concentration contours at time t=10 s for the initial-boundary-value problem discussed in Subsection 4.3 using the standard Galerkin and an SUPG based scheme respectively.
Figure
Figure 7: Spatial oscillations observed in the solute concentration contours (at time t=10 s) obtained while solving Example 4.3 using the standard Galerkin scheme.
Figure
Figure 8: Smooth, stabilised solutions observed for the solute concentration (at time t=10 s) obtained while solving Example 4.3 using an SUPG scheme.

3.4  Nature of the sources

There exists a large body of literature, Cowin and Hegedus [1976],Epstein and Maugin [2000],Ambrosi and Mollica [2002], that addresses growth in biological tissue mainly based upon a single species undergoing transport and production/annihilation. However, when chemistry is accounted for, it is apparent that growth depends on cascades of complex biochemical reactions involving several species, and additionally involves intimate coupling between mass transfer, biochemistry and mechanics. An example of this chemo-mechanical coupling is described in Provenzano et al., [2003].
The modelling approach followed in this work is to select appropriate functional forms of the source terms for collagen, Πc , and the solutes, Πs , that abstract the complexity of the biochemistry. In our earlier exposition [Garikipati et al., 2004], we utilised simple first order chemical kinetics to define Πc . We now replace it with a source that has greater relevance from the standpoint of biochemistry.

3.4.1  Enzyme kinetics

Michaelis-Menten enzyme kinetics (see, for e.g., Sengers et al., [2004]) involves a two-step reaction with the collagen and solute production terms given by

Πs = ( kmax ρs ) ( ρm s + ρs ) ρcell ,    Πc =- Πs , (22)
where ρcell is the concentration of fibroblasts, kmax is the maximum value of the solute production reaction rate constant, and ρm s is half the solute concentration corresponding to kmax . Denoting the two-stage reaction, involving an enzyme (E), substrate (S) and a product (P), by
it can be shown that ρm s = ( k2 + k-1 ) k1 . See Bromberg and Dill [2002] for details on chemistry modelled by the Michaelis-Menten model.

3.4.2  Strain energy dependent collagen production

A strain energy dependent source term was originally proposed in the context of bone growth [Harrigan and Hamilton, 1993] and induces growth at a point when the energy density deviates from a basal value, suitably weighted by a relative density ratio. Written for collagen, it has the form

Πc =( ρ0 c ρ 0ini c )-m ψF - ψF * , (23)
where ψF is the mass-specific strain energy function, and ψF * is a reference value of this strain energy density. Equation (30) models collagen production when the strain energy density (weighted by a relative density ratio) at a point exceeds this reference value, and models annihilation otherwise.

4  Numerical examples

The theory presented in the preceding sections results in a system of non-linear coupled partial differential equations. A finite element formulation employing a staggered scheme based upon operator splits Armero, [1999,Garikipati and Rao, [2001] has been implemented in FEAP [Taylor, 1999] to solve the coupled problem. In the biphasic-fluid and collagen-problem for instance, the basic solution scheme involves keeping the displacement field fixed while solving for the concentration field using the mass transport equation. The resulting concentration field is then fixed to solve the mechanics problem. This procedure is repeated until the resulting fields satisfy the differential equations within some suitable magnitude of an error norm.
The mechanics problem is solved quasi-statically. Backward Euler is used as the time-stepping algorithm for mass transport. Non-linear projection methods [Simo et al., 1985] are used to treat the near-incompressibility imposed by water. Mixed methods, as described in Garikipati and Rao, [2001], are used for stress (and strain) gradient driven fluxes.
The following examples aim to demonstrate the mathematical formulation and aspects of the coupled phenomena as the tissue grows. The model geometry, based on our engineered tendon constructs (see Figure 9), is a cylinder 12 mm in length and 1 mm 2 in cross-sectional area. The finite element mesh used for the calculations is depicted in Figure 10.

4.1  The constriction induced growth problem

Only two phases-fluid and collagen-are included for the mass transport and mechanics. The collagen is represented by the anisotropic worm-like chain model outlined previously (see Section 3.1) and the fluid phase is modelled as ideal and nearly incompressible. The parameters used in the analysis are as presented in Table 2.1. The values chosen are representative of the kinds of biological systems we are working with. The classes of initial and boundary conditions imposed are also based on physical experiments.
Since we only have two species and we want to demonstrate growth, an "artificial" fluid sink Πf is introduced following simple first order kinetics. The collagen source will be the negative of the fluid sink: Πf =- kf ( ρ0 f - ρ 0ini f );    Πc =- Πf , where kf is the reaction rate, and ρ 0ini f is the initial concentration of fluid. When ρ0 f > ρ 0ini f , this acts as a source for collagen. The mixing entropy of fluid in the mixture with collagen is written as ηmix f =- k Mf log( ρ0 f ρ0 ), where Mf is the molecular weight of the fluid.
The boundary conditions simply corresponding to immersing the tendon in a nutrient rich bath. The initial collagen concentration is 500 kg/m 3 everywhere and the fluid concentration is 400 kg/m 3 everywhere. This is exposed to a bath where the fluid concentration is 500 kg/m 3 , so with these concentration boundary conditions set, nutrient rich fluid rushes into the tissue, and growth occurs to form more collagen. The following plots present a few results from the analysis.

4.2  A swelling problem

Figure
Figure 9: Engineered tendon constructs Calve et al., [2004].
Figure
Figure 10: The finite element mesh used.
Figure
Figure 11: The collagen concentration (kg/m 3 ) after 1800 seconds.
Figure
Figure 12: The stress (Pa) vs extension (m) curves before and after growth.
Figure
Figure 13: The volume of the tendon (m 3 ) evolving with time.
Figure shows the initial collagen concentration in the tendon. After it has been immersed in a nutrient rich bath for half an hour, the tendon has grown and the collagen concentration is now higher as seen in Figure 11. On performing a simple uniaxial tension test on the tendon before and after growth, it is observed that the grown tissue is stiffer and stronger as seen in Figure 12. Additionally, the swelling of the tendon as it is immersed in the bath takes place in two clear regimes as seen in Figure 13. There as an initial rapid swelling in a diffusion dominated regime, and a slower growth dominated swelling later on.

4.3  An enzyme-kinetics based multiphasic problem

5  Conclusion

In this paper, we have discussed a number of enhancements to our original growth formulation presented in Garikipati et al. [2004]. That formulation has served as a platform for posing a very wide range of questions on the biophysics of growth. Some issues, such as saturation, incompressibility of the fluid species and its influence upon the tissue response, and the roles of biochemical and strain energy-dependent source terms are specific to soft biological tissues. We note, however, that other issues are also applicable to a number of systems with a porous solid, transported fluid and reacting solutes. Included in these are issues of current versus reference configurations for mass transport, swelling, Fickean diffusion, fluid response in compression and tension, cavitation, and roles of permeabilities and mobilities.
These issues have been resolved using arguments posed easily in the framework derived in Garikipati et al. [2004]. The interactions engendered in the coupled reaction-transport-mechanics system are complex, as borne out by the numerical examples in Section 4. We are currently examining combinations of sources defined in Section 3.4, and aim to calibrate our choices from tendon growth experiments. The treatment of these issues has led to a formulation more suited to the biophysics of growing soft tissue, making progress toward our broader goal of applying it to applications such as the study of wound healing, pathological hypertrophy and atrophy, as well as drug efficacy and interaction.

References

[Ambrosi and Mollica 2002]
Ambrosi D, Mollica F (2002) On the mechanics of a growing tumor. Int J Engr Sci 40:1297-1316
[Armero 1999]
Armero F (1999) Formulation and finite element implementation of a multiplicative model of coupled poro-pplasticity at finite strains under fully-saturated conditions. Comp Methods in Applied Mech Engrg 171:205-241
[Bischoff et al. 2002a]
Bischoff JE, Arruda EM, Grosh K (2002a) A microstructurally based orthotropic hyperelastic constitutive law. J Applied Mechanics 69:570-579
[Bischoff et al. 2002b]
Bischoff JE, Arruda EM, Grosh K (2002b) Orthotropic elasticity in terms of an arbitrary molecular chain model. J Applied Mechanics 69:198-201
[Brennen 1995]
Brennen CE (1995) Cavitation and Bubble Dynamics. Oxford University Press
[Bromberg and Dill 2002]
Bromberg S, Dill KA (2002) Molecular Driving Forces: Statistical Thermodynamics in Chemistry and Biology. Garland
[Brooks and Hughes 1982]
Brooks A, Hughes T (1982) Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comp Methods in Applied Mech Engrg 32:199-259
[Bustamante et al. 2003]
Bustamante C, Bryant Z, Smith SB (2003) Ten years of tension: Single-molecule DNA mechanics. Nature 421:423-427
[Calve et al. 2004]
Calve S, Dennis R, Kosnik P, Baar K, Grosh K, Arruda E (2004) Engineering of functional tendon. Tissue Engineering 10:755-761
[Cowin and Hegedus 1976]
Cowin SC, Hegedus DH (1976) Bone remodeling I: A theory of adaptive elasticity. Journal of Elasticity 6:313-325
[Epstein and Maugin 2000]
Epstein M, Maugin GA (2000) Thermomechanics of volumetric growth in uniform bodies. International Journal of Plasticity 16:951-978
[Fung 1993]
Fung YC (1993) Biomechanics: Mechanical properties of living tissues, 2nd edn. Springer-Verlag, New York
[Garikipati and Rao 2001]
Garikipati K, Rao VS (2001) Recent advances in models for thermal oxidation of silicon. Journal of Computational Physics 174:138-170
[Garikipati et al. 2004]
Garikipati K, Arruda EM, Grosh K, Narayanan H, Calve S (2004) A continuum treatment of growth in biological tissue: Mass transport coupled with mechanics. Journal of Mechanics and Physics of Solids 52:1595-1625
[Graham et al. 2004]
Graham JS, Vomund AN, Phillips CL, Grandbois M (2004) Structural changes in human type I collagen fibrils investigated by force spectroscopy. Experimental Cell Research 299:335-342
[Han et al. 2000]
Han S, Gemmell SJ, Helmer KG, Grigg P, Wellen JW, Hoffman AH, Sotak CH (2000) Changes in ADC caused by tensile loading of rabbit achilles tendon: Evidence for water transport. Journal of Magnetic Resonance 144:217-227
[Harrigan and Hamilton 1993]
Harrigan TP, Hamilton JJ (1993) Finite element simulation of adaptive bone remodelling: A stability criterion and a time stepping method. Int J Numer Methods Engrg 36:837-854
[Hughes et al. 1987]
Hughes T, Franca L, Mallet M (1987) A new finite element formulation for computational fluid dynamics: VII. Convergence analysis of the generalized SUPG formulation for linear time-dependent multidimensional advective-diffusive systems. Comp Methods in Applied Mech Engrg 63(1):97-112
[Hughes 1987]
Hughes TJR (1987) Recent progress in the development and understanding of SUPG methods with special reference to the compressible Euler and Navier-Stokes equations. International Journal for Numerical Methods in Fluids 7:1261-1275
[Kratky and Porod 1949]
Kratky O, Porod G (1949) Röntgenuntersuchungen gelöster Fadenmoleküle. Recueil Trav Chim 68:1106-1122
[Kuhl and Steinmann 2003]
Kuhl E, Steinmann P (2003) Theory and numerics of geometrically-nonlinear open system mechanics. Int J Numer Methods Engrg 58:1593-1615
[Kuhl et al. 2005]
Kuhl E, Garikipati K, Arruda E, Grosh K (2005) Remodeling of biological tissue: Mechanically induced reorientation of a transversely isotropic chain network. Journal of the Mechanics and Physics of Solids 53(7):1552 - 73
[Landau and Lifshitz 1951]
Landau LD, Lifshitz EM (1951) A Course on Theoretical Physics, Volume 5, Statistical Physics, Part I. Butterworth Heinemann (reprint)
[Lee 1969]
Lee EH (1969) Elastic-Plastic Deformation at Finite Strains. J Applied Mechanics 36:1-6
[Marko and Siggia 1995]
Marko JF, Siggia ED (1995) Stretching DNA. Macromolecules 28:8759-8770
[Nordin et al. 2001]
Nordin M, Lorenz T, Campello M (2001) Biomechanics of tendons and ligaments. In: Nordin M, Frankel VH (eds) Basic Biomechanics of the Musculoskeletal System, Lippincott Williams and Wilkins, N.Y., pp 102-125
[Provenzano et al. 2003]
Provenzano PP, Martinez DA, Grindeland RE, Dwyver KW, Turner J, Vailas AC, Vanderby R (2003) Hindlimb unloading alters ligament healing. Journal of Applied Physiology 94:314-324
[Rief et al. 1997]
Rief M, Oesterhelt F, Heymann B, Gaub HE (1997) Single Molecule Force Spectroscopy o Polysaccharides by Atomic Force Microscopy. Science 275:1295-1297
[Sengers et al. 2004]
Sengers BG, Oomens CWJ, Baaijens FPT (2004) An integrated finite-element approach to mechanics, transport and biosynthesis in tissue engineering. J Bio Mech Engrg 126:82-91
[Simo et al. 1985]
Simo JC, Taylor RL, Pister KS (1985) Variational and projection methods for the volume constraint in finite deformation elasto-plasticity. Comp Methods in Applied Mech Engrg 51:177-208
[Sun et al. 2002]
Sun YL, Luo ZP, Fertala A, An KN (2002) Direct quantification of the flexibility of type I collagen monomer. Biochemical and Biophysical Research Communications 295:382-386
[Swartz et al. 1999]
Swartz M, Kaipainen A, Netti PE, Brekken C, Boucher Y, Grodzinsky AJ, Jain RK (1999) Mechanics of interstitial-lymphatic fluid transport: Theoretical foundation and experimental validation. J Bio Mech 32:1297-1307
[Taber and Humphrey 2001]
Taber LA, Humphrey JD (2001) Stress-modulated growth, residual stress and vascular heterogeneity. J Bio Mech Engrg 123:528-535
[Taylor 1999]
Taylor RL (1999) FEAP - A Finite Element Analysis Program. University of California at Berkeley, Berkeley, CA
[Tezduyar and Sathe 2003]
Tezduyar T, Sathe S (2003) Stabilization parameters in SUPG and PSPG formulations. Journal of Computational and Applied Mechanics 4:71-88
[Truesdell and Noll 1965]
Truesdell C, Noll W (1965) The Non-linear Field Theories (Handbuch der Physik, band III). Springer, Berlin

Footnotes:

1 At this point, we do not distinguish the solid species further. This is a good approximation to the physiological setting for tendons, which are relatively acellular and whose dry mass consists of up to 75% collagen [Nordin et al., 2001].
2 Currently, we do not consider physiological processes which involve migration of the cells or surrounding matrix within the tissue-such as the migration of fibroblasts within the extra-cellular matrix during wound healing.
3 The amino acids, nutrients and regulators are in solution at low concentrations, and do not bear any appreciable stress.
4 Under isothermal conditions, the only contribution to ψc is from the strain energy.
5 Where we are referring to the fluid being subject to pure tension, not a reduction in fluid compressive stress from an initial ambient pressure, manifesting itself as tensile stress.
6 The argument above works for all cases other than the one where the system starts off unsaturated and the applied deformation has determinant <1, necessitating a collapse of voids. The fix is probably right here. Something like det( F )< 1 1+ v ~ ini v .
7 and not in terms of acceleration relative to fixed stars for e.g., as discussed in Truesdell and Noll, [1965].
8 where every observer has an implicit knowledge of the directionality of the field relative to a fixed frame, allowing it to transform objectively.

(Will be populated when I have the time.)

Material forces and remodelling in biological tissue

Diffusion studies in Silicon

Emission studies in engines modified to LPG