The interplay between tissue growth and scaffold degradation in engineered tissue constructs

In vitro tissue engineering is emerging as a potential tool to meet the high demand for replacement tissue, caused by the increased incidence of tissue degeneration and damage. A key challenge in this field is ensuring that the mechanical properties of the engineered tissue are appropriate for the in vivo environment. Achieving this goal will require detailed understanding of the interplay between cell proliferation, extracellular matrix (ECM) deposition and scaffold degradation. In this paper, we use a mathematical model (based upon a multiphase continuum framework) to investigate the interplay between tissue growth and scaffold degradation during tissue construct evolution in vitro. Our model accommodates a cell population and culture medium, modelled as viscous fluids, together with a porous scaffold and ECM deposited by the cells, represented as rigid porous materials. We focus on tissue growth within a perfusion bioreactor system, and investigate how the predicted tissue composition is altered under the influence of (1) differential interactions between cells and the supporting scaffold and their associated ECM, (2) scaffold degradation, and (3) mechanotransduction-regulated cell proliferation and ECM deposition. Numerical simulation of the model equations reveals that scaffold heterogeneity typical of that obtained from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}CT scans of tissue engineering scaffolds can lead to significant variation in the flow-induced mechanical stimuli experienced by cells seeded in the scaffold. This leads to strong heterogeneity in the deposition of ECM. Furthermore, preferential adherence of cells to the ECM in favour of the artificial scaffold appears to have no significant influence on the eventual construct composition; adherence of cells to these supporting structures does, however, lead to cell and ECM distributions which mimic and exaggerate the heterogeneity of the underlying scaffold. Such phenomena have important ramifications for the mechanical integrity of engineered tissue constructs and their suitability for implantation in vivo.

deposited by the cells, represented as rigid porous materials. We focus on tissue growth within a perfusion bioreactor system, and investigate how the predicted tissue composition is altered under the influence of (1) differential interactions between cells and the supporting scaffold and their associated ECM, (2) scaffold degradation, and (3) mechanotransduction-regulated cell proliferation and ECM deposition. Numerical simulation of the model equations reveals that scaffold heterogeneity typical of that obtained from μCT scans of tissue engineering scaffolds can lead to significant variation in the flow-induced mechanical stimuli experienced by cells seeded in the scaffold. This leads to strong heterogeneity in the deposition of ECM. Furthermore, preferential adherence of cells to the ECM in favour of the artificial scaffold appears to have no significant influence on the eventual construct composition; adherence of cells to these supporting structures does, however, lead to cell and ECM distributions which mimic and exaggerate the heterogeneity of the underlying scaffold. Such phenomena have important ramifications for the mechanical integrity of engineered tissue constructs and their suitability for implantation in vivo.

Introduction
Mathematical modelling of tissue growth is a wide field of research, aiming to provide a more complete understanding of the myriad biological and biophysical processes that contribute to tissue growth. Such theoretical models underpin the emerging field of in vitro tissue engineering, which, by the creation of replacement tissue in the laboratory, has the potential to alleviate the shortage of replacement tissue available for implantation into patients. A typical method for generating such implants entails seeding a biodegradable porous scaffold with cells; subsequent incubation in a bioreactor allows the cells to colonise the porous scaffold (termed a tissue construct). On implantation, the degrading scaffold is replaced by extracellular materials such as collagen and proteoglycans, which are laid down by the cells (Freed et al. 1994). Ensuring that the rates of nascent tissue growth and scaffold degradation (e.g., due to hydrolysis) are appropriately matched is therefore crucial in maintaining the mechanical integrity of the construct, a factor of especial importance for load-bearing constructs, such as bone implants (Wu and Ding 2004). The biological processes which contribute to tissue construct growth operate on disparate spatio-temporal scales and range from intracellular gene networks to tissue-level mechanics; reviews are given by Curtis and Riehle (2001), Cowin (2000), Cowin (2004), Sipe (2002) and Burdick and Mauck (2010). In this paper, we concentrate on a tissue-scale description of tissue growth, and employ a continuum model to focus on the way in which the properties of the supporting scaffold influence the structure of the resulting tissue construct.
In addition to the scaffold's mechanical properties, its chemical features are of great importance. For example, most cell types are anchorage-dependent, their growth being affected by interactions between a substrate or deposited extracellular matrix (ECM); the surface chemistry of the scaffold crucially affects such interactions. Adherence to polymer scaffolds commonly employed in tissue engineering applications is mediated by adsorption of deposited ECM molecules onto the scaffold surface or by, for instance, artificially embedded cell recognition molecules (Freed and Vunjak-Novakovic 1998;Nikolovski and Mooney 2000). The type and density of such molecules may vary dramatically throughout the scaffold due to, e.g., inhomogeneous ECM deposition, leading to spatial variations in cell adhesion characteristics, or preferential adherence to ECM or other deposited materials over artificial scaffolds.
Controlling the biochemical environment of cells (both with respect to nutrient/oxygen delivery and the provision of growth factors and other cell-signalling molecules) is key to producing constructs of a size appropriate for implant, while minimising the necrotic core which often forms at the centre of tissue constructs in static culture. To this end perfusion bioreactors are frequently employed, which exploit advection of culture medium to enhance mass transport. This approach also provides mechanical stimulation to cells contained within a porous scaffold via culture medium flow-induced shear stress, and can allow for the addition of cyclical compressive loads. For example, El Haj et al. (1990) have developed a perfusion/compression bioreactor system which comprises a poly(l-lactic acid) (PLLA) scaffold, through which culture medium is perfused via a peristaltic pump; macroscale compression of the scaffold may also be effected by the addition of a piston. (The layout of the bioreactor system in the absence of macroscale compression is depicted in Fig. 1.) Cells contained within a porous scaffold are therefore subjected to culture medium flow-induced shear stress and macroscale strain, in addition to mechanical interactions which exist between adjacent cells and between cells and scaffold/ECM. Other strategies for applying mechanical stimulation to tissue engineered constructs are reviewed by Martin et al. (2004) and Cartmell and El Haj (2005). The process by which such stimuli are integrated into the cellular response, for example, in terms of its proliferative behaviour, is known as mechanotransduction. The biochemical and biomechanical environment required for optimum growth is specific to the tissue under consideration; bespoke bioreactors are therefore required to provide appropriate cues for different tissue engineering applications. Well-studied examples include osteocytes (terminallydifferentiated bone cells), which are known to be sensitive to fluid shear stress (Bakker et al. 2004a); and chondrocytes, whose metabolism and maintenance of ECM integrity are regulated by mechanical stress (Urban 1994;Wang et al. 2010). Indeed, in El Haj et al. (1990 attention focussed on the influence of such mechanical stimulation on bone tissue growth. In what follows, we employ a mathematical model relevant to perfusion bioreactor systems in which cells are cultivated within a porous scaffold. Our formulation accommodates the cells' progression from a proliferative to an apoptotic phenotype, via an ECM depositing phase, in response to changes in the local cell volume fraction (a methodology for accomodating mechanotransduction-mediated cell proliferation in response to a range of mechanical stimuli is given in O' Dea et al. (2008Dea et al. ( , 2010 and ), as well as considering in detail the interactions between the cells and their supporting structures.
A variety of approaches has been employed to model tissue growth, their respective benefits depending on the specific application under consideration. Here, we study an Cell seeded scaffold Culture medium reservoir Peristaltic pump Fig. 1 Layout of the bioreactor system of El Haj et al. (1990) extension to a recently-developed continuum model  in which we consider the evolution of the spatial distribution of PLLA scaffold and ECM density. We employ a multiphase formulation which enables us to incorporate interactions between the many constituent materials which comprise biological tissue; we model explicitly cell-cell and cell-scaffold/ECM interactions as well as mass transfer between phases (representing cell proliferation, ECM deposition and scaffold degradation). Such multiphase approaches have been widely employed in industrial applied mathematics (Drew and Segel 1971) and, more recently, modelling of tumour growth and in vitro tissue engineering processes; examples include Breward et al. (2002), Byrne and Preziosi (2003), Franks and King (2003), Araujo and McElwain (2005), Lemon et al. (2006), Lemon and King (2007), Wilson et al. (2007), O'Dea et al. (2008, , and references therein. Reviews are given by, e.g., Tosin (2009) andO'Dea et al. (2012).
In O'Dea et al. (2010), tissue growth within a perfusion bioreactor was modelled, using a three-phase continuum model in a 2D channel geometry. In common with multiphase models of similar biological systems (Landman and Please 2001;Byrne and Preziosi 2003;Franks and King 2003;Lemon et al. 2006), the cells and associated ECM were represented as a viscous fluid phase that is distinct from the culture medium; and the porous scaffold was modelled as a rigid porous medium. Two factors of key importance to modelling the growth and adaptation of engineered tissue constructs were investigated: (1) cell-cell and cell-scaffold interactions and, (2) mechanotransduction mechanisms. The formulation was simplified via the long-wavelength limit (in which the bioreactor's aspect ratio is assumed to be small) and by considering constant, spatially-homogeneous scaffold porosity. Numerical simulation of the model equations (validated by analytic solutions obtained in the limit of asymptotically-small cell volume fraction), revealed that inclusion of cell-cell and cell-scaffold interactions leads to significant differences in the extent to which the cell population colonises the scaffold, depending upon the relative importance of cell aggregation and repulsion. It was further shown that the composition of the resulting construct was strongly influenced by whether cell proliferation and ECM deposition were regulated by mechanical stimulation related to the cell population density, pressure or shear stress. Employing two-dimensional finite element simulations,  demonstrated that, when considering total tissue yield, the long-wavelength limit of O'Dea et al. (2010) provides an excellent approximation to the full two-dimensional model, even for relatively large values of bioreactor aspect ratio. However, this work further demonstrated that mechanotransduction-mediated tissue growth can lead to significant twodimensional spatial variation of tissue density, a feature which is not captured by the long-wavelength limit. The authors concluded that, while spatial effects in two-or three-dimensions cannot be ignored in comprehensive models of tissue growth, its relative simplicity makes the long-wavelength model a natural framework with which to estimate parameters relevant to specific bioreactor systems, for subsequent use in more complex two-or three-dimensional models.
Whilst the inclusion of the scaffold phase is a significant departure from two-fluid models [see, for example, Franks and King (2003)] since interactions between the rigid scaffold and cell phases are explicitly modelled, a key assumption of O' Dea et al. (2010) and  is that the scaffold phase is spatially-homogeneous and constant in time, leading to significant simplification of the model formulation. Lemon and King (2007) considered non-uniform porosity only near the edge of the scaffold; however, μCT scans of porous scaffolds typically employed in tissue engineering applications indicate significant heterogeneity throughout the scaffold, as well as inhomogeneous deposition of extracellular materials (see, e.g., Yang and El Haj (2006) for a discussion of tissue engineering scaffolds). Figure 2a shows experimental data indicating the cross section-averaged scaffold volume fraction along the axial length of a typical scaffold of the type employed in El Haj et al. (1990), together with that following culture, indicating the level of mineralisation by osteocytes in such a scaffold. Such heterogeneity in scaffold density is likely to have implications for the mechanical properties of the resulting tissue construct, these considerations being especially important where the tissue is to be load-bearing, as is the case for cartilage or bone tissue: areas of weakness in implanted tissues may fail under physiological loading. For this reason, we aim to use our model to study how experimentally-relevant spatial variations in scaffold porosity may influence construct composition, and to indicate its importance in the generation of viable replacement tissue.
A simple study considering scaffold degradation and ECM deposition is given by Haider et al. (2010), in which the (spatially-independent) evolution of the total scaffold and ECM density is calculated using a phenomenological mixture model; specific consideration is given to scaffold-ECM linkage. The study concludes that the initial scaffold density will affect significantly the resulting construct's material properties. More complex models, considering spatial dependence, include Kelly and Prendergast (2003), in which a core of underdeveloped tissue in a poroelastic model of cartilage tissue was considered, showing that construct inhomogeneity dramatically reduces its mechanical integrity. Optimal design of porous scaffolds was discussed in Adachi et al. (2006) where the interplay between tissue growth and scaffold degradation, as well as the scaffold microstructure, was considered. Sanz-Herrera et al. (2008) developed a multiscale model to investigate the interplay between scaffold design parameters and bone tissue regeneration. By constructing solutions using finite element methods, the authors indicated that bone regeneration increases with scaffold stiffness and mean pore size. Byrne et al. (2007) considered the influence of scaffold porosity and degradation rate on in vitro bone tissue growth; mechanotransduction-regulated differentia- tion of stem cells to fibroblasts, chondrocytes and osteoblasts was also accommodated, each phenotype displaying different migration and material properties. A random walk model for cell movement within a poroelastic scaffold (whose deformation was simulated via a finite element method) was used to compute the tissue composition. The study concluded that under low load, high porosity and stiffness, together with an intermediate scaffold degradation rate, stimulate increased bone tissue generation; to prevent collapse of the scaffold under high load, a reduced degradation rate is required. Due to the nature of the problems investigated (three-dimensional scaffolds with specific pore geometry, deformation, flow, scaffold material properties), many of the studies mentioned above give rise to complex systems of coupled PDEs, which are heavily reliant on numerical investigation. In this study, we demonstrate that considerations relevant to biological tissue growth may be accommodated within a continuum model, amenable to asymptotic simplification; we extend our earlier work ) by relaxing the assumption of constant, homogeneous scaffold porosity, and employ the resulting model to investigate the interplay between scaffold degradation and nascent tissue growth in a perfusion bioreactor. Our model incorporates a cell population and culture medium, each represented as a viscous fluid, as well as both a PLLA scaffold and ECM deposited by the cells, modelled as rigid porous phases. This approach allows spatio-temporal variations in cell-scaffold and cell-ECM interactions, realistic scaffold porosity distributions (informed by experimental data), scaffold degradation, and mechanotransduction-regulated cell proliferation and ECM deposition to be accommodated.
In the limit of asymptotically-small bioreactor aspect ratio, the resulting model simplifies to three nonlinear differential equations. Via this formulation, we seek to provide a more comprehensive description of in vitro tissue growth, while remaining within a simplified modelling framework. Our investigations reveal that spatial inhomogeneity in scaffold volume fraction strongly influences the cells' mechanical environment, leading to inhomogeneous cell proliferation and ECM deposition. Further, our model suggests that preferential adherence to ECM in favour of the PLLA scaffold has no significant influence on the eventual construct composition; we therefore conclude that such additional mathematical complexity is unnecessary, so that simplified models, in which cells interact uniformly with their supporting structures, may be employed to describe biological tissue growth. Lastly, we indicate that careful manipulation of the rate of PLLA scaffold degradation is required in order to maintain the mechanical integrity of constructs.
The remainder of the paper is organised as follows. In §2, the multiphase model of O'Dea et al. (2010) is summarised, and extended by the addition of spatial and temporal variation in scaffold and ECM volume fractions, and the resulting governing equations and boundary conditions are stated (a detailed derivation is provided in the Appendix). In §3, numerical simulations of the model equations are presented and the importance of scaffold degradation, ECM deposition and heterogeneity in scaffold porosity on construct composition is investigated. §4 provides a summary of the results contained in the preceding sections and a discussion of their implications for in vitro tissue engineering, together with suggestions for future avenues of investigation.

Model formulation
In this section, we present a multiphase model which describes the growth of a tissue construct within a nutrient-rich perfusion bioreactor. The bioreactor under consideration is illustrated in Fig. 1 and comprises a cell-seeded porous scaffold within a culture medium-filled cylinder, through which a flow is driven (see El Haj et al. (1990)  For simplicity, we view the perfusion bioreactor as a two-dimensional channel (we expect results for an axisymmetric cylinder to be qualitatively similar) containing a mixture of four interacting phases. The cell population and culture medium are modelled as distinct viscous fluids, and rigid porous phases represent the PLLA scaffold and the ECM. The interplay between cell proliferation, ECM deposition and scaffold degradation is captured by mass exchange between the relevant phases, effected by the specification of mass transfer functions which account for the influence of mechanotransduction on cell proliferation and ECM deposition.
The mechanical interactions between phases comprise interphase viscous drag (proportional to differences in phase velocity) and active forces. The latter enter the governing equations via prescribed contributions to the cell phase pressure, arising due to cell-cell interactions and traction between the cell and scaffold or ECM phases, respectively. Interactions between the culture medium and scaffold/ECM phases are assumed to involve only viscous drag.
The mechanical interactions in this four phase formulation are simplified by lumping the rigid PLLA scaffold and ECM together into a single phase, referred to as the 'substrate' [a similar approach is used by Lubkin and Jackson (2002)]. We neverthe-less track individually the evolution of the PLLA scaffold and ECM volume fractions, which allows us to distinguish between cell-scaffold and cell-ECM interactions. From a mechanical point of view, this model may therefore be thought of as a three phase system, differing from O' Dea et al. (2010) in the sense that the scaffold and ECM phases are mechanically identical, but chemically distinct. We further assume that the bioreactor has a small aspect ratio, so that significant simplification of the governing equations can be achieved. We remark that the dimensions of the bioreactor system of El Haj et al. (1990) are inconsistent with this simplifying limit; however, our previous work ) indicates that such a limit provides an excellent approximation to the full two-dimensional model. We employ this limit in preference to an a priori assumption of one-dimensional flow for consistency with our previous work, and for wider applicability to other bioreactor systems.

Dimensionless model equations and boundary conditions
We consider a 2D Cartesian coordinate system x = (x, y), in which the bioreactor is assumed to occupy the dimensionless region 0 x 1, 0 y h 1, and within which the PLLA scaffold phase is localised in the region a We assume that all dependent variables are functions of x and dimensionless time, t.
The volume fractions of the cell, culture medium, PLLA scaffold and ECM phases are denoted by θ n , θ w , θ s and θ e , respectively; and the substrate phase is defined by (1) The dimensionless velocities and pressures of the cell and culture medium phases are denoted u i = (u i , v i ) and p i (i = n, w). The rigidity of the PLLA scaffold and ECM implies u s = u e = 0; the solid phase pressures p s , p e are not required in this analysis and remain undetermined. Tissue growth, scaffold degradation and ECM deposition are captured via material transfer functions S i which we specify below. The governing equations are stated below in dimensionless form.
The model is constructed by considering mass and momentum balances for each phase, assuming that the fluid phases are incompressible with equal density, and by neglecting inertial effects (details of the model derivation and nondimensionalisation are provided in the Appendix). Assuming that the bioreactor aspect ratio is asymptotically small, and employing the momentum balance equations, together with the no-voids condition, i θ i = 1, it is straightforward to show that the flow is unidirectional and that the pressure and volume fraction of each phase are functions of x and t. By eliminating dependent variables, the system may be reduced to the following differential equations for θ s , θ e , θ n and p w : in which μ n is the relative viscosity of the cell and culture medium phases, and and ψ capture active forces that exist between adjacent cells, and between cells and the substrate. We note that appropriate choice of the interaction functions and ψ and mass transfer terms S i ensures that the volume fractions obey θ i ∈ [0, 1] (see e.g. Preziosi and Tosin (2009) for more details). Equations (2)-(4) are the mass conservation equations for the PLLA scaffold, ECM and cell phases. The latter states that the rate of change of the cell volume fraction θ n is due to advection and the proliferation of cells; PLLA scaffold and ECM evolution is due to deposition of ECM and scaffold degradation only, and therefore θ e , θ s have an implicit spatial dependence. Equation (5) embodies conservation of mass for the multiphase mixture. In view of Eq. (1), Eqs.
(2)-(5) may be recast as a set of three nonlinear differential equations for , θ n and p w , with (2) and (3) replaced with reflecting the pseudo-three phase nature of the model. Equations (4) and (5) are closed by imposing the following boundary conditions: where P U and P D are the imposed up-and downstream pressures; initial conditions for θ n , θ s and θ e are specified subsequently.
In this paper, we consider the influence of scaffold properties on tissue construct growth and composition by examining the interplay between scaffold degradation and nascent tissue growth under the influence of biologically-relevant cell-substrate interactions, and mechanotransduction-regulated cell proliferation and ECM deposition. These phenomena are captured by the interaction functions , ψ and the material transfer functions S n , S e and S s which are specified in terms of the dependent variables below.

Cell-cell and cell-scaffold interactions
The functions and ψ describe respectively mechanical interactions between cells and between cells and their substrate; examples of relevant interactions include cell-cell and cell-substrate adhesion, and, in the case of motile cells such as fibroblasts, tractions between cells and their substrate.
Following Lemon et al. (2006) and O'Dea et al. (2010), we prescribe and ψ as follows: wherein ν, χ, δ a and δ b account for the cells' tendency to aggregate, their affinity for the substrate and the strength of the repulsive forces between cells and between cells and substrate. We remark that in this simplified formulation, cell-scaffold and cell-ECM interactions are lumped together in (10); however, as detailed below, we exploit our knowledge of the separate evolution of θ s and θ e to distinguish between cell-scaffold and cell-ECM interactions. As discussed in §1, cell adhesion to polymer scaffolds is mediated by adsorption of ECM molecules, or the embedding of specific cell recognition molecules, into the scaffold surface (Freed et al. 1994;Nikolovski and Mooney 2000). Since the type and density of such molecules can vary markedly, depending on the surface chemistry of the polymer, it is reasonable to expect that cells preferentially adhere to their extracellular matrices rather than synthetic substitutes. In this study, we aim to understand how such a disparity in cell-scaffold and cell-ECM interactions may influence the composition of the developing construct. We achieve this by comparing simulation results for the case χ = constant with those obtained for χ = χ(θ e ), so that the affinity of the cells for the scaffold phase depends on how much ECM has been deposited on the PLLA substrate; for simplicity, the remaining interaction terms ν, δ a and δ b are treated as fixed parameters.
We specify χ(θ e ) as follows: For suitable g 1 = constant, this represents a smoothed switch between the values χ 0 (the affinity between cells and scaffold in the absence of deposited ECM) and χ 1 > χ 0 (the elevated affinity due to ECM accumulation); θ e 0 is the threshold level of ECM deposition at which this change occurs. In the limit g 1 → ∞, the progression between χ 0 and χ 1 approximates a step function.

Mechanotransduction-mediated growth
The cells' response to biomechanical stimuli is accounted for by appropriate specification of the mass transfer rates S i (i = n, e, s); S w is chosen to ensure conservation of mass. Stimuli relevant to tissue engineering applications include contact inhibition, residual stress caused by tissue growth (Fung 1991;Skalak et al. 1996;Roose et al. 2003;Chaplain et al. 2006;Holzapfel and Ogden 2006) and the local fluid dynamics, such as the local hydrostatic pressure ; Klein-Nulend et al. represent the rate of cell proliferation and ECM deposition whilst k n d represents the rate of cell death. The thresholds are θ n 1 and θ n 2 1995) or fluid shear stress (You et al. 2000(You et al. , 2001Bakker et al. 2004b;Han et al. 2004;Yourek et al. 2004).
The focus of the current study is the interplay between scaffold degradation, tissue growth and cell-cell/cell-substrate interactions; for brevity, we restrict attention to cases in which cell proliferation and ECM deposition are regulated by contact inhibition. (Consideration of the influence of a wider range of mechanical stimuli is given in O' Dea et al. (2010) and .) We represent such regulation in our model by introducing net cell proliferation and ECM deposition rates κ n , κ e which depend upon the local cell volume fraction; the PLLA scaffold is assumed to degrade at a constant rate. We identify three distinct cellular responses to the local cell volume fraction: a proliferative phenotype, an ECM-depositing phenotype and an apoptotic phenotype, for each of which κ n , κ e take different values. The progression from one phenotype to the next occurs at threshold densities θ n 1 and θ n 2 (0 < θ n1 < θ n2 ). Functional forms for κ n and κ e are specified below and depicted in Fig. 3.
Equations (12)-(14) embody the following assumptions: (1) cells proliferate at a constant rate k n m until the cell volume fraction exceeds θ n 2 , when they become apoptotic (with death rate k n d ); (2) at intermediate values, the cells also deposit ECM at a rate k e m . As in Eq. (11), in the limit g 2 → ∞, κ n and κ e are piecewise-constant and the progression between phenotypes obeys a step function; in what follows, we choose g 2 = g 1 = 100.

Initial conditions
Tissue engineers employ a number of different techniques to seed porous scaffolds with cells. Here we consider static seeding, in which a suspension of cells is injected onto the surface of, or into, the scaffold, leading to an initial cell population which is localised near the point of injection. Following O'Dea et al. (2010) and , we prescribe: wheren s is the maximum initial cell volume fraction, and x = α and x = β represent the left and right hand edges of the localised cell distribution within the scaffold region (i.e. 0 < a < α < β < b < 1), near which g 3 governs the spatial gradient of θ n . In what follows, we fix a = 0.25, b = 0.75,n = 0.2, α = 0.4375, β = 0.5625, and g 3 = 50 without loss of generality.
Assumptions of uniform porosity have been employed in previous studies of in vitro tissue growth (Lemon et al. 2006;O'Dea et al. 2010;; however, the structure of scaffolds typically employed in such tissue engineering systems is highly heterogeneous (see Fig. 2). To determine the influence of such heterogeneity on construct evolution, we compare the predicted construct composition resulting from two separate initial conditions for the PLLA scaffold: We set θ ideal s = 0.0928, so that the former initial condition represents a scaffold of width (b − a) and uniform porosity ≈ 91%, this being the average initial porosity of the data shown in Fig. 2. The choice θ μCT s (x) denotes the spatially-varying initial experimental data shown in Fig. 2. In all cases, the initial ECM distribution is specified as

Numerical results
In this section, we present numerical simulations of Eqs.
The numerical scheme that we use and its validation are described in  and  and so we do not include details here. To summarise, Eq. (5) is solved for p w using θ n , θ s and θ e from the previous timestep (with a linear finite element approximation for p w ); employing this solution for p w , the remaining dependent variables, θ n , θ e and θ s are updated by applying to (2)-(4) an implicit Euler method, together with linear finite element approximations. For numerical convenience, we include a diffusive term, with constant diffusivity D (here we use D = 0.001), in Eq. (4) Table 1 and repeated in the relevant figure captions.

PLLA scaffold heterogeneity
To illustrate the model behaviour, in Fig. 4 we present simulations corresponding to the case for which cells are seeded within a PLLA scaffold of uniform porosity and their interactions with the PLLA scaffold and deposited ECM are identical; that is, we employ Eq. 16a to initialise θ s and choose χ = constant in place of (11). Figure 5 shows corresponding simulation results in the case for which θ s (x, 0) is specified via Eq. 16b.

(a) (b) (c) (d)
(e) (f) Fig. 4 Illustrative plots of the evolution of a the cell volume fraction (θ n ), b the scaffold phase volume fraction (θ s ), c the ECM volume fraction (θ e ), d the culture medium pressure ( p w ), e the axial culture medium velocity (u w ), and f the axial cell phase velocity (u n ) in the regime of cell volume fraction dependent growth and degradation and uniform cell-scaffold interaction properties at times t = 0 − 3 (in steps of t=0.5). The model parameters are as described in Table 1 except χ = χ 1 in place of (11). Initial conditions are given by (15), (16a) and (17). Vertical dotted lines indicate the end points of the scaffold at x = a, b Figures 4a and b reveal that the cells proliferate and deposit ECM asymmetrically under the influence of perfusion. The PLLA scaffold degrades uniformly so that θ s = θ s (t) in a x b. Near x = a, b where the cell volume fraction is small, PLLA degradation dominates ECM deposition and a decrease in total substrate fraction is observed; in regions of high cell volume fraction, ECM deposition ensures that the substrate density is maintained or increased. Figure 4a indicates that the cell volume fraction tends to a spatially-uniform value θ n = θ n 2 . As θ n reaches θ n 2 , the transition to the apoptotic phenotype precludes any further increase in cell density. Re-entry to the proliferative phenotype on subsequent reduction of cell density ensures a uniform distribution is maintained; the ECM deposition is greatest upstream due to perfusion causing advection of cells downstream, so that θ n attains the upper threshold at later times. Figure 4 indicates that, in the regime shown, the addition of spatial variation to the substrate phase via ECM deposition does not affect significantly the dynamics of the other model variables shown in Figs. 4d-f. This behaviour has been discussed in our previous work (see Figs. 4-7 of O'Dea et al. (2010)). Here, it suffices to note that the interplay between aggregative and repulsive cell behaviour [see Eq. (10)] is reflected in the evolution of the pressures and velocities of each phase (see Figs. 4d-f; the cell phase pressure is omitted for brevity). At low cell volume fraction (early times), aggregation dominates and the cell phase velocity indicates cell movement towards the centre of the population to form a dense aggregate; at high cell volume fraction (later times) repulsive effects dominate, leading to migration away from the central region; aggregative effects dominate at the edges of the population where the cell population remains low. Combined, these effects generate a dense cell aggregate, with steep spatial gradients of cell volume fraction near its up-and down-stream periphery (shown by the final line in Fig. 4a). To conserve mass, the culture medium phase moves in the opposite direction.
Comparison of Figs. 4 and 5 shows that, while the global features of the predicted cell and ECM volume fractions remain similar, the introduction of PLLA scaffold heterogeneity induces large, short-range variations in cell and ECM volume fractions and, additionally, has a significant influence on the cells' fluid-mechanical environment, with large variations in culture medium velocity evident over the short lengthscale associated with the heterogeneity of the substrate (Fig. 5e).
The importance of the scaffold (and ECM) distribution on the culture medium and cell phase velocities is evident from Equation (5), which shows that spatial gradients in influence the culture medium pressure gradient, leading to significant changes in culture medium flow. In addition, the axial cell phase velocity, obtained by integrating the axial component of the cell phase momentum equation (Equation (22), Appendix 4, in the long-wavelength limit), is defined as follows: which indicates that, under the influence of the final two terms, the large variations in substrate porosity present in Figs. 2 and 5b will lead to significant changes in cell velocity. (We also note that the singularity in θ s at x = a, b therefore leads to the spike in cell and culture medium velocities, illustrated in Figs. 4e, f and 5e, f. Additionally, inspection of Eqs. (10) and (18) Table 1, except χ = χ 1 in place of (11). Initial conditions are given by (15), (16b) and (17). Vertical dotted lines indicate the end points of the scaffold at x = a, b Spatial gradients of the cell and substrate phases also influence this behaviour: when ν and χ are constant, aggregation and attachment modify the cells' velocity via: so that when aggregation or attachment dominates repulsion, cells move up spatial gradients of θ n and ; in view of Eqs. (5), (10) and (18), spatial gradients of additionally modulate cell advection via the culture medium pressure. Given the importance of culture medium flow-induced mechanical stimulation to the growth and differentiation of various cell types (see §1 and §2.3), such variations in cell and culture medium velocity are likely to have a significant effect on local cell behaviour. We pause to remark that Figs. 4a and 5a show that as the cell population expands to colonise the scaffold, our fluid-based model allows egress of cells (and, eventually, the subsequently deposited ECM) from the scaffold into the up-and downstream regions (x < a and x > b). Lemon and King (2007) demonstrate that, due to cell-scaffold adhesion, such behaviour is minimised in scaffolds whose density distributions decay to zero at the boundaries; we have, however, employed experimentally-relevant data to initialise θ s here. Cell egress is not desirable in the current context of a perfusion bioreactor (and is in any case minimal in the seeding protocol employed here: the cell flux θ n u n is small); however, it is of biological relevance to modelling tissue invasion after implantation.

Cell-substrate interactions
In this subsection we investigate the influence of the properties of PLLA scaffold on the eventual construct composition, both in terms of its spatial distribution and its interactions with the cells.
First, we compare simulation results for which χ = constant with those for which cell-scaffold adherence is governed by (11) and, in so doing, investigate how a disparity between cell-scaffold and cell-ECM attachment strength (due to the surface chemistry of the polymer scaffold) may influence the composition of the developing construct. For clarity, we employ Eq. (16a) to initialise θ s , corresponding to a PLLA scaffold of initially uniform porosity; the combined influence of cell-substrate interactions and scaffold heterogeneity on construct composition is investigated subsequently.
In Figs. 6a and b, we assume that the cells' affinities for the PLLA scaffold and the deposited ECM are identical; the figures show cell and scaffold phase distributions at illustrative time points for different values of χ . Comparison of the construct composition when χ = 0 and χ = 5, suggests that cell migration up spatial gradients of θ s leads to a more sharply-defined cell volume fraction profile; however, the maximal volume fraction is reduced. This is because cell movement is confined to the edge of the cell aggregate (at the centre ∂ /∂x ≈ ∂θ n /∂ x ≈ 0).
Inspection of the cell volume fraction distribution at later times reveals more interesting behaviour. As θ n increases, advection of cells becomes more significant, leading to profiles which are skewed in the downstream direction. For large cell volume fraction, the combination of advection and cell-cell and cell-substrate repulsion leads to cell migration away from the aggregate's centre. This outward drift is balanced by inward movement of cells at the periphery of the densely populated region. In the  (15), (16a) and (17). Except as stated, all parameters are as given in Table 1. Vertical dotted lines indicate the scaffold periphery at x = a and x = b case of strongly adherent cells, increased (inward) movement up spatial gradients of causes the cell volume fraction to peak at the periphery of the aggregate, with a flatter profile near the centre (see the final profile in Fig. 6a in the case χ = 5). For small values of χ , a profile similar to that shown in Fig. 4a is obtained. Figures 6c and d show how the construct composition is influenced by differential adhesion between the porous scaffold and the deposited ECM. The cells' adherent behaviour is modelled by Eq. (11) so that χ varies between χ 0 = 0 and χ 1 = 5 in response to the ECM volume fraction. These simulations indicate that in the case of a uniform initial scaffold, preferential adherence to ECM leads to the creation of a tissue construct whose composition is indistinguishable from that obtained when the adherence properties of the substrate are uniform. Corresponding results were obtained when the initial scaffold porosity was spatially-nonuniform, but are not included. In view of this somewhat surprising outcome, we conclude that the mathematical complexity (a) (b) (c) (d) Fig. 7 The evolution of the cell volume fraction (t = 0, 1, 2, 3) and the ECM volume fraction (t = 2, 3) in the case of cell volume fraction dependent cell proliferation and ECM deposition defined by Eqs. (12) and (14); arrows indicate the direction of increasing time. Plots are for χ = 1 (solid line) and χ = 0 (dashed line). In a, b initial conditions are given by (15), (16a) and (17) and in c, d Eq. (16b) is employed to specify θ s (x, 0). Except as stated, all parameters are as given in Table 1. Vertical dotted lines indicate the end points of the scaffold at x = a, b associated with the biologically-inspired cell-substrate interactions [and embodied by Eq. (11)] adds no additional predictive capability to the model. Henceforth, we return to the simplified model in which χ = constant. We note that in this case, inspection of Eqs. (5) and (10) reveals that the influence of the cell-scaffold affinity strength on the model behaviour may only be studied in the presence of a spatially-heterogeneous substrate volume fraction; our previous work (O'Dea et al. 2010) therefore neglected its influence. In Fig. 7 we demonstrate how the combined effects of cell-substrate interactions (χ = constant) and the heterogeneity of PLLA scaffold volume fraction influences the construct composition. As indicated by Eq. (18), strong cell-scaffold adherence enhances cell movement up scaffold and ECM gradients. Large, short-range spatial gradients of scaffold porosity exist in the experimental data we have employed to initialise θ s . We therefore observe large deviations in cell and ECM distributions, which mirror, and exaggerate, the underlying PLLA scaffold porosity distribution. Since the deposition of ECM (and other associated extracellular materials) enables the maintenance of the mechanical properties of the degrading PLLA scaffold, such  (15), (16b) and (17). Except as stated, all parameters are as given in Table 1 heterogeneous ECM distributions have important implications regarding the structural suitability and suitability for implant of the resulting construct. With this in mind, in the following subsection we focus on the interplay between PLLA degradation and ECM deposition and the maintenance of substrate porosity during this process.

Scaffold degradation
The maintenance of tissue construct material properties is crucially affected by achieving a match between the rates of scaffold degradation and deposition of ECM and other extracellular materials. In this section we indicate how our model may be employed to determine how the evolution of the substrate volume fraction depends on the model parameters. For clarity, in preference to the spatio-temporal distributions presented in § §3.1,3.2, we consider the evolution of the total substrate mass (t), and its PLLA and ECM components θ s (t), θ e (t), which are defined as follows: The simulations presented in Fig. 8 demonstrate that close control of scaffold degradation and ECM deposition is required in order to maintain substrate density. We do not present the corresonding cell volume fraction evolution since our focus here is on the tissue construct's rigid, load-bearing components. In Fig. 8a, we consider a nondegrading scaffold, and the substrate volume fraction is therefore determined by the evolution of the ECM phase. With the addition of scaffold degradation (Fig. 8b), we observe an initial decrease in substrate due to scaffold degradation; as the seeded cell population increases, Eq. (14) leads to progression to an ECM-depositing phenotype and the substrate volume fraction increases. Our model predicts that, eventually, the cell volume fraction will increase to a point at which all cells in the scaffold enter apoptosis, at which point the substrate volume fraction achieves an equilibrium level. However, the timescale over which we perform our simulations is restricted by our wish to restrict tissue egress from the scaffold into the up-and downstream regions x < a, x > b (see §3) and we therefore do not observe such equilibrium behaviour. Nevertheless, our simulation results indicate that the rate of scaffold degradation is a key experimental variable, and suggest that there is a threshold time, before which the construct's mechanical properties are likely to be unsuitable for implantation.

Discussion
In this paper we have presented a multiphase model describing tissue growth within a perfusion bioreactor, modelled as a two-dimensional channel containing cells, culture medium, a porous PLLA scaffold and deposited ECM. The formulation employed is based on the general multiphase formulation proposed in Lemon et al. (2006) and extends the three phase model of O'Dea et al. (2010), where the PLLA scaffold was assumed to be spatially-homogeneous and inert. Many similar studies of tissue growth have employed this simplifying assumption, tacitly assuming that the importance of spatial variation of scaffold/ECM volume fraction is negligible. Here, we include additional mass conservation equations for the scaffold and ECM phases with which to model scaffold degradation and ECM deposition. This allows us to incorporate PLLA scaffold phase heterogeneity, inhomogeneous deposition of ECM, and to consider explicitly the interactions between cells and their different supporting structures, while remaining within a simplified modelling framework.
Comparison of our simulation results with our previous work  indicates that the predicted cell volume fraction and variables related to the mechanical environment are strongly affected by heterogeneous scaffold volume fraction distributions; for instance, the culture medium velocity shows strong short-range variation, which is likely to have a profound effect on the mechanical environment of the cells. Additionally, we have demonstrated that ECM deposition by the cells is highly localised in the regions of elevated cell volume fraction and that spatial variation in the PLLA scaffold volume fraction leads to large deviations in cell and ECM distributions.
In addition to spatial variation of the scaffold volume fraction, and motivated by the possible variation in cell binding sites in tissue engineering scaffolds, we considered the heterogeneous interactions between cells and their supporting scaffold and the deposited ECM, postulating that cells have a greater affinity for ECM (and other deposited materials) than the scaffold onto which they are deposited. Employing a smoothed switch between low and high affinity in response to the local ECM volume fraction, we showed that such a disparity has no significant influence on eventual construct composition. This conclusion indicates that simplified models in which cells interact uniformly with their supporting structures may be employed without affecting their biological relevance (in the present study, this corresponds to χ = constant). We showed that in such a simplified model cell-scaffold/ECM interactions can have dramatic affects on the construct composition, leading to significantly enhanced migration of cells up spatial gradients of substrate volume fraction. In the case of experimentallyrelevant initial PLLA scaffold porosity distributions, which display large gradients in scaffold distribution, this leads to cells and ECM profiles which mirror and exaggerate the underlying scaffold porosity. The heterogeneity of the resulting tissue construct has important ramifications for its structural stability and suitability for implant.
Our model associates scaffold degradation with a reduction in scaffold phase volume fraction, from which we infer deleterious effects on material properties, which are ameliorated by deposition of ECM (and other extracellular materials). The interplay between scaffold degradation and ECM deposition is therefore crucial in determining tissue construct material properties. The results of our simulations discussed above indicate that the production of scaffolds with uniform porosity may play an important role in producing tissue constructs with mechanical properties appropriate for implant.
We have extended our previous studies (O'Dea et al. 2008, by accommodating spatial non-uniformity in scaffold and ECM volume fraction, demonstrating the importance of such a consideration; however, we have made a number of simplifying assumptions to enable analysis. We have restricted attention to a rigid scaffold and ECM phase (the remaining phases being modelled as viscous fluids) and so our formulation applies to those constructs whose solid characteristics are dominated by the rigidity of the scaffold and/or deposited materials. We note also that our simplified treatment of the PLLA scaffold and ECM means that the mechanics of our model are in fact accommodated within a three phase framework. Our model formulation is further simplified by exploiting the long-wavelength limit; our previous work ) has indicated that, while two-dimensional variation of mechanical stimulation has an important effect on construct growth, the long-wavelength limit provides a good approximation to the averaged behaviour of the two-dimensional model, even for bioreactors with large aspect ratio, such as that illustrated in Fig. 1. Nevertheless, validation of our model results within a two-dimensional framework remains important future work. We have assumed that the scaffold degrades uniformly; however, it is known that bi-products of tissue growth can influence the mechanical and structural nature of tissue engineering scaffolds (see, e.g., Ahearne et al. (2010) and references therein). Consideration of such effects represents another interesting extension to our study.
We have employed our simplified formulation to make inferences regarding the likely composition and mechanical integrity of engineered tissue constructs. Important further work includes explicit modelling of the mechanical properties of the scaffold and ECM, considering, for example, a poroelastic (Roose et al. 2003) or poroviscoelastic (Byrne and Preziosi 2003) model. Additionally, the robustness of our conclusions should be investigated by considering the addition of nutrient-limited growth to our formulation; methodologies for such investigations are given by Lewis et al. (2005) and Lemon and King (2007). Tissue engineers employ a number of different techniques to seed porous scaffolds with cells. We have restricted attention to an initial seeding which is highly localised near the centre of the scaffold (via, e.g. injection of a cell suspension), to illustrate the model behaviour. With the inclusion of nutrient-limited growth, it is natural to investigate the influence of different experimentally-relevant cell seedings on the eventual construct morphology. Such studies allow for further incorporation of experimental data (via the initial seeding density, in addition to the scaffold porosity distribution) and lend themselves for validation against histological samples from relevant cell culture experiments.
Acknowledgments This work was supported by funding from the EPSRC in the form of a Ph.D studentship (RDO) and an Advanced Research Fellowship (SLW). It was also supported by the EPSRC/BBSRC-funded OCISB project BB/D020190/1 (JMO), and based on work supported in part by Award KUK-013-04 made by King Abdullah University of Science and Technology (HMB). We are also grateful to E. Baas, ISTM, Keele University for the provision of experimental data.

Appendix: Model derivation
We consider a bioreactor of length L * and width h * , modelled as a two-dimensional channel containing a mixture of four interacting phases, representing cells, culture medium, PLLA scaffold and ECM and denote these via a subscript i = n, w, s, e, respectively. The viscosity of any fluid phase is denoted μ * i , and the typical timescale for tissue growth (comprising both cell proliferation and ECM deposition) is denoted K * . Asterisks distinguish dimensional quantities from their dimensionless equivalents.
We introduce a Cartesian coordinate system L * x = L * (x, y) and time K * t and the channel occupies the dimensionless region 0 x 1, 0 y h = h * /L * . The volume fraction of each phase is denoted θ i , while the dimensionless volumeaveraged velocities, pressures and stress tensors of the each phase are denoted K * L * u i = K * L * (u i , v i ), K * μ * w p i and K * μ * w σ i . Tissue growth, scaffold degradation and ECM deposition are captured via material transfer functions K * S i . We assume that all dimensionless dependent variables are functions of x and t.
The model is constructed by considering mass and momentum balances for each phase, assuming that each phase is incompressible, with equal density, and neglecting inertial effects; the equations governing the ith phase (with volume fraction θ i ) are as follows [see Lemon et al. (2006) Additional conservation conditions may be obtained by summing over all phases and exploiting the no-voids condition i θ i = 1. In Eq. (21) K * S i is the net material production term associated with phase i (mass conservation demands that S i = 0); in (22), K * μ * w /L * F i j is the interphase force exerted by phase j on phase i, obeying F i j = −F ji . These interphase forces comprise interphase viscous drag (with drag coefficient μ * w /L * 2 k) and active forces, the latter being embodied within extra pressures which arise due to cell-cell, cell-ECM and cellscaffold interactions; interactions between the culture medium and scaffold phases are assumed to involve only viscous drag. The mechanics of this four phase formulation is simplified by lumping the scaffold and ECM components into a single 'substrate' phase, denoted θ S = θ s + θ e , and modelled as a rigid porous material. For notational convenience, in this Appendix, we employ the subscript S to denote the substrate, in preference to . Separate mass conservation equations are nevertheless employed for θ s and θ e to track their individual evolution. Fig. 9 The two-dimensional domain, outward-pointing normaln and associated boundary conditions. The arrows indicate the perfusion direction in the case of dimensionless up-and downstream pressures P U and P D obeying P U > P D The cell population and culture medium are represented as distinct viscous fluids, modelled by standard viscous stress tensors; the rigidity of the substrate implies u S = 0. These constitutive assumptions are embodied in the following equations.
p n = p w + θ 2 n n + θ n ψ nS , wherein μ i are the dimensionless viscosities of each phase, and n and ψ nS are defined n = −ν + δ a θ n θ w and ψ nS = −χ + δ b θ n θ w .
In Eq. (26) ν, χ, δ a , δ b > 0 dictate the cells' tendency to aggregate, their affinity for the scaffold/ECM and the strength of cell-cell/cell-scaffold repulsion. In a more general formulation, the coefficient of viscous drag k between two phases i and j varies depending upon the phases under consideration and may depend upon their respective volume fractions or other state variables. A suitable representation is to replace k by, say, k i j (θ i , θ j ) (obeying k i j = k ji ). Full details and discussion of the above choice of interphase interaction terms may be found in Lemon et al. (2006). Figure 9 depicts the two-dimensional model of the bioreactor, together with appropriate boundary conditions. These correspond to no-slip and no-penetration of cells or culture medium through the channel walls, a pressure-driven flow imposed via up-and downstream pressures K * μ * w P U and K * μ * w P D , partitioned normal stress conditions and fully-developed flow at x = 0, 1.
We now simplify the two-dimensional equations by considering the limit for which the aspect ratio of the bioreactor is asymptotically small (h 1). We remark that, since the culture medium volume fraction may be eliminated via θ w = 1 −θ n −θ s −θ e and the substrate is rigid, we need consider momentum conservation equations for the fluid (cell and culture medium) phases, only.