PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Soft Matter. Author manuscript; available in PMC 2017 September 28.
Published in final edited form as:
PMCID: PMC5341105
NIHMSID: NIHMS812169

Tuning Tissue Growth with Scaffold Degradation in Enzyme-Sensitive Hydrogels: a Mathematical Model

Abstract

Despite tremendous advances in the field of tissue engineering, a number of obstacles remain that hinder its successful translation to the clinic. One challenge that relates to the use of cells encapsulated in a hydrogel is identifying a hydrogel design that can provide an appropriate environment for cells to successfully synthesize and deposit new matrix molecules while providing a mechanical support that can resist physiological loads at the early stage of implementation. A solution to this problem has been to balance tissue growth and hydrogel degradation. However, identifying this balance is difficult due to the complexity of coupling diffusion, deposition, and degradation mechanisms. Very little is known about the complex behavior of these mechanisms, emphasizing the need for a rigorous mathematical approach that can assist and guide experimental advances. To address this issue, this paper discusses a model for interstitial growth based on mixture theory, that can capture the coupling between cell-mediated hydrogel degradation (i.e., hydrogels containing enzyme-sensitive crosslinks) and the transport of extracellular matrix (ECM) molecules released by encapsulated cells within a hydrogel. Taking cartilage tissue engineering as an example, the model investigates the role of enzymatic degradation on ECM diffusion and its impact on two important outcomes: the extent of ECM transport (and deposition) and the evolution of the hydrogel's mechanical integrity. Numerical results based on finite element analysis show that if properly tuned, enzymatic degradation yields the appearance of a highly localized degradation front propagating away from the cell, which can be immediately followed by a front of growing neotissue. We show that this situation is key to maintaining mechanical properties (e.g., stiffness) while allowing for deposition of new ECM molecules. Overall, our study suggests a hydrogel design that could enable successful tissue engineering (e.g., of cartilage, bone, etc.) where mechanical integrity is important.

1. INTRODUCTION

The degeneration of tissue such as cartilage and bone due to injury, aging or disease is a major source of disability, pain and economic burden in the U.S. and throughout the world. Current solutions range from prosthetic replacement to the autograft/allografts that replace the affected tissues, and thus alleviating the pain for a limited time. Tissue engineering presents a promising alternative whereby the damaged tissue is injected with a population of cells carried in a three-dimensional material that can stimulate repair mechanisms and restore normal function. Early successes have achieved skin regrowth (1), and recent progress in stem cell translation raises hope for the development of personalized strategies that uses the patients’ own stem cells without requiring a tissue biopsy to obtain tissue specific cells (2). Despite these early successes, a number of obstacles still hinders applications to a broader population with consistent outcomes. A reason is that tissue growth requires multi-functional scaffolds that can simultaneously provide a mechanical support to the nascent tissue, convey adequate biological signals to implanted cells, enable mass transfer and degrade in time to provide sufficient space for tissue development (3).

Hydrogels have proven to be promising systems due to their controllable properties mimicking that of native tissues (4,5), their ease of (non-invasive) implantation by injection, as well as their ability to support cell encapsulation and promote extracellular matrix (ECM) synthesis (6,7). They also have an enormous potential for functionalization in terms of adhesion, microarchitecture and degradation to guide the behavior of encapsulated cells. Tuning the design of these complex hydrogels has been a significant challenge because of our limited understanding, and thus control, of how hydrogel structure and cell response evolve during growth. In particular, contrary to more classical tissue engineering in porous scaffolds (8), a challenging issue for injectable hydrogel-based techniques is to ensure sustained tissue growth while maintaining the scaffold's structural integrity during the development of engineered cartilage (3). Indeed, on the one hand, nondegradable hydrogels possess a network of cross-links that inhibit the diffusion of most ECM molecules (9,10) and restrict tissue development to the immediate region surrounding the cell (1113). One the other hand, degradable cross-links may solve the problem in the short-term, but ultimately leads to the loss of the hydrogel's load carrying capacity (and thus construct failure) before the neotissue is formed. Solutions have been suggested to address this issue, such as introducing cell-mediated degradation (14,15) but they often make the design more complex and harder to predict without theoretical guidance.

In contrast to most engineering materials, a major hindrance in hydrogel design has been the lack of theoretical and computational developments to assist and guide experimental efforts. One reason is that in tissue engineering, hydrogels are not designed for their instantaneous properties, but for the temporal evolution of their structure, which interact with a biologically active component (cells). A fundamental question can therefore be cast as follows: Is it possible to finely tune the structure and degradation properties of a cell-laden hydrogel in order to enable tissue growth and continuous mechanical integrity during its transition to native tissue? To explore this problem, an interdisciplinary approach was taken; first, theoretical models of tissue growth, accounting for mass transport and reaction kinetic, must be integrated with physical model of hydrogels and their degradations. Second, these models must be combined with micro-macro analysis (sometimes called multi-scale analysis), in order to identify how phenomena occurring at the scale of the cell and polymer structure yield an emerging behavior at the scale of the construct (scaffold + tissue). Third, these models must be validated and integrated with experiments in the field of tissue engineering. While the literature on growth models is rich (1624), research on the interaction of growth mechanisms and the mechanics and degradation of a polymer scaffold has been more elusive. Generally, the scaffold's function is to provide a mechanical support for cells and a material onto which new material can be deposited. On the other hand, degradation not only is necessary to ensure the transition from construct to tissue but it can also be used to facilitate the transport and macromolecules that later becomes the new tissue. Hydrogel design can therefore direct the temporal and spatial evolution of swelling, degradation and therefore controls the way by which new tissue grows and acquires its mechanical properties (25,26). Few researchers have investigated these processes, through modeling (25, 2729) at the cell-hydrogel level. In this context, Dhote et al. (26) built a single cell model under the centro-symmetry assumption and showed that localized degradation of the encapsulating scaffold helps maintain the mechanical integrity of the construct. Sengers et al. (29) investigated the competition among ECM degradation or deposition and transport using a 2D model and its effect on the overall construct stiffness. Trewenack et al. (27) proposed a multispecies formulation of cell-mediated growth in cartilage constructs, pointing out the distinct roles of advection of diffusion fluxes at the microscopic level. Finally, Haider et al. (28) extended the phenomenological model of Wilson et al. (30) to incorporate experimentally measurable quantities, including the apparent densities of the scaffold, the deposited and the unlinked ECM. Although these models provide a qualitative understanding hydrogel supported-growth, efforts are required to make simulation more realistic and closer to experimental systems so that they become predictive and guide the current efforts in designing appropriate scaffolds that can be translated to the clinical.

In this paper, we investigated the possibility of growing a tissue in a cell-laden hydrogel through mathematical and numerical modeling. Without lack of generality, we focus our study on the case of cartilage cells (chondrocytes) embedded in a PEG-hydrogel. Due to their small mesh-size (few nanometers) compared to ECM molecules, growth in these systems has typically been hindered when degradation does not occur. The introduction of hydrolytic (or bulk) degradation has however not solved the problem since the scaffold loses its mechanical integrity before growth can occur. We investigate here how both issues can be resolved by using a scaffold that can be locally degraded by enzymes produced by cells. For this, we build a multiphasic model of the growing construct that eventually materializes through a system of nonlinear reaction-diffusion equations whose key parameters can be tuned with hydrogel design. We show, through three-dimensional numerical simulation of tissue development and scaffold degradation around a cells that growth and sustained mechanical integrity can be achieved in a fairly narrow region of the gel design parameters. A discussion of these results, their impact on the field of tissue engineering and the motivation for new research direction are then provided.

2. GROWTH IN ENZYME DEGRABLE HYDROGEL SCAFFOLD

Tissue engineering starts with the encapsulation of chondrocytes (cartilage cells) in a hydrogel, followed by their local synthesis, diffusion and deposition of ECM molecules; in most cases, these processes occur under physiological loads so that their evolution into a fully functional tissue depends on the continued mechanical integrity of the scaffold. A better understanding of how scaffold design affects mechanics can be gained by investigating the evolution of the ECM/scaffold architecture at the cellular scale, and particularly the percolation of ECM into a viable neo-tissue. We concentrate here on a particular model of degradation, or enzymatic degradation, in which cell-mediated enzymes can diffuse and locally disrupt cross-links in the gel. More specifically, we concentrate on one example system, an enzyme-sensitive poly(ethylene glycol) hydrogel made by photo-polymerization of 8-arm PEG functionalized with norbornene and enzyme sensitive peptide crosslinks (Fig. 1). The hydrogel's elastic modulus can be adjusted by changing the molecular weight of the monomers or formulation (31), while its degradation kinetics can be controlled by changing the amino acids in the peptide (32). As shown in Fig. 2, cell-mediated scaffold degradation and tissue growth can be summarized in a few steps. (a) The encapsulated cells release both enzyme and ECM molecules through their cell membrane. (b) Due to their small size, enzymes can diffuse through the polymer network and cleave existing cross-links following Michaelis-Menten kinetics (32). The change of cross-link density ρ, defined as the number of cross-links per volume of hydrogel in the dry state (i.e. when no solvent is present) can therefore be written:

DρDt=kceρρ+k
(2.1)

where D/Dt represents the material time derivative, ce is the enzyme concentration while k′ = kcat(J0 – 1) is the degradation rate constant, k″ = Km (J0 – 1) is the cross-link density at which the reaction rate is at half-maximum, while kcat and Km are Michaelis-Menten kinetic constants (32,33), and J0 is the initial swelling ratio after cell encapsulation, but before degradation occurs. Cross-link density will eventually reach the critical (or reverse-gelation) point ρc at which the polymer network loses its overall connectivity. In other words, the polymer locally loses its solid characteristics and its resistance to molecular diffusion. (c) The removal of cross-links enables the transport of large (unlinked) ECM molecules that can eventually link to form a new tissue. (d) A new tissue is formed when localized pockets of linked ECM percolate and are able to sustain mechanical loads without the support of the hydrogel. It is important to note that the above scenario is difficult to achieve experimentally since the point of reverse gelation must be closely matched to point when ECM percolation occurs. When reverse gelation occurs before ECM percolation, the result is a mechanical collapse (or dissolution) of the hydrogel and an overall loss in mechanical competence of the construct (34).

Fig. 1
Schematic of a enzymatically-degradable hydrogel made by photo-polymerization of monomers with enzyme sensitive peptide crosslinkers.
Fig. 2
Cell mediated tissue growth in a porous scaffold. (a) The construct is composed of cells and the mixture of solid and fluid phases which are either secreted by cell or key elements of a hydrogel. (b) Tissue growth in a hydrogel at the microscale is shown. ...

As stated above, one of the key mechanism underlying tissue growth in a hydrogel is molecular transport, and particularly the transport of enzyme and ECM molecules. The diffusivity of a molecule (denoted here by the index α) depends on the ratio rα of its hydrodynamic radius rα and the hydrogel mesh size ξ where the latter is defined as the distance between two opposite polymer chains in the swollen ideal network (Fig. 1). One well accepted model to describe this relationship was derived by Lustig and Peppas (35) and displays a linear relationship between the diffusivity and the relative particle radius rα of the form:

Dα(rα,ρ,J)=Dαf(J,J0)(1rαξ(ρ,J));rα<ξ<ξc
(2.2)

where the mesh size ξ implicitly depends on cross-link density ρ and equilibrium volumetric swelling ratio J and ξc = ξ(ρc, J) is the mesh size when the gel reaches reverse gelation. The function f(J, J0) represents the probability for a solute molecule to find a free volume for diffusion. Its expression can be derived from Eyring's equation (37) in which the diffusivity in a swollen network differs from that in a pure solvent due to the entropic contributions of polymer chains. Cohen and Turnbull (38) and Peppas and Reinhart (39) particularly showed that a good approximation of this probability is given by f(J, J0) = exp(−1/(JJ0 – 1)). Note that once the hydrogel reaches reverse gelation, polymer chains are free to diffuse and assumed to quickly leave the construct. This means that as ρρc, chains eventually disappear, and since this process is relatively faster than the dynamics of growth, one assumes that the function f(J, J0) becomes unity as soon as ρ = ρc. We further note that relation (2.2) is only valid over the defined range rα < ξ < ξc; indeed, when ξ < rα, the small mesh size restrict diffusion and Dα → 0. When ξ > ξc, however, the mesh disappears (ξ) and the diffusivity becomes that of a particle in a pure solvent, i.e. Dα=Dα as given by the Stokes-Einstein relation (40). While the radii of molecules cannot be changed, it is possible to tune the hydrogel mesh size via its cross-link density through the relation (35):

ξ(ρ,J)=J13ρ0ρ;ρc<ρ
(2.3)

where J is the swelling ratio of the gel from its dry state and the lumped length parameter [ell] (unit of Angstrom) is given by

=l3ρpCnMrρ0
(2.4)

Here, Cn is the characteristic ratio that defines the restricted rotation of a single polymer chain, l is the average bond length in the polymer, Mr is the molar weight of the polymer repeat unit and ρp is the density of the polymer. We note that equation (2.3) is only valid before reverse gelation and diverges to infinity as ρρc. Two observations can be made based on equations (2.2) and (2.3). First, enzymes, due to their small hydrodynamic radius (~60-85 Å (33)) may enter the hydrogel space but their diffusion can be strongly hindered for tightly cross-linked gels. ECM molecules (in particular collagen and aggrecans) are however characterized by relatively large sizes, on the order of 200 Å and larger, and are therefore unable to be transported within the gel. Second, as the gel reaches reverse gelation, the diffusivity of ECM suddenly reaches a finite value, allowing for matrix transport. This means that diffusion of ECM molecules is only possible once the gel is fully degraded, a phenomenon that has been a major hindrance to successful tissue growth. Using equations (2.2) and (2.3), Fig. 3a shows the dependency of the normalized diffusion constant DαDα on cross-link density and its evolution during degradation for both enzyme and ECM. One can clearly see the sharp rise in diffusivity at normal gelation for large ECM particles in contrast to smaller enzymes. Fig. 3b confirms this trend by showing histological images of ECM (aggrecan and collagen) around the cartilage cells in two non-degradable gels characterized by their polymer fraction (10 and 20%, respectively) at the time of hydrogel formation. It is evident that due to their inability to diffuse, these molecules remain confined in the pericellular space. One notes that a stiffer gel (20% polymer fraction) provides a stronger confinement and thus a smaller pericellular matrix.

Fig. 3
(a) Relationship between particle diffusivity and cross-link density for two different particle sizes (enzyme and ECM molecules), characterized by the ratio r/ξ. (b) Restricted diffusion of two ECM building blocks (aggrecan and collagen) around ...

3. MODEL FORMULATION

To explore the interactions between molecular transport, polymer degradation and tissue growth, we propose here to develop a mathematical model based on (multiphasic) mixture theory (42,43). More precisely, the construct is described as a continuous mixture of fluid and solid constituents whose respective volume fraction are denoted, respectively by ϕ and ϕ~. Solid constituents including the polymer (ϕ~p) and newly deposited matrix (ϕ~m) primarily participate to the mechanical integrity of the construct. Fluid-like constituents, including water (ϕw), enzymes (ϕe) and freshly synthetized and unlinked matrix molecules (ϕm), play an important role in degradation and transport. We note here that the latter two constituents are typically found in very small proportions, such that ϕm << 1 and ϕe << 1. We further assume that the above constituents constitute the majority of the inter-cellular space, yielding the saturation condition:

ϕ~m+ϕ~p=ϕ~,ϕwϕ,ϕ+ϕ~=1
(3.1)

where we omitted to specify the argument (x, t) for the sake of clarity. It is however clear that these quantities vary both in time and space during the growth process. We show next that the evolution of these quantities in time and space, together with the mechanical integrity of the mixture can be described in terms of a coupled system of reaction-diffusion equations and force balance over the physical domain of the construct.

3.1 Mass transfer

As the construct composition evolves, motion and mass transfer occurs in the inter-cellular space. In this dynamic process, all solid constituents typically follow the same displacement, denoted here as u (x, t) from their reference location x during cell encapsulation. This concept allows us to consider a more amenable mathematical treatment known as the constrained mixture formulation (20). Mass transfer within the construct may thus be described in terms of the velocities vα(x, t), (α = w, e, m) of each fluid-like constituent relative to the solid constituents. This description therefore attaches our point of observation to the solid phase, such that the total velocity of a fluid phase in a fixed frame is v + u, where u denotes the solid velocity. Assuming that all constituents are incompressible, the continuity equation for the solvent takes the form (44):

u.(ϕvw)=0
(3.2)

where [nabla] · u is to be interpreted as the rate of solid swelling. Noting that the length scale of our analysis (on the order of microns) is significantly larger than the pore scale (on the order of nanometers), the fluid flow through the polymer network typically falls within the Darcy regime (4547). This means that the fluid velocity can be related to the water pressure π and the osmotic pressure π by ϕvw = −κ [nabla] (pπ) (48,49). The continuity equation for cell secreted molecules further takes the form of a convection-diffusion-reaction system of the form:

DcαDt=(Dαcα)(vwcα)cαu.Γαα=e,m
(3.3)

We recall here that the diffusivity Dα is a nonlinear function of cross-link density and swelling ratio as shown by equations (2.2) and (2.3). Furthermore, the molar concentrations are related to volume fractions by cα = ϕα/vα with vα the molar volume of constituent α, all of which remain constant due to the incompressibility assumption. The first term is purely diffusive, the second is the convective term from the motion of the solvent, the third describes a change of concentration due to polymer swelling and the last is a reaction term that accounts for the rate of production or consumption of the constituent. We assume here that enzyme only affects gel degradation and does not transform during the process, i.e. Γe = 0. By contrast, unlinked ECM molecules do transform into linked ECM molecules and this process is described by the consumption rate Γm. Mass balance of solid constituents (polymer and linked matrix) may similarly be derived by cancelling the diffusion and convection terms in equation (3.3), which leads to:

Dc~αDt=c~αu.+Γα,α=p,m
(3.4)

where again, the only non-vanishing reaction term is the transformation rate Γm between unlinked and linked ECM. Motivated by the concept of product inhibition hypothesis (50), the state that a biological system tends to evolve towards a specific target, we propose a model in which the rate of ECM linkage is regulated by the level of existing linked ECM as follows:

Γm=km(1c~mc~m0)cm
(3.5)

where c~m0 is the target ECM concentration and km is the intrinsic rate of linkage (i.e., the rate constant in the absence of deposited ECM (cm = 0).)

3.2 Evolution of construct mechanics

As discussed above, an important aspect of hydrogel assisted growth is the evolution of mechanical properties, and especially the transition from a cell-laden hydrogel to a tissue. Since the construct possesses a heterogeneous structure during growth, this transition depends on both the spatial distribution and the load carrying capacity of each phase (hydrogel and linked ECM) in time. To characterize this evolution, let us consider that at a given time, the construct is subjected to a macroscopic deformation, which at the level of the intercellular space, yields a heterogeneous displacement field u(x, t) around cells. Material deformation can therefore be characterized by the deformation gradient F = I + [partial differential]u/[partial differential]x, or alternatively the nonlinear strain ϵ = (1/2)(FTFI), where I is the identity tensor. The change in material volume (or swelling ratio) from its initial equilibrium value prior to any degradation is then measured by J = det(F). For an elastic material, strains give rise to the appearance of stresses, represented by the tensor σ that must verify the balance of linear momentum:

(σpI)=0andσ=ψFFTJ
(3.6)

More precisely, σ is the elastic stress supported by the hydrogel and the newly deposited ECM while p is the internal pressure in the solvent. When no external loads are applied to the construct, this pressure is equal to the osmotic pressure π in the hydrogel (24,51):

p=RTvw(ln(ϕ)+(1ϕ)+χ(1ϕ)2)
(3.7)

The elastic stress is further related to the deformation gradient through a strain energy function, which we assume to be additively decomposed (52) into polymer and ECM components as:

ψ(ρ,cm,ϵ)=ψp(ρ,ϵ)+ψm(C~m,ϵ)
(3.8)

The energy ψp stored in the polymer is expressed in terms of the classical Rubber elasticity theory as (51):

ψp(ρ,F)=12ρRT[J023tr(FTF)3ln(J0J)],ρ>ρc
(3.9)

where R and T are the gas constant and the absolute temperature, respectively, and J0 is the initial equilibrium swelling ratio of the hydrogel prior to degradation. The latter term implies that the polymer is in a stressed state, balanced by the osmotic pressure of the solvent, equation (3.7), after cell encapsulation and equilibration with the solvent. The consequence of this phenomenon is that the stored elastic energy remains positive, even in the absence of deformation F. One can note from equation (3.9) that a drop in cross-link density ρ tends to decrease the elastic energy stored in the gel and hence its load carrying capacity. Importantly, when the gel reaches reverse gelation (ρ = ρc), it loses all of its elasticity and:

ψp(ρ,F)=0,ρ<ρc
(3.10)

This behavior creates a sharp discontinuity in the material's response as the hydrogel undergoes reverse gelation; this feature is most likely responsible for the sudden failure of the polymer network before a functional tissue can be obtained. Unlike the hydrogel, the elastic energy stored in the newly deposited matrix has not been well characterized (20,22,24). However, for convenience and consistency with previous work (24), we choose here to express it in terms of a Saint-Venant strain energy function, that is a generalization of Hooke's law in the context of finite deformation. The corresponding strain energy function reads:

ψm(ϵ)=C~m[λ2[tr(ϵ)]2+μtr(ϵ2)]
(3.11)

where λ and μ are the Lame constants, expressed in units of energy per mole of linked ECM. This definition of material parameters implies that the stiffness of the material increases the apparent density of linked ECM (28). In this study, we introduce the constraint λ = 0.8μ in order to satisfy the porous network behavior (53, 54) of linked ECM similar to the that of polymer. Furthermore, since the deformation of the hydrogel and neo-tissue remains relatively small during the growth process investigated in this work, the stress-strain relations in (3.6) can be simplified according to the small strain theory for which one can use the approximation FTF, ≈ FT + FI, ϵ = 1/2(FT + F) – I and J = tr(ϵ + I). The total elastic stress at any point within the extra-cellular space σ = σp + σm can therefore be written as:

σp=ρH(ρρc)RTtr(ϵ+I)(J023(2ϵ+I)12I)
(3.12)
σm=C~m(λ[tr(ϵ)]I+2μϵ)
(3.13)

In which the Heaviside function H(ρρc) ensures that the polymer stress vanishes when ρ = ρc.

4. MODEL ANALYSIS AND SOLUTION

In this study, we concentrate on the conditions by which the competition between hydrogel degradation and ECM deposition enables a smooth transition between scaffold and tissue properties. As shown in Fig. 4, for a homogeneous distribution of cells at small density, the analysis simplifies to the problem of a periodic representative volume of characteristic length L, centered on a single cell of radius Rc. The average spacing L between cells is further related to their volume fraction f in the construct as:

L=Rc4π3f3
(4.1)

To specify proper assumption and boundary conditions, one must first specify the external conditions surrounding the construct during the growth process. We first consider that the specimen is free of external loads (neither constant nor cyclic) at any time; consequence of this is that (a) external stresses on the external boundary of the represented domain vanish and (b) solvent permeation becomes negligible, i.e. vw = 0. This means that the convective term in equation (3.3) is negligible and the transport of enzyme and ECM is described by a coupled system of diffusion-reaction equations. Further assuming that the hydrogel reaches its new equilibrium significantly faster than the dynamics of degradation and enzyme diffusion, the osmotic pressure can be shown to follow expression (3.7). Based on the periodicity assumption of the representative domain, boundary conditions (B.C.) for g = {u, p, ce, Cm} become (55):

gS+gS=Lg
(4.2)

where S+and S are two opposite boundaries of the representative domain while g is macroscopic gradient of the field g over the domain. In particular, the absence of loading imposes that p=Cα=0, while u is directly related to the overall strain experienced (such as swelling) by the domain during growth. This quantity is calculated by ensuring stress-free boundary conditions as described in (56). Boundary conditions on the cell surface Scell can further be specified in terms of enzyme and ECM productions as:

qαScell=qα0(1cαcα0)
(4.3)

where qα0 is the production rate of constituent α (= e, m). This simple relation indicates that cells produce enzyme and ECM molecule until they reach a homeostatic state, given by a target concentration cα0 (28).

Fig. 4
Model domain and boundary conditions.

To analyze the solution, we first apply the technique of non-dimensionalization that consists of rescaling all variables to reference dimensions and times as:

x=xL,t=tDmL2
(4.4)

This operation aims to simplify the original equations and reduce the number of physically meaningful parameters for a better interpretation of the system's response. Using the above scaling relations, the new dimensionless variables become:

u=uL,cα=cαcα0,p=pRT,ρ=ρρ0
(4.5)

Where the characteristic length scale L is determined by the average spacing between cells, energies are scaled with the thermal energy RT and all concentrations are defined in relation to their threshold value cα0 defined in equation (4.2). Substituting these variables in the original formulation, our final system preserves the structure of the original transport and mechanics equations (3.3) and (3.6), for which the critical parameters are expressed as follows. For the transport phenomena, normalized permeabilities are written as a fraction of the maximum value Dm for ECM molecules and are expressed in terms of swelling ratio, molecule radius and cross-link density as:

Dα(J,rα,ρ)=ΔαJ13exp(1JJ01)(J13=rαρ)
(4.6)

with Δα=DαDm and

rα=rα
(4.7)

It is clear here that rα is interpreted as the relative size of an enzyme or ECM molecule compared to the hydrogel mesh size. Similarly, the normalized rate of degradation and rate of ECM deposition read:

κe=kc0eDmρ0L2,κm=kmDmL2
(4.8)

They characterize, respectively, the competition between enzyme degradation and diffusion and the competition between matrix deposition and diffusion. The flux boundary condition at the cell boundary is similarly normalized as qα0=qα0L2(Dmcm0) and the properties of the ECM are written in relation to the polymer stiffness and RT as:

μ=Cm0μRT,λ=λμ
(4.9)

where Poisson's ratio is v = 0.5(1 + λ*). Note here that the non-dimensional critical cross-link density ρc is also known as the network connectivity parameter β(33,57). This parameter is important to the degradation process as decrease in its value would imply an increase in the number of crosslink to be cleaved before reverse gelation. The parameters used in the non-dimensionalization, along with their value and corresponding references are listed in Table 1. These values are key as they enable the mapping of all non-dimensional results obtained in the next section into physical quantities.

Table 1
Parameters used in non-dimensionalization

5. RESULTS

The above equations were solved numerically using a nonlinear finite element scheme, whose details are given in appendix. Our approach relies on investigating the spatial and temporal evolution of the hydrogel and ECM locally around cells and relating this to the temporal evolution in the overall construct mechanical properties. For clarity, we decompose our approach in three steps; in the first example, we investigate the role of hydrogel design on the local hydrogel degradation dynamics, without transport of ECM molecules. The second example then investigates how hydrogel degradation and transport of ECM molecules interact to enable localized growth. We finally explore how various ECM growth/hydrogel degradation dynamics influence the mechanical integrity of the construct in time. This exercise allows us to identify a region in the hydrogel design that enables an optimized combination of ECM growth and mechanical integrity in time.

5.1. Degradation around a Single Cell

We first investigate the coupled mechanisms of enzyme transport-hydrogel degradation-enzyme diffusion without ECM production. For this problem, the relevant non-dimensional parameters can be reduced to the relative enzyme size r~e, and the enzymatic relative degradation κe. The effect of polymer mesh size and cross-link density are implicitly contained in these quantities (equations (4.6) and (4.7)). For the sake of illustration, we consider a construct with low cell density (L/rc ≈ 20), β = 0.8 (see Table - 1) and a constant enzyme production equal to Pe0=6×1016 mole/(cell · week) (58).

Degradation front

We show in Fig. 5 that, according to the value of parameters re and κe, the nonlinear diffusion-reaction equation yields an enzyme concentration profile that ranges from a diffusion-like appearance to a more narrow and propagating wave. The corresponding evolution in cross-link density exhibits similar features (Fig. 5 (a,b,c)) between the fully degraded region (ρ* = β) and the intact region (ρ* = 1). Importantly, we note that this transition region or “degradation front” travels away from the cell surface in time. In this study, this front is characterized in two ways: its speed is defined as that of the boundary (ρ* = β) that separates the completely degraded and non-degraded gel regions while its width is defined as the distance between the boundary (ρ = ρc) and the point at which ρ* = 0.99. Similar observations are reported both numerically and experimentally in (33).

Fig. 5
Characteristics of degradation dynamics as a function of the normalized enzyme size re and degradation rate κe. Plots (a,b,c) show the distribution of cross-link density in terms of the distance × from the cell surface ...

Diffusion-dominated and diffusion-limited systems

Figs. 5 (a,b,c) clearly shows that the width and speed of the degradation front, can be controlled by varying κe and r~e. For large enzyme radius (or small hydrogel mesh size), the diffusivity becomes so small that it is the rate-limiting step. Such a diffusion-limited system typically exhibits a very sharp degradation front followed by a region of intact polymer (Fig. 5c). By contrast, when the enzyme size becomes small (or alternatively, when the hydrogel mesh size becomes large), the diffusivity is close to that experienced in a pure solvent and hydrogel degradation becomes the rate limiting step. In this case, the system is diffusion-dominated and the enzyme concentration and cross-link density both display a diffusion-like profile away from the cell surface (Fig. 5a).

Characterization of the degradation dynamics

At first sight, our results seem to imply that sharp fronts are faster than their wider counterpart. This can be seen in Fig. 5a-c by comparing the evolution of the point ρ* = 0 (shaded regions) for three gels exhibiting wide (I), intermediate (II) and sharp (III) fronts, respectively. To better understand this trend, we further investigated the role of the relative enzyme size and degradation on the width and speed of the degradation front. This was done by numerically scanning the space (r~e,κe) and estimating the values of front width (w) and velocity (v) for each simulation. Results, recorded in the form of maps (Fig. 5d,e) clearly show nonlinear relationships between the hydrogel design and the degradation behavior, but generally suggest that increasing r~e and κe results in a sharper and faster degradation front. We note here that these results cannot be generalized as different trends may be observed for different values of the network connectivity β (33). A more general understanding of this system may be found by studying the role of β on the width and speed of the degradation front, yielding three-dimensional maps in the (re, κe, β) space. This is however, beyond the scope of this study. Overall these results indicate that it is possible to tune the hydrogel design and specifically its mesh size, cross-link density as well as the sensitivity of degradable links to enzymes to yield a variety of localized degradation dynamics around cells. The way by which such dynamics influence growth is discussed next.

5.2. The role of hydrogel degradation on the nature of ECM deposition

In this work, tissue growth is defined as cell-mediated (a) synthesis, (b) release, (c) transport and (d) deposition of ECM molecules within a hydrogel. Since the size of these macro-molecules (whether it is collagen, GAGs or other key constituents of native tissue) is usually a few times larger that of the hydrogel mesh size, their diffusion is strongly hindered in intact hydrogels as can be seen in Fig. 3. This means that growth is limited in non-degradable hydrogels. However, degradable hydrogels may enable transport via the depletion of their cross-links and the increase in mesh-size. This effect is especially important as the polymer crosses the reverse gelation point, due to the sudden disappearance of a mesh and the sudden increase in diffusivity.

Here, we study the interplay between hydrogel degradation and diffusion of ECM molecules and their deposition in an enzymatically degradable hydrogel through the solution of two coupled reaction-diffusion equations (2.1) and (3.3) of enzyme and ECM transport/deposition. We concentrate on the effect of three key parameters, re, κe and κm. Note that since the relative ECM size [r with tilde]m is typically much larger than unity in applications (i.e., larger than the mesh size of the hydrogel at any point during degradation), we keep it fixed in all simulations. Furthermore, for the sake of clarity, we consider a pair of parameters re and κe, which display, respectively, a wide and sharp degradation and for each case, investigate the effect of small and large rates of ECM deposition.

Results are shown in Fig. 6. subfigures (a) and (b) illustrate the case of a slow rate of ECM deposition coupled with a (a) wide and (b) sharp degradation front. In the first case, one can observe that the hydrogel undergoes significant degradation before ECM can be deposited. In the second case, free ECM can easily diffuse in the cavity left by the degradation front around a cell, but because of its slow conversion to solid ECM, growth remains concentrated in the peri-cellular region. By contrast, Fig. 6 (c) and (d) display the situation where the rate of ECM deposition is fast compared to the rate of degradation. When the degradation front is wide (Fig. 6 (c)), degradation is fairly uniform around cells, and the point of reverse gelation propagates slowly away from the cell boundary. Because of its relatively slow speed, this degradation front is immediately followed by a reservoir of free ECM molecules that quickly convert into a solid phase. Eventually, the construct is made up of a composite of deposited ECM surrounded by a region of non-degraded gel. In the case of a sharp degradation front, however, the region between intact and fully degraded gel is extremely thin, and the degradation profile may be thought of as an expanding sphere centered around a cell. Within this sphere, ECM can freely diffuse and deposit at a high rate. This yields a situation where the degradation front is immediately followed by a deposition front and a construct that resembles a hydrogel matrix filled with expanding inclusions of ECM. Interestingly, this case maximizes the growth to degradation ratio, that is, we observe a significant amount of growth for a minimal level of gel degradation. Although all of the above situations allow for growth in time, one may ask whether they could maintain a continuous mechanical integrity in time.

Fig. 6
Profile of polymer cross-link density, free ECM concentration and linked ECM concentration as a function of distance from the cell surface. Results are shown for three characteristic times during growth t=3.5, 7 and 10.5 day. (a) Slow ECM deposition in ...

5.3. Evolution of the construct's mechanical integrity during combined hydrogel degradation and ECM growth

The primary function of a tissue engineering scaffold is to provide a temporary mechanical support to cells as the new tissue grows. In vivo, this means that the construct must be able to resist physiological load at all time during its transition from scaffold to tissue. As discussed earlier, current hydrogel-based strategies that are used for cell encapsulations often lead to dissolution and failure of the construct before the ECM can bear any load. Here, we propose to numerically investigate whether the concept of localized degradation has the potential to switch this paradigm and allow both tissue growth and continuity of the construct's mechanical integrity. For this, we use techniques of numerical homogenization that allows us to estimate the overall stress-strain response of the unit cell (Fig. 4) in time as the spatial composition of the construct evolves. More precisely, an overall state of deformation [F with macron] is applied to the domain Ω and the corresponding displacement field u(x) and stress σ are calculated as (60):

u(x)=F¯x+u~(x),σ¯=(1V0¯Γ¯tXdΩ)F¯TJ¯
(5.1)

where ũ(x) are fluctuations in displacement that remain periodic on the domain's boundary, t is the traction force vector on Γ and J = det [F with macron]. The construct's Young's modulus may then be estimated by the secant modulus at a 5% strain for a uniaxial deformation:

E(t)=σ¯xx(t)ϵ¯xx=σ¯yy(t)ϵ¯yy=σ¯zz(t)ϵ¯zz
(5.2)

where ϵxx=ϵyy=ϵzz=5% and the last two equalities stem from the symmetry of the unit cell. To understand how this modulus correlates with the evolution of the hydrogel and ECM composition around a cell, we further introduce a measure of connectivity as follows. A solid phase, polymer or ECM, is defined as connected if any two points within that space can be connected by a continuous path. We note here that polymer and ECM are considered mechanically non-existent if ρρc and cm=0, respectively. Furthermore, since the unit cell is periodic, a continuous path can cross an external boundary and reappear on its opposite side.

To compare the stiffness of the gel to that of the linked ECM, we define a new dimensionless parameter E* = EECM/Egel where Egel=ρRT(1.05J013(2.1J0)1) is the secant modulus for 5% strain (31) and EECM = μ(3λ + 2μ)/(λ + μ). The evolution of the construct's properties is now affected by the three non-dimensional parameters re, κe, κm and the ECM properties μ* and λ*. However, to clarify the analysis, we show the four characteristic cases considered in Fig. 6 and for each, explore the effect of different ECM participation to the overall stiffness; specifically, we consider E* = 0.5 (the ECM is weaker than the gel), E* = 1 (the ECM stiffness is comparable to the gel) and E* = 1.5 (the ECM is stronger than the gel). Key results are presented in Fig. 7 for (a) a wide degradation front & slow deposition rate, (b) a sharp degradation front & slow deposition rate, (c) a wide degradation front & fast deposition rate and (d) a sharp degradation front & fast deposition rate. For each case, the figure shows the evolution of the Young's modulus (curve) and the connectivity of each phase (shaded blue for hydrogel and shaded red for ECM). Three dimensional contours of polymer cross-link density (blue) and linked ECM concentration (red) are also depicted at three characteristic times during the construct's evolution.

Fig. 7
Evolution of construct modulus over time. The mechanical integrity is preserved in the case of sharp degradation front and fast deposition relative to the front speed.

When deposition is slow compared to degradation (Fig. 7a and b), results show that the mechanical integrity of the construct monotonically drops with time until it completely dissolves (Ē = 0). A closer look at the ultrastructure evolution clearly shows that regardless of the sharpness and speed of the degradation front, ECM deposition lags behind and is unable to produce a well-connected phase before the hydrogel fully degrades. Interestingly, we note that a sharper degradation front yields a faster loss in construct's stiffness; a phenomenon that can be attributed to the fact that sharper fronts move relatively faster than wide ones (Fig. 5 b and d). When ECM deposits fast (relative to degradation) one predicts that it can reach connectivity before the hydrogel is completely degraded. As a consequence, even in the case of a wide degradation front, the model suggests that the construct does not completely lose its mechanical integrity. Indeed, while there exists a time interval for which none of the phases are connected, their mechanical interactions allow for a load transfer between them and an overall non-zero (although small) elastic modulus. In the situation of a sharp degradation front however, ECM can effectively grow within the empty interstices left by the propagating front. This eventually leads to a situation in which both polymer and ECM are connected and an optimized continuity of the overall construct's modulus.

6. DISCUSSION AND CONCLUDING REMARKS

In summary, we have constructed a multiphasic mixture model to represent the combined cell-mediated hydrogel degradation and tissue growth. The resulting model has the form of a coupled system of two reaction-diffusion equations corresponding to enzyme diffusion/hydrogel degradation and ECM diffusion/linkage. We have shown that according to the design of the hydrogel (in particular, its initial cross-link density and degradation kinetics), the system's behavior ranges from diffusion-dominated to diffusion-limited. The latter situation is associated with the appearance of a localized, spherical degradation front propagating outwards from each cell. This space creates pockets of unhindered space enabling the diffusion of large ECM molecules. If the rate of ECM linkage is significantly faster than the rate of hydrogel degradation (and therefore the degradation front velocity), the model suggests that spherical bodies of solid ECM may grow within cavities left by degradation and eventually connect in time. In this case, the construct structure displays a double connected network of ECM and hydrogel, which eventually allows a smooth transition between hydrogel and tissue and continuous mechanical integrity. In a nutshell, the model therefore points out that continuous mechanical integrity of the construct can be achieved by tuning the hydrogel design to achieve both a sharp and slow moving degradation front.

Although it does not appear explicitly, the model also captures the role of cell density during growth. Indeed, due to the non-dimensionalization procedure, the typical width w = w/L of the front is measured relative to the cell spacing L. Since L increases with decreasing cell density f as shown in equation (4.1), a degradation front appears sharper, and thus more favorable to the growth process for low cell densities. This also means that the predicted mechanism may be difficult to achieve for high cell density systems. Interestingly, experimental studies (61) have shown that strategies based on hydrolytic (bulk) degradation show a better potential when higher cell densities are employed. This suggests that optimized growth and degradation can be achieved with a controlled combination of bulk (hydrolytic) and localized (enzymatic) degradation kinetics. This dosage would move towards preferably hydrolytic in high cell density system to enzymatic in low cell-density systems. A quantitative analysis of these dynamics will be the object of future studies.

On a final note, it is important to mention that the scenario highlighted by our study may not be straightforward to reproduce experimentally due to the number of uncertainties and imperfections that characterizes both the polymer structure and the behavior of embedded cells. A number of effects may indeed affect the above model predictions. First of all, the homogeneous spatial distribution of cells is often not verified experimentally. Histology images (Fig. 3 b) indeed typically display a heterogeneous distribution containing local clusters of high cell density separated by low density regions. This distribution can further evolve in time via cell division, migration and death (62). This may affect the predicted mechanism by only allowing tissue growth in localized regions. To capture this mechanisms, the presented model would need to consider a large number of cells and their distribution, making the computational problem prohibitively costly. This issue can potentially be avoided by using multi-scale techniques based on homogenization as discussed in (63,64). Another uncertainty not captured by our model pertains to the behavior of cells, particularly in terms of enzyme and ECM production rate in time. Tissue-producing cells isolated from either tissue or derived from stem cells typically constitute a heterogeneous population that highly reactive to changes (chemical and mechanical) in their environment. The mechanisms discussed here could therefore display variations in space and time and be affected by the type of cell used. In addition to these, the enzymatic activity, i.e. degradation rate, can vary spatially due to the chemical and mechanical changes during the remodeling of the construct (32,65). Finally, cell proliferation is not explicitly included in the model. Although cell proliferation significantly influences the overall growth of neo-tissue, the increase in tissue growth on a per cell basis is less dramatic (66,67). Furthermore, the effect of cell proliferation on tissue synthesis will be captured within the rate of ECM synthesis per cell that used as an input parameter in the model. Therefore, the model, which describes ECM growth per cell, should provide a reasonable assessment of ECM growth. Further experimental characterizations are needed to address these questions.

To summarize, promoting tissue growth in a hydrogel is a complex problem that relies on a deep understanding of the physical, chemical and biological interaction between cells and hydrogels. In this complexity, an integrated approach merging modeling and experimental research will be crucial to identify the dominating mechanisms that control growth, degradation and continuous mechanical integrity of the tissue construct. This study is a step towards identifying these mechanisms. As such, we point to a direction in which experimental efforts should be oriented.

ACKNOWLEDGEMENT

Research reported in this publication was supported by the National Institute of Arthritis and Musculoskeletal and Skin Diseases of the National Institutes of Health under Award Number 1R01AR065441 to S.J.B and F.J.V. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. F.J.V. acknowledge the support of the National Science Foundation under the CAREER award 1350090. Work by UA was partially supported by the Department of Civil, Environmental & Architectural Engineering at the University of Colorado Boulder through a Doctoral Assistantship for Completion of Dissertation.

APPENDIX

A. METHODS: FINITE ELEMENT IMPLEMENTATION

In order to describe the tissue growth in porous scaffolds we use three balance law; balance of linear momentum (equation (3.6)), balance of mass (equation (3.3)) and the equation of osmotic pressure (equation (3.6)):

(σp+σmpI)=0
(A.1)
DcmDt=(Dmcm)cmu.κm(1c~mc~m0)cm
(A.2)
DceDt=(Dece)ceu.
(A.3)
lnϕw+(1ϕw)+χ(1ϕw)2+pvw=0
(A.4)

These governing equations are solved in three-dimension for the displacement field u, the pore pressure p, the concentration of enzyme ce and ECM cm. Moreover we treat the crosslinking density ρ and concentration of deposited ECM cm as internal state variables for which we use degradation/deposition kinetics law (equations (2.1) and (3.5)). We use finite element method to solve the system of equations for which the weak form is written as follows:

Vw1T(σp+σmpI)dV=S0w1Tt0dS
(A.5)
V(w2Tcmt+w2T(u.cm))dV+Vw2T(Dm)cmdV+Vw2Tκm(1c~mc~m)cmdV=S0w2TQmdS
(A.6)
V(w2Tcet+w2T(u.ce))dV+Vw2T(De)cedV=S0w2TQedS
(A.7)
Vw2T(lnϕw+(1ϕw)+χ(1ϕw)2+pvw)dV=0
(A.8)

Note that these equations are coupled and nonlinear in terms of the field variables. Also note that if the deformation field is small (5-10% strain) then the nonlinearities due to deformation and material behavior vanishes which leads the linear the balance of of linear momentum equation. The linearized forms are as follows:

Vw1T[(σpϵ+σmϵ)ϵuδu+Iδp]dV=S0w1Tt0dS
(A.9)
V(w2Tδcmt+w2T(δu.cm+u.δcm))dV+Vw2T(Dm)δcmdV+Vw2Tκm(1c~mc~m)δcmdV=S0w2TQmdS
(A.10)
V(w2Tδcet+w2T(δu.ce+u.δce))dV+Vw2T(De)δcedV=S0w2TQedS
(A.11)
Vw2T[(1ϕwϕw2χ(1ϕw))(1J2)Juδu+vwδp]dV=0
(A.12)

where the left and right hand sides are called internal Fint and external Fext force vectors respectively. After linearization our solution scheme becomes iterative in which we are solving for δd = [δu* δp* δce δcm] with the stopping criterion |R| ≈ 0 where R = FextFint. As the iterative method we chose Newton-Raphson, and in order to capture the transient behavior we implemented this method in backward Euler time integration scheme:

δ.d.i(t+Δt)=[Ci1(t+Δt)+Ki1(t+Δt)Δt]1Ri1(t+Δt)
(A.13)
d.i(t+Δt)=d.i1(t+Δt)+δd.i(t+Δt)
(A.14)
d(t+Δt)=d.(t+Δt)Δt+d(t)
(A.15)

where C, K, R, are the damping matrix stiffness matrix and residual vector respectively. In order to calculate these matrices and vector we use the FE discretization as follows. We used a mixed up formulation, that is 27-node element for u* and 8-node element for p* (68). Moreover for the concentration fields cm and ce 8-node element formulation is used. The test functions w1 and w2 and the interpolation of the fields then becomes:

w1=N27w1andw2=N8w2u=N27u,p=N8p,cm=N8cmandce=N8ce
(A.16)

Note that same interpolation rules are used for the the rates of these fields. The stiffness and damping matrices become:

K=[KuuKup00KpuKpp0000Kmm0000Kee]
(A.17)
C=[00000000Cmu0Cmm0Ceu00Cee]
(A.18)

where

Kuu=1#elV0B1T(Cp+Cm)B1dV
(A.19)
Kup=1#elV0B1T(I)N8dV
(A.20)
Kpu=1#elV0N8T(1ϕwϕw2χ(1ϕw))(1J2)adj(F)B1dV
(A.21)
Kpp=1#elV0N8T(vw)N8dV
(A.22)
Kmm=1#elV0[N8T(u.+κm(1c~mc~m))N8+B2T(Dm)B2]dV
(A.23)
Kee=1#elV0[N8T(u.)N8+B2T(De)B2]dV
(A.24)
Cmu=1#elV0[N8T(cm)N27]dV
(A.25)
Cmm=1#elV0[N8TN8]dV
(A.26)
Ceu=1#elV0[N8T(ce)N27]dV
(A.27)
Cee=1#elV0[N8TN8]dV
(A.28)

where

Cp=2ρH(ρρc)RTJ023J(C)
(A.29)
Cm=Cm~(λ(II)+μ(C)
(A.30)
B1=[N27x000N27x000N27x]
(A.31)
B2=[N8x]
(A.32)

Verification of the model is performed by comparing the solution with that of 1D degradation-diffusion model presented in (33) (Figure 8c). In order to check the consistency of the model, we investigated the spatial and temporal change of crosslinking density ρ along the diffusion direction. Note from equation (2.1) that any numerical error in the enzyme transport affects the calculated crosslinking density. The error in crosslinking density depending on the mesh size (Figure 8 a) and time step (Figure 8 b) is calculated using L2 norm:

e2ρ=[V(ρk+1ρk)2dV]12
(A.33)

The convergence rate of the error is determined by fitting a least square curve to the calculated errors. The convergence rate with respect to the mesh size is quadratic, in agreement with theoretical predictions (68). Moreover the convergence of the Euler time integration scheme is linear as expected (69). For the analysis in this paper we chose a mesh size and time step given by h* = 0.08 (shown in red circle in Figure 8 a) and Δt* = 0.006 (shown in red circle in Figure 8 b), which provides reasonable accuracy (Figure 8c).

Fig. 8
Convergence study. (a) Convergence of the error measure in (A.34) as a function of mesh size h. (b) Convergence of the error measure in (A.34) as a function of time step. (c) Comparison of the model prediction with the one presented in (33) for the evolution ...

REFERENCES

1. Yannas I, Lee E, Orgil DP, Skrabut EM, Murphy GF. Synthesis and characterization of a model extracellular matrix that induces partial regeneration of adult mammalian skin. Proc Natl Acad Sci. 1989;86(3):933–7. [PubMed]
2. Howard D, Buttery LD, Shakesheff KM, Roberts SJ. Tissue engineering: strategies, stem cells and scaffolds. J Anat. 2008;213(1):66–72. [PubMed]
3. Ikada Y. Challenges in tissue engineering. J R Soc Interface. 2006;3:589–601. [PMC free article] [PubMed]
4. Drury JL, Mooney DJ. Hydrogels for tissue engineering: scaffold design variables and applications. Biomaterials. 2003;24(24):4337–51. [PubMed]
5. Armstrong CG, Mow VC. Variations in the intrinsic mechanical properties of human articular cartilage with age, degeneration, and water content. J Bone Jt Surg. 1982;64(1):88–94. [PubMed]
6. Villanueva I, Hauschulz DS, Mejic D, Bryant SJ. Static and dynamic compressive strains influence nitric oxide production and chondrocyte bioactivity when encapsulated in PEG hydrogels of different crosslinking densities. Osteoarthritis Cartilage. 16(8):909–18. [PMC free article] [PubMed]
7. Nicodemus GD, Bryant SJ. Cell encapsulation in biodegradable hydrogels for tissue engineering applications. Tissue Eng Part B Rev. 2008;14(2):149–65. [PMC free article] [PubMed]
8. Ratner BD, Bryant SJ. BIOMATERIALS: Where We Have Been and Where We are Going. Annu Rev Biomed Eng. 2004;6:41–75. [PubMed]
9. Nicodemus GD, Bryant SJ. The role of hydrogel structure and dynamic loading on chondrocyte gene expression and matrix formation. J Biomech. 2008;41(7):1528–36. [PMC free article] [PubMed]
10. Bryant SJ, Anseth KS. Controlling the spatial distribution of ECM components in degradable PEG hydrogels for tissue engineering cartilage. J Biomed Mater Res A. 2003;64(1):70–9. [PubMed]
11. Bryant SJ, Anseth KS. Hydrogel properties influence ECM production by chondrocytes photoencapsulated in poly(ethylene glycol) hydrogels. J Biomed Mater Res. 2002;59(1):63–72. [PubMed]
12. Chung C, Mesa J, Randolph MA, Yaremchuk M, Burdick JA. Influence of Gel Properties on Neocartilage Formation by Auricular Chondrocytes Photoencapsulated in Hyaluronic Acid Networks. J Biomed Mater Res A. 2006;77(3):518–25. [PMC free article] [PubMed]
13. Park H, Guo X, Temenoff JS, Tabata Y, Caplan AI, Kasper FK, et al. Effect of swelling ratio of injectable hydrogel composites on chondrogenic differentiation of encapsulated rabbit marrow mesenchymal stem cells in vitro. Biomacromolecules. 2009;10(3):541–6. [PMC free article] [PubMed]
14. Hanjaya-Putra D, Wong KT, Hirotsu K, Khetan S, Burdick JA, Gerecht S. Spatial Control of Cell-Mediated Degradation to Regulate Vasculogenesis and Angiogenesis in Hyaluronan Hydrogels. Biomaterials. 33(26):6123–31. [PMC free article] [PubMed]
15. Levesque SG, Shoichet MS. Synthesis of Enzyme-Degradable, Peptide-Cross-Linked Dextran Hydrogels. Bioconjug Chem. 2007;18(3):874–85. [PubMed]
16. Skalak R, Dasgupta G, Moss M. Analytical Description of Growth. J Theor Biol. 1982;94:555–577. [PubMed]
17. Rodriguez EK, Hoger A, Mcculloch AD. Stress-dependent finite growth in soft elastic tissues. J Biomech. 1994;27(4):455–67. [PubMed]
18. Klisch SM, Chen SS, Sah RL, Hoger A. A growth mixture theory for cartilage with application to growth-related experiments on cartilage explants. J Biomech Eng. 2003;125(2):169–79. [PubMed]
19. Klisch SM, Van Dyke TJ, Hoger A. A theory of volumetric growth for compressible elastic biological materials. Math Mech Solids. 2001;6(6):551–75.
20. Humphrey JD, Rajagopal KR. A Constrained Mixture Model for Growth and Remodeling of Soft Tissues. Math Models Methods Apllied Sci. 2002;12(3):407–30.
21. Ateshian GA, Costa KD, Azeloglu EU, Morrison B, Hung CT. Continuum modeling of biological tissue growth by cell division and alteration of intracellular osmolytes and extracellular fixed charge density. J Biomech Eng. 2009;131:101001. [PMC free article] [PubMed]
22. Ateshian GA, Ricken T. Multigenerational interstitial growth of biological tissues. Biomech Model Mechanobiol. 2010;9(6):689–702. [PMC free article] [PubMed]
23. Garikipati K, Arruda EM, Grosh K, Narayanan H, Calve S. A continuum treatment of the growth in biological tissue: the coupling of mass transport and mechanics. J Mech Phys Solids. 2004;52:1595–625.
24. Vernerey FJ. A mixture approach to investigate interstitial growth in engineering scaffolds. Biomech Model Mechanobiol. 2015:1–20. [PMC free article] [PubMed]
25. Dhote V, Skaalure SC, Akalp U, Roberts J, Bryant SJ, Vernerey FJ. On the role of hydrogel structure and degradation in controlling the transport of cell-secreted matrix molecules for engineered cartilage. Jounarl Mech Behav Biomadeical Mater. 2013;19:61–74. [PMC free article] [PubMed]
26. Dhote V, Vernerey FJ. Mathematical model of the role of degradation on matrix development in hydrogel scaffold. Biomech Model Mechanobiol. 2014;13:167–83. [PMC free article] [PubMed]
27. Trewenack A, Please C, Landman K. A continuum model for the development of tissue-engineered cartilage around a chondrocyte. Math Med Biol. 26(3):241–262. [PubMed]
28. Haider MA, Olander JE, Arnold RF, Marous DR, McLamb AJ, Thompson KC, et al. A phenomenological mixture model for biosynthesis and linking of cartilage extracellular matrix in scaffolds seeded with chondrocytes. Biomech Model Mechanobiol. 2011;10:915–24. [PMC free article] [PubMed]
29. Sengers B, Taylor M, Please CP, Oreffo ROC. Computational modelling of cell spreading and tissue regeneration in porous scaffolds. Biomaterials. 2007;28:1926–40. [PubMed]
30. Wilson CG, Bonassar LJ, Kohles SS. Modeling the dynamic composition of engineered cartilage. Arch Biochem Biophys. 2002;408(2):246–54. [PubMed]
31. Akalp U, Chu S, Skaalure S, Bryant SJ, Doostan A, Vernerey FJ. Determination of the polymer-solvent interaction parameter for PEG hydrogels in water: Application of a self learning algorithm. Polymer. 2015 [PMC free article] [PubMed]
32. Lutolf MP, Lauer-Fields JL, Schmoekel HG, Metters AT, Weber FE, Fields GB, et al. Synthetic matrix metalloproteinase-sensitive hydrogels for the conduction of tissue regeneration: Engineering cell-invasion characteristics. Proc Natl Acad Sci. 2003;100(9):5413–8. [PubMed]
33. Skaalure SC, Akalp U, Vernerey FJ, Bryant SJ. Tuning Reaction and Diffusion Mediated Degradation of Enzyme-Sensitive Hydrogels. Adv Healthc Mater. 2016;5(4):432–8. [PMC free article] [PubMed]
34. Roberts JJ, Nicodemus GD, Greenwald EC, Bryant SJ. Degradation improves tissue formation in (un)loaded chondrocyte-laden hydrogels. Clin Orthop. 2011;469(10):2725–34. [PMC free article] [PubMed]
35. Lustig SR, Peppas NA. Solute Diffusion in Swollen Membranes. IX Scaling Laws for Solute Diffusion in Gels. J Appl Polym Sci. 1988;36(4):735–47.
36. Skaalure SC, Chu S, Bryant SJ. An Enzyme-Sensitive PEG Hydrogel Based on Aggrecan Catabolism for Cartilage Tissue Engineering. Adv Healthc Mater. 2015 Feb 1;4(3):420–31. [PMC free article] [PubMed]
37. Glasstone S, L K. J., Eyring H. The Theory of Rate Processes The Kinetics of Chemical Reactions, Viscosity, Diffusion and Electrochemical Phenomena. Mcgraw-Hill Book Compagny; New York, N.Y.: 1941.
38. Cohen MH, Turnbull D. Molecular Transport in Liquids and Glasses. J Chem Phys. 1959 Nov 1;31(5):1164–9.
39. Peppas NA, Reinhart CT. Solute diffusion in swollen membranes. Part I. A new theory. J Membr Sci. 1983 Oct 1;15(3):275–87.
40. Miller CC. The Stokes-Einstein Law for Diffusion in Solution. Proc R Soc Lond Math Phys Eng Sci. 1924 Dec 1;106(740):724–49.
41. Nicodemus GD, Skaalure SC, Bryant SJ. Gel structure impacts pericellular and extracellular matrix deposition which subsequently alters metabolic activities in chondrocyte-laden PEG hydrogels. Acta Biomater. 2011;2:492–504. [PMC free article] [PubMed]
42. Bowen RM. Incompressible porous media models by use of the theory of mixtures. Int J Eng Sci. 1980;18(9):1129–48.
43. Biot MA. Mechanics of Deformation and Acoustic Propagation in Porous Media. J Appl Phys. 1962;33(4):1482–98.
44. Li C, Borja RI, Regueiro RA. Dynamics of porous media at finite strain. Comput Methods Appl Mech Eng. 2004;193:3837–70.
45. Brinkman HC. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Appl Sci Res. 1949;1(1):27–34.
46. Whitaker S. Flow in porous media I: A theoretical derivation of Darcy's law. Transp Porous Media. 1986 Mar;1(1):3–25.
47. Holmes MH, Mow VC. The nonlinear characteristics of soft gels and hydrated connective tissues in ultrafiltration. J Biomech. 1990 Jan 1;23(11):1145–56. [PubMed]
48. Vernerey FJ. A theoretical treatment on the mechanics of interfaces in deformable porous media. Int J Solids Struct. 2011 Nov;48(22–23):3129–41.
49. Vernerey FJ. The Effective Permeability of Cracks and Interfaces in Porous Media. Transp Porous Media. 2012 Mar 22;93(3):815–29.
50. Cannon WB. ORGANIZATION FOR PHYSIOLOGICAL HOMEOSTASIS. Physiol Rev. 1929;9(3):399–431.
51. Flory PJ. Principles of polymer chemistry. Cornell University Press; 1953.
52. Ateshian GA. On the theory of reactive mixtures for modeling biological growth. Biomech Model Mechanobiol. 2007;6(6):423–45. [PMC free article] [PubMed]
53. Mow VC, Gibbs MC, Lai WM, Zhu WB, Athanasiou KA. Biphasic Indentation of Articular Cartilage - II A Numerical Algorithm and An Experimental Study. J Biomech. 1989;22(8/9):853–61. [PubMed]
54. Wong M, Ponticiello M, Kovanen V, Jurvelin JS. Volumetric changes of articular cartilage during stress relaxation in unconfined compression. J Biomech. 2000;33:1049–54. [PubMed]
55. Vernerey FJ, Liu WK, Moran B, Olson GB. A Micromorphic Model for the Multiple Scale Failure of Heterogeneous Materials. J Mechancis Phys Solids. 2008;56(4):1320–47.
56. Tyrus JM, Gosz M, DeSantiago E. A local finite element implementation for imposing periodic boundary conditions on composite micromechanical models. Int J Solids Struct. 2007;44:2972–89.
57. Miller DR, Macosko CW. A New Derivation of Post Gel Properties of Network Polymers. Macromolecules. 1976;9(2):206–11.
58. Wilhelm SM, Eisen AZ, Teter M, Clark SD, Kronberger A, Goldberg G. Human fibroblast collagenase: Glycosylation and tissue-specific levels of enzyme synthesis. Proc Natl Acad Sci. 1986;83:3756–60. [PubMed]
59. Gribbon P, Hardingham TE. Macromolecular Diffusion of Biological Polymers Measured by Confocal Fluorescence Recovery after Photobleaching. Biophys J. 75:1032–9. [PubMed]
60. Nemat-Nasser S, Hori M. Applied Mathematics and Mechanics. Vol. 37. Elsevier Science Publishers; North-Holland: 1993. Micromechanics: Overall Properties of Heterogeneous Materials.
61. Rice MA, Anseth KS. Controlling Cartilaginous Matrix Evolution in Hydrogels with Degradation Triggered by Exogenous Addition of an Enzyme. Tissue Eng. 2007;13(4):683–91. [PubMed]
62. Roberts JJ, Bryant SJ. Comparison of photopolymerizable thiol-ene PEG and acrylate-based PEG hydrogels for cartilage development. Biomaterials. 2013;34(38):9969–79. [PMC free article] [PubMed]
63. Kouznetsova VG, Geers MGD, Brekelmans WAM. Multiscale second-order computational homogenization of multi-phase materials: a nested finite element strategy. Comput Methods Appl Mech Eng. 2004;193:5525–50.
64. Vernerey FJ. A Microstructure-Based Continuum Model for Multiphase Solids. Mech Adv Mater Struct. 2014;21(6):441–56.
65. Nicodemus GD, Bryant SJ. Mechanical loading regimes affect the anabolic and catabolic activities by chondrocytes encapsulated in PEG hydrogels. Osteoarthritis Cartilage. 2010;18(1):126–37. [PubMed]
66. Huang AH, Yeger-McKeever M, Stein A, Mauck RL. Tensile properties of engineered cartilage formed from chondrocyte- and MSC-laden hydrogels. Osteoarthritis Cartilage. 2008 Sep;16(9):1074–82. [PMC free article] [PubMed]
67. Yodmuang S, McNamara SL, Nover AB, Mandal BB, Agarwal M, Kelly T-AN, et al. Silk microfiber-reinforced silk hydrogel composites for functional cartilage tissue repair. Acta Biomater. 2015 Jan 1;11:27–36. [PMC free article] [PubMed]
68. de Borst R, Crisfield MA, Remmers JJC, Verhoosel CV. Computational Mechanics. 2nd ed. Wiley; United Kingdom: 2012. Non-linear Finite Element Analysis of Solid and Structures.
69. Mathews JH, Fink KD. Numerical Methods Using Matlab. 3rd ed. Prentice Hall; New Jersey: 1999.