A two-fluid model for tissue growth within a dynamic flow environment

We study the growth of a tissue construct in a perfusion bioreactor, focussing on its response to the mechanical environment. The bioreactor system is modelled as a two-dimensional channel containing a tissue construct through which a flow of culture medium is driven. We employ a multiphase formulation of the type presented by G. Lemon, J. King, H. Byrne, O. Jensen and K. Shakesheff in their study (Multiphase modelling of tissue growth using the theory of mixtures. J. Math. Biol. 52(2), 2006, 571–594) restricted to two interacting fluid phases, representing a cell population (and attendant extracellular matrix) and a culture medium, and employ the simplifying limit of large interphase viscous drag after S. Franks in her study (Mathematical Modelling of Tumour Growth and Stability. Ph.D. Thesis, University of Nottingham, UK, 2002) and S. Franks and J. King in their study (Interactions between a uniformly proliferating tumour and its surrounding: Uniform material properties. Math. Med. Biol. 20, 2003, 47–89). The novel aspects of this study are: (i) the investigation of the effect of an imposed flow on the growth of the tissue construct, and (ii) the inclusion of a mechanotransduction mechanism regulating the response of the cells to the local mechanical environment. Specifically, we consider the response of the cells to their local density and the culture medium pressure. As such, this study forms the first step towards a general multiphase formulation that incorporates the effect of mechanotransduction on the growth and morphology of a tissue construct. The model is analysed using analytic and numerical techniques, the results of which illustrate the potential use of the model to predict the dominant regulatory stimuli in a cell population.


Introduction
In vitro tissue engineering is the logical extension of transplant surgery and involves the growth of replacement tissue outside the body to alleviate the chronic shortage of tissue available from donors [15]. Significant research has been dedicated to the study of tissue engineering processes relevant to liver, skin, collagen (for cartilage/tendon replacement) and bone (see [15,60] for interesting reviews) and progress has been made in a number of areas. These include: maintenance of tissue-specific functions via mimicry of in vivo conditions through appropriate cell co-culture and/or three-dimensional spheroidal culture [3,26,28,41,55]; and understanding the influence of artificial scaffolds which provide mechanical support for construct growth and whose surface chemistry and pore size may be altered to encourage cell anchorage or increased population of the scaffold [30,58,65].
Three-dimensional tissue culture is desirable to maintain cell-specific function; however, limitations in the diffusion of nutrients and waste products often result in the formation of a construct with viable, proliferating cells at its periphery but a necrotic core: early studies have shown that cellular spheroids of diameter greater than 1 mm generally develop a necrotic core [61]. Depending upon the application, engineered constructs must be relatively large to serve as grafts for tissue replacement; it is clear, therefore, that mass-transfer limitations represent challenges for in vitro tissue engineers [45]. To rectify this, bioreactors are widely used. A bioreactor is an advanced tissue-culture apparatus which enables control of the culture environment. Nutrient supply and metabolite removal may be enhanced by using the perfusion or circulation/mixing strategies; furthermore, bioreactors allow monitoring and control of factors such as pH, and the provision of growth factors and other cell-signalling molecules.
Bioreactors also allow controlling the biochemical environment, many bioreactors are designed specifically to provide mechanical stimulation to cell cultures via, for instance, fluid shear stress or tensile or compressive forces applied either on the macroscale or via magnetic particles embedded in the cell membrane (see [8,45] for a review). These stimuli are integrated into the cellular response via a process known as mechanotransduction. The importance of mechanical stimuli to tissue function is noted by many authors including Fung [24], who asserts that the correct function of organs is dependent on the level of internal stress. A great many "proof of principle" studies have illustrated the beneficial effect of mechanical conditioning on the structural and functional properties of engineered tissues; however, little is known about either the manner in which these forces should be applied for specific tissues or how these stimuli are interpreted by the cells [45]. It is clear, however, that the mechanical environment required for optimum growth will be peculiar to the tissue under consideration; bespoke bioreactors are therefore required to provide appropriate physical (and biochemical) cues for different tissue engineering applications.
Mathematical modelling of tissue growth is a subject area which has received a great deal of attention; the clinical applications of such models are myriad. In particular, much work has been dedicated to modelling tumour growth [1], angiogenesis [9,11], wound healing [59] and, more recently, in vitro tissue engineering processes [13,14]. Within the context of tissue engineering, the goals of such models are to explain how observed tissue engineering problems arise, to suggest mechanisms to resolve them and, thereby, to predict optimal protocols for tissue growth. Furthermore, idealised tissue growth studies can provide insights useful in the design of bespoke bioreactor systems.
A variety of approaches has been used to model tissue growth. Discrete cellular automata models are employed in, for example, [20,49,50], while in [31,32] chemotaxisbased continuum models are exploited (see [40] for a review). Numerous studies have considered the effect of limited nutrient availability and/or the presence of inhibitors on tissue growth (e.g. in the context of tumour growth [5,25,42]; for a review see [1]). The stresses experienced by cells at the microscopic level are calculated in [29,46,62]; the influence of the cells' mechanical environment on growth is considered by Roose [57] in the context of growth-induced stresses and, scaffold adhesion-guided behaviour is investigated by Powers et al. [52,53].
The preceding studies analyse tissue growth processes in general terms with the results applied to (for instance) tumour growth; studies which consider specifically tissue growth in porous scaffolds include Malda et al. [43] in which the development of oxygen gradients in the absence of perfusion is investigated using a simple diffusion-consumption model. Parameter estimation is achieved via comparison with experimental evidence. Three-dimensional fluid flow through porous scaffolds in a perfusion bioreactor is studied by Porter et al. [51] in which a detailed model of a porous scaffold is obtained via micro-computed tomography imaging and the flow profile calculated using the Lattice-Boltzmann method. Relating simulation results to experimental results, it is concluded that a mean pore-surface shear stress of 5 × 10 −5 Pa corresponds to an increased cell proliferation and viability. Raimondi et al. [54] demonstrate that the material properties and cell viability of constructs resulting from perfusion shows a two-fold improvement compared to surface perfusion or static culture; computational modelling is used to quantify the fluid-dynamical environment at the microscopic level. Modelling of both the cell growth and fluid flow within a three-dimensional scaffold in a perfusion bioreactor is considered by Colletti et al. [12]. The flow through the scaffold is governed by Brinkman's equation and the nutrient distribution is governed by a reaction-advection-diffusion equation. Cell growth is assumed to depend upon local nutrient availability via an ordinary differential equation.
In this paper, we present a continuum mathematical model relevant to tissue growth processes in a perfusion bioreactor. In contrast to the above analyses of perfusion bioreactors [12,43,51,54], we employ a multiphase formulation after Lemon et al. [37], restricted to two interacting viscous fluid phases, representing a cell population and attendant extracellular matrix (ECM) and a culture medium (fluid-based models for biological tissue growth have been widely exploited, especially in modelling tumour growth; see, for example, [6,7,21,23]). We consider this model to be a first step towards a general multiphase formulation allowing examination of the effect of mechanotransduction on the growth and subsequent morphology of a tissue. Our continuum macroscale model allows explicit consideration of the complex interactions involved in tissue growth without considering the precise microscopic detail; however, the averaging process involved in deriving models of this type ensures that the terms present in the model equations arise from appropriate microscopic considerations. The derivation and analysis of multiphase models (for biological and other applications) has been treated in great detail in the literature; see, for example, [2,35,44,63] and for two phase flows [16][17][18].
The perfusion bioreactor under consideration is based upon that employed by El-Haj et al. (see, for example, [19]) which comprises a tissue construct within a culture mediumfilled cylinder along which a flow is driven; however, the mathematical techniques and modelling approaches employed here are readily transferrable to other tissue culture systems. Our model accommodates mechanotransduction-affected cell proliferation, ECM deposition and cell death. Nutrient-limited growth is not considered since we assume that perfusion provides an abundant supply of nutrient. The investigation of the effect of the imposed flow on the response of the cell phase and the inclusion of a simple mechanotransduction mechanism are the novel aspects of this study.
In §2, we introduce the perfusion bioreactor system and present the equations governing our two phase model. Guided by a parameter estimation, in §2.1, we exploit the limit of large interphase viscous drag to recover the model equations of [21,23], which we then employ for the remainder of this paper. In §3, we investigate the impact of an ambient flow on the growth of a one-dimensional tissue. Analytic solutions are presented in the limit for which the tissue is defined by two sharp interfaces, and the stability of these interfaces to transverse perturbations is investigated ( §3.1). Numerical solutions are obtained for a diffuse one-dimensional cell population growing at a constant rate ( §3.2) and in §4.1, corresponding two-dimensional numerical simulations are presented. Lastly, in §4.2 we extend the model by coupling the cells' proliferative response to the local mechanical environment, allowing consideration of a simple mechanotransduction mechanism. A discussion of our model and its applications within a tissue engineering context are presented in §5, together with suggestions for future work.

A two-phase model for tissue growth
We apply the general multiphase formulation given in [37] to develop a simple model relevant to tissue engineering processes. For brevity, we do not recapitulate the multiphase formulation here: the interested reader is directed to [17,35,37] for more details.
We consider the growth of a tissue construct within a nutrient-rich fluid culture medium and investigate the effect of an imposed flow on the response of the cells. We neglect the solid characteristics of the tissue construct (including, for instance, the presence of a scaffold) and employ a two-fluid model of the type presented in [4,6,17,21,23,33]. This system is representative of tissue culture within a dynamic flow environment, for example, a perfusion-type bioreactor system. Our idealised perfusion bioreactor model is based upon that employed by El-Haj et al. and we represent this system as a two-dimensional channel containing a two-phase mixture of interacting viscous fluids. A two-dimensional Cartesian geometry is chosen for simplicity; however, generalisation to a cylindrical geometry is straightforward. The cells and ECM are modelled as a single phase (henceforth denoted the "cell phase"); the second phase represents the culture medium. Perfusion is represented by an imposed flow of culture medium. In the following we will employ the term "tissue construct" to distinguish the region occupied by the interacting system of cell and culture medium phases from the remainder of the channel which contains only culture medium. At certain stages in this paper, the interface between the construct and the surrounding culture medium will be assumed to be sharp or diffuse to allow different analyses to be undertaken. A Cartesian coordinate system x * = (x * , y * ) is chosen with corresponding coordinate directions (x,ŷ) and the channel occupies 0 6 x * 6 L * , 0 6 y * 6 h * (see Figure 1). In this paper, asterisks distinguish dimensional quantities from their dimensionless equivalents. We associate with each of the cell and culture medium phases a volume fraction, denoted by n and w, respectively, and a volume-averaged velocity, u * i = (u * i , v * i ), pressure, p * i , stress tensor, σ * i and density, ρ * i (where i = n, w denotes variables associated with the cell and culture medium phases, respectively) and assume that these are functions of x * and t * , where t * represents time. Table 1 summarises the variables employed in the model, together with their units (where appropriate).  Assuming that each phase is incompressible with the same density, and neglecting inertial effects, we obtain the following governing equations (see [37]): conservation of mass: ∂n ∂t * + ∇ * · (nu * n ) = S * n + D * ∇ * 2 n, (2.1) conservation of momentum: no-voids: In equations (2.1) and (2.2), S * n , S * w are the net rates of material production for each phase and D * is the diffusion coefficient for the cell phase. In (2.3), F * wn represents the force exerted by the cell phase on the culture medium at the interface between these phases, and we assume F * wn = −F * nw as implied by (2.4). As in Drew [17], the influence of these interfaces, especially with respect to thermodynamic properties, is not considered; see Kolev [35] for a thorough discussion of these terms.
We remark that a diffusive term has been added to the mass conservation equation (2.1); whilst cells do exhibit random motion, the tissue growth and perfusion-induced flow fields are the dominant mechanisms leading to cell movement, with diffusive effects assumed to be negligible [22,33]. However, we retain diffusive terms for numerical convenience since they eliminate the moving boundaries between the tissue construct and the culture medium, ensuring that we need not track explicitly the sharp interface which is evident when D * = 0.
To close the model equations (2.1)-(2.5), we now specify constitutive laws for the interphase forces (F * ij ), stress tensors (σ * i ) and material production rates (S * i ) which describe the behaviour of our system.
Following [4,6,37], we represent the interphase force as follows: wherein the pressure in the culture medium phase is denoted p * w , and k * is the coefficient of interphase viscous drag, which we assume is constant. The pressure in the cell phase is related to that in the culture medium by the following relation: where Σ * n is a prescribed intraphase pressure resulting from interactions within the cell phase such as osmotic stresses or surface tension within cell membranes. Equation (2.7) will be employed to eliminate the cell phase pressure, p * n , from the governing equations and carries the tacit assumption that the interphase tractions are negligible. We do not specify a functional form for the intraphase pressure function, Σ * n , since in the analysis presented in this paper, we employ a lumped pressure field which encapsulates this function. Appropriate functional forms for intraphase pressures and interphase tractions are given in [6,37]. Each of the phases are modelled as viscous fluids and the associated volume-averaged stress tensors are therefore where μ * i and λ * i are the dynamic shear and bulk viscosity coefficients of the ith phase and I is the identity matrix. Functional forms for the material production terms S * i will be chosen to reflect different physical processes at various points in this paper and we delay specification of appropriate boundary conditions until the derivation of a model in the limit of large interphase viscous drag.
We non-dimensionalise as follows: wherein U * w is a velocity scale which will be determined by the choice of pressure drop or upstream flux. A viscous scaling is employed for the culture medium and intraphase pressures (p w , Σ n ) since we assume that viscous effects dominate inertia in the momentum equations. In view of the non-dimensionalisation t * = L * t/U * w , the timescale of interest is the time taken for a fluid particle to travel along the length of the bioreactor. The flow rate for the bioreactor system of El-Haj et al. is 0.1 ml/min and the construct has a diameter of 9 mm and length 4 mm; assuming that this is of the same order as the bioreactor length, L * , this gives a timescale of approximately 2.5 minutes. We note that this is very short in comparison to the timescale over which tissue growth occurs; however, in this paper, we consider the case for which the ratio of the growth and flow timescales is O(1) (which corresponds to employing fast growth rates) to minimise computation time and to illustrate features of the system.
The dimensionless form of equations (2.1)-(2.4) is as follows: and equation (2.5) is used to eliminate w.
In equations (2.10)-(2.13), the dimensionless parameters D, μ n , k, γ w and γ n are defined as follows: (2.14) and the channel now occupies 0 6 x 6 1, 0 6 y 6 h = h * /L * . The physical interpretation of the dimensionless diffusion coefficient (or inverse Peclet number) D, relative viscosity μ n and drag coefficient k is obvious. The parameter γ i describes the relative importance of the viscosity associated with the rate of change of volume of the ith phase compared to that associated with fluid shear. It is usual to take λ * i = −2μ * i /3 (so that in equation (2.8), we have p i = −σ i,kk /3; i = n, w) implying γ w = −2/3 and γ n = −2μ * n /3μ * w [21,23,33,37]. We expect μ * n > μ * w and therefore assume γ n 6 −2/3.

Large drag limit
By considering the Stokes drag due to water flowing past a spherical cell, Lubkin & Jackson [39] estimated the drag coefficient as: k * = 4.5 × 10 7 N · m −4 · s. An appropriate lengthscale for the system is the bioreactor length, and assuming that this is of the same order as the scaffold length, we choose L * = 4×10 −3 m. Assuming that the culture medium viscosity is equal to that of water, μ * w = 1 × 10 −3 N · m −2 · s [36], the dimensionless drag coefficient may then be calculated from equation (2.14) to give Motivated by this, we now consider the limit in which the interphase viscous drag is large (k 1) and derive appropriate simplified versions of the governing equations and boundary conditions. We consider a power series expansion of the dependent variables as follows: with similar expansions for u n , u w , p w and Σ n ; at leading order, equation (2.12) gives u n0 = u w0 = u 0 so that we may associate with each phase a common velocity. Following [21,23], we assume that each phase has the same material properties (γ i = −2/3, μ n = 1) and that at leading order the mass production terms S n0 , S w0 are given by in which k m = (k * m L * )/U * w represents the dimensionless rate of cell mitosis and ECM deposition (in the interests of brevity we henceforth refer to this as the "growth rate"), k d = (k * d L * )/U * w represents the dimensionless rate of cell death and ECM degradation (which we will refer to as the "death rate") and k * m and k * d are the corresponding dimensional rates. In general, these will depend upon the cells' mechanochemical environment (e.g. nutrient availability, growth factors, local cell density or stress) and in the following derivation we therefore allow spatio-temporal variation: k m (x, t), k d (x, t).
Dropping the subscript notation for brevity, from (2.11)-(2.13), we obtain the following equations at leading order: where the lumped pressure field, p, is related to p w0 via p = p w0 + n 0 Σ n0 . Taking the divergence of equation (2.19) gives an expression for ∇ 2 (∇ · u), which, on substitution into the Laplacian of (2.18), yields an equation for the pressure in terms of the cell volume fraction only. We may then restate our governing equations in the more natural form given in [21,23]: This model breaks down in a boundary layer of thickness O(1/ √ k) near the channel walls when intraphase viscous effects become important (see equation (2.12)); appropriate boundary conditions on the outer problem are given in Franks [21] as follows: which imply slip along, but no-penetration of cells or culture medium through y = 0, h. Axial boundary conditions which ensure that the tissue construct does not extend along the channel's length and, additionally, set an axial pressure drop which drives a flow are as follows: wherein P u and P d are the dimensionless up-and downstream pressures, related to the corresponding dimensional pressures via (P u , P d ) = (P * u , P * d )L * /(U * w μ * w ). In the following sections, we present solutions to the model equations (2.20)-(2.25) in various limits.
In an extension to [21,23], in §3 and §4.1 we investigate the effect of an ambient flow on a uniformly proliferating tissue construct, for which k m and k d are constant. In §3, solutions are obtained in the one-dimensional limit: in §3.1 analytic solutions are presented in the regime for which the tissue construct is delineated by two sharp interfaces (D = 0) and its stability to transverse perturbations is determined; in §3.2, numerical solutions for a diffuse cell population (D 0) are presented and compared with the analytic solutions from §3.1. In §4.1, corresponding two-dimensional numerical solutions are presented. In §4.2, we further extend the model by postulating functional forms for the growth rate, k m (x, t), which allow the influence of a range of mechanical stimuli on the growth response of the cells to be accommodated. We remark here that in all of the subsequent numerical simulations, the parameter values are selected to illustrate the behaviour of the model under a particular growth regime. The parameter values employed in each of the following sections are summarised in Table 2.

One-dimensional growth
We now assume that the tissue undergoes one-dimensional growth parallel to the x-axis and that the associated pressure and velocity fields are functions of x and t only. For mathematical convenience, we initially consider growth within a channel of semi-infinite length. In this case, equations (2.20)-(2.22) reduce to give and we emphasise that we consider uniform growth, for which k m and k d are constant.
The axial boundary conditions presented in §2.1 require modification since the channel is now semi-infinite in extent. Integration of equation ( where α(t), β(t) are arbitrary functions of time; application of the boundary conditions p = P u , n = 0 at x = 0, p = P d , n = 0 as x → ∞ indicates that we require α = 0, β = P u = P d . A pressure drop-induced flow may therefore not be imposed; instead, we impose an upstream flow, U * , corresponding to u(x = 0, t) = U, where U = U * /U * w is the dimensionless upstream flowspeed. In the following, we choose U = 1, which sets the velocity scale to be U * w = U * . The value of the up-and downstream pressure is arbitrary and we choose P u = 0 = P d without loss of generality. We remark that the lengthscale, L * , is now defined by our choice of growth rate according to k m = k * m L * /U * w , so that the lengthscale of interest is the distance travelled by a fluid particle over the growth timescale. Appropriate boundary conditions are 3) (with α = β = 0) allows us to obtain a reduced model in terms of the cell volume fraction (n) and axial velocity (u) only (details omitted).

Sharp interface limit: D = 0
We now consider the regime in which the interfaces between the tissue construct and surrounding culture medium are sharp, corresponding to D = 0. The domain is then decomposed into three distinct regions by planar interfaces. We denote the interfacial positions by x = L(t), R(t), across which we impose continuity of velocity and normal stress In (3.6), we have adopted the notation [..] + − to denote the jump across an interface, the superscript + indicating the limiting value x = L(t) or x = R(t) from within L(t) 6 x 6 R(t). The evolution of the interfacial positions x = L(t), R(t) is determined from the following kinematic conditions which ensure that particles on the interfaces remain there: For simplicity we consider a growing construct of uniform density, represented as follows: It is then straightforward to integrate equations (3.2a, b) to determine The evolution of the cell distribution is determined from equation (3.1), which yields the following well-known logistic growth behaviour: , (3.12) in which n(0) is the initial cell density, r = k m −k d is the net growth rate and K = 1−k d /k m is the carrying capacity. Using (3.7), we deduce that the interfacial positions of the growing tissue construct can be written where L(0) and R(0) are the initial interfacial positions. These solutions show that the construct is advected with the imposed flow (u(x = 0, t) = 1); furthermore, its width, given by R(t) − L(t), increases exponentially due to growth of the cell phase. We note here that our model predicts axially asymmetric tissue growth: equations (3.13) show that the upstream interface is advected at the speed of the imposed flow, whilst advection of the downstream interface is augmented by tissue growth. This growth asymmetry is evident for both static (in which L(t) = L(0)) and dynamic culture conditions.

Linear stability analysis
The stability properties of such one-dimensional tissues have been analysed by a number of authors, especially in the context of solid tumour growth, see e.g. [21,23], in which the effects of material properties and limited nutrient availability on the stability of tumours of constant density were considered and the results were used to characterise the malignancy and fingering instability of such tumours. Here, we consider the effect of an imposed flow on the stability of a growing tissue to disturbances in the transverse direction. Specifically, we consider perturbations to the interfaces which define the construct size and to the cell density within this construct. We perturb the planar interfaces L(t), R(t) as follows: L(y, t) = L 0 (t) + L 1 (y, t) + · · · , (3.14) where 0 < 1/ √ k 1 and we have adopted the following notation: L 0 , R 0 are the planar interfaces defined by (3.13) and L 1 , R 1 are perturbations. Correspondingly, we seek solutions to the governing equations of the form n(x, y, t) = n 0 (x, t) + n 1 (x, y, t) + · · · , (3.16) where n 0 is the one-dimensional solution (3.8), and we consider similar expansions for u, v and p. We remark that the subscript notation which previously indicated the terms in the large drag expansion (2.16) has been replaced by subscripts denoting the solutions associated with the planar interfaces and perturbations. This calculation follows the methodology presented in [21,23], and so much of the detail is omitted.
Returning to the two-dimensional system given by equations (2.18), (2.20) and (2.22), we find that the perturbations to the cell volume fraction, pressure and velocity satisfy Similar conditions apply at x = R 0 . The perturbations to the interfaces are governed by the following kinematic conditions: 24) which are applied at x = L 0 (t) and x = R 0 (t), respectively. To summarise, we now determine the stability of the one-dimensional tissue defined by equations (3.8), (3.12) and (3.13) from the solution of the system (3.17)-(3.24). We proceed by assuming that the perturbations to the interfacial positions are separable, of the form L 1 (y, t) = l 1 (t) cos(λy), R 1 (y, t) = r 1 (t) cos(λy), (3.25) for arbitrary wave number λ, and seek solutions for n 1 of a like form n 1 (x, y, t) = n(t) cos(λy) L 0 6 x 6 R 0 , 0 o t h e r w i s e . (3.26) In view of the boundary condition (3.21), we find λ = qπ/h, for integer q. Using (3.18), we deduce in which G is arbitrary (we set G = 1 in the following without loss of generality). We now introduce a subscript i to denote the region in which the solution is valid; the regions i = 1, 2, 3 correspond to 0 6 x < L 0 , L 0 6 x 6 R 0 , x > R 0 , respectively. Following [21], we write u 1i = f i (x, t) cos(λy) and from equations (3.17) we obtain After some algebra, the twelve functions, A i (t)-D i (t) may be specified in terms of the planar interfaces L 0 (t), R 0 (t), the interface perturbation amplitudes l 1 (t), r 1 (t) and the cell distributions n(t), n(t) (details omitted for concision). Equations (3.24) then yield where L 0 , R 0 are defined by (3.13) and A 1 and C 3 are given by Equations (3.30) are solved numerically using a MATLAB initial value problem solver, which employs an explicit Runge-Kutta formula to compute the solution to non-stiff problems. Figure 2 shows how the perturbations evolve over time for different values of the constant growth rate, k m . Figure 2 shows that there are marked differences in the behaviour of the perturbations to each interface: the amplitude of the upstream perturbation (l 1 ) decreases monotonically with time, passing through zero; conversely, the amplitude of the downstream perturbation (r 1 ) increases monotonically. The effect of this phenomenon on the behaviour of the interfaces x = L(y, t), R(y, t) for the q = 2 mode is illustrated in Figure 3, showing how the reversal of sign of the upstream perturbation amplitude corresponds to a dramatic difference in the behaviour of the up-and downstream interfaces. This behaviour is  due to the perturbation to the cell volume fraction, n 1 = n cos(λy) whose transverse variation results in invagination in the sparse regions and protrusion in the dense regions, as depicted in Figure 3. In the absence of transverse perturbations to the cell volume fraction (n 1 = 0), the tissue density is spatially uniform and the perturbations to each interface do not exhibit this behaviour, both l 1 and r 1 increasing monotonically at exactly the same rate (as reported in [21]); we note that in this case the growth rate of the upand downstream perturbation amplitudes is lower than that observed for the downstream interface in the regime for which n 1 0 (details omitted). Inspection of equations (3.31) and (3.32) indicates that the influence of n 1 on the stability of the interfaces diminishes with increasing wave number, λ(q): as λ increases, the growth of the perturbation amplitudes l 1 , r 1 tends to that observed in the absence of perturbations to the cell volume fraction (n 1 = 0). Furthermore, equation (3.27) shows that for large time (t r −1 ), n 1 decays to zero (in order that we remain within the linear regime we require e rt 1); indeed, for large time, equations (3.30) reduce to  and the one-dimensional solution for which n(x, t) = n(t), L 0 (t) 6 x 6 R 0 (t) is therefore unstable to small transverse perturbations if the net growth rate is positive (we remark that both the amplitude of the perturbations and the construct width increase exponentially with time). It is interesting to note that the stability of the interfaces x = L 0 , R 0 is largely unaffected by the presence of the ambient flow, which only enters the analysis through the advection of the interfaces (see equation (3.13)). Qualitatively similar results are obtained in the zero-flow regime (results omitted).
In the preceding analysis, we have considered the growth and stability of a tissue defined by two planar interfaces within a perfusion bioreactor. For convenience, the bioreactor was modelled as a channel of infinite length; the influence of considering a finite-length bioreactor on the stability properties of this tissue is determined by performing a linear stability analysis in an identical fashion to that given above, with the conditions at infinity now imposed on the truncated boundary, x = X (details omitted). As discussed previously, the lengthscale L * is defined in terms of the growth timescale; the domain length, X, is the number of multiples of this. Equations (3.24) are integrated numerically to yield the behaviour of the perturbations; we find that the up-and downstream planar interfaces are again unstable to transverse perturbations (provided that the net growth rate is positive); the growth rates of the perturbations, l 1 and r 1 match those predicted by the linear stability analysis on an infinite domain provided that the interfaces do not approach x = X (as reported in [21,23] and illustrated in Figure 4).

Numerical solution: D 0
We now present numerical solutions of the one-dimensional equations (3.1) and (3.2a, b) on the truncated domain 0 6 x 6 1 subject to (3.4) and (3.5); the conditions specified as x → ∞ are now imposed at x = 1. We emphasise that, in contrast to the previous section, this system represents a diffuse tissue construct. To illustrate the behaviour of our one-dimensional model, we consider a small cell population intially situated near the upstream end of the bioreactor as follows:  Solutions are calculated with the NAG routine D03PCF which uses a backward differentiation method to solve the system of ordinary differential equations that emerge when the original PDEs are spatially discretised using finite differences. The evolution of n and u is shown in Figure 5. Results for the pressure, p are omitted since it is directly proportional to n (see equation (3.3)). To validate these numerical simulations, in Figure 6, we compare the positions of the up-and downstream construct boundaries predicted by the sharp interface analysis (equations (3.13)) and the numericallycalculated results for a diffuse construct. The positions of the boundaries of the diffuse construct are taken to be the up-and downstream half-maximal values of n; in equations (3.13), n(0) is approximated by the maximum value of n(x, 0). In Figure 7, we compare the evolution of the numerically-computed maximum cell volume fraction and the logistic growth predicted by (3.12).
Inspection of Figures 5-7 shows that the analytical solutions (3.9)-(3.13), corresponding to a growing tissue construct of uniform density defined by two sharp interfaces, capture much of the qualitative behaviour of the full one-dimensional model (3.1), (3.2a, b), (3.4) and (3.5). Figure 5 shows how the cell population is advected downstream by the flow: the diffuse upstream interface moves with constant speed, while the downstream interface is advected at a rate which increases with increasing cell density, and the construct domain elongates accordingly. Figure 6 shows that the positions of these diffuse interfaces are in good agreement with the sharp interfaces defined by equations (3.13). We highlight that, as indicated by Figure 5(a), as t increases towards t = 0.48, the construct approaches the downstream domain boundary and meaningful comparison cannot be made. In addition, Figure 5(a) indicates that the evolution of the peak cell density approximates the logistic growth predicted by (3.12) which is sigmoidal for n(x, 0) < K/2 [47] (where K = 1− k d /k m is the carrying-capacity): following the initial fast growth phase, the peak cell density increases more slowly, tending towards K. In Figure 7 we show that the numericallycalculated peak cell density is in good agreement with that predicted by equation (3.12). We remark that the imposed flow advects the construct to the end of the domain before the steady-state value n = K = 0.99 can be attained. The numerically-computed velocity profile, shown in Figure 5(b), is constant prior to (and after) the tissue, increasing approximately linearly within, with gradient k m n(x, t) (see equation (3.10)) as demanded by the continuity equation (2.18); the downstream velocity increases approximately exponentially with time, as predicted by equations (3.11)-(3.13) (details omitted).
In this section, we have studied the effect of an ambient flow on a one-dimensional tissue construct whose rates of growth and death remain constant. Analytic solutions constructed in the limit for which the interfaces between the growing cell phase and the surrounding culture medium are sharp, were shown to be in good qualitative agreement with numerical simulations for a diffuse cell population. In each case, the effect of the flow is to advect the cells downstream. We find that the stability of the sharp interfaces to transverse perturbations is largely unaffected by the imposed flow; however, the early-time behaviour of these interfaces is dramatically altered by perturbations to the cell volume fraction and we observe markedly different behaviour to that reported in [21].
In the context of the bioreactor system described in §1, this model predicts that the cells and ECM will be advected through the bioreactor at the speed of the imposed perfusion. A long bioreactor or a low flow rate is therefore required to prevent tissue from being flushed out of the bioreactor before tissue growth can be achieved. This prediction is due to the simplifying limit of large interphase drag, in which each phase moves with common velocity. Spatial discretisation is performed using central finite differences; an upwind scheme is used for the convective terms in equation (2.20). Numerical accuracy is maintained via spatial and temporal mesh refinement. For brevity, the full discretisation is omitted here.
To illustrate the behaviour of the model, we focus on how an initially y-independent cell population (corresponding to a uniform cell seeding across the scaffold width) subject to constant cell proliferation and death evolves in response to the imposed flow. Other initial conditions relevant to cell growth on porous scaffolds (such as the more uniform distributions achieved via dynamic seeding on a cortical shaker [64] or cells seeded around the scaffold's periphery) may easily be investigated and will form part of a subsequent study. We choose a similar initial condition to (3.34) as follows: (4.1) Plots of the cell volume fraction, pressure and axial and transverse velocity components are shown in Figures 8-12. Figures 8 and 9 show how the initially sparse, y-independent construct given by (4.1) grows and spreads through the channel. The evolution of the cell volume fraction at the channel centreline is shown in Figure 9 indicating the advection of the tissue construct along the bioreactor. The cells are advected along the channel by the axial velocity, u, which is parabolic in y. This pressure-induced parabolic flow increases advection in the channel centre where flow speed is maximal, introducing significant transverse variation. A parabolic flow profile is obtained since the model does not account for the influence of the porous scaffold on the fluid dynamics of the bioreactor which leads to a "plug flow"   Figure 8. Parameter values as given in Figure 8. Figure 11. A surface plot of the axial velocity profile corresponding to the tissue construct in Figure 8. Parameter values as given in Figure 8.
profile within the scaffold. Our model therefore overestimates the amount of transverse variation induced in the cell population by the flow. The influence of the scaffold will be accounted for in a subsequent study. We remark that the solution for the axial velocity, u, is unique up to the addition of an arbitrary one-dimensional solution,û; we have chosen u = 0 so that no-slip is assured prior to the tissue; within and downstream from it the slip velocity is significant. The advection/growth behaviour of the tissue construct warrants further discussion. Comparison of the profiles shown in Figure 8 and the initial condition (4.1) shows that  Figure 8. Parameter values as given in Figure 8. advection of the upstream periphery of the diffuse construct is greatest at the channel centreline (y = h/2) where the parabolic flow is maximal. Near the channel walls, where the flow speed is low, limited upstream movement is observed due to the presence of diffusion (for a larger diffusion coefficient than that employed in §3.2, this behaviour is, of course, also observed in the case of a one-dimensional diffuse tissue). Corresponding transverse variation is induced in the downstream periphery; however, it displays significantly greater axial advection due to the greater flow speed there. In the absence of an imposed flow (P u = 0 = P d ), the resulting cell population remains independent of the transverse coordinate. However, the population does not grow symmetrically about the midpoint of the initial distribution: the upstream periphery remains almost stationary while the downstream periphery spreads downstream (results omitted). In each case, this axially asymmetric advection behaviour is predicted by the one-dimensional, sharp interface analysis which indicates that the upstream interface moves at the speed of the imposed flow, the advection of the downstream interface being augmented by tissue growth (see equation (3.13)). Figure 10 shows how the pressure distribution is affected by the presence of the cells. Up-and downstream from the tissue the pressure field decreases linearly with x; an increase in pressure is observed as the fluid flows through the area in which cells are present. As in the one-dimensional case ( §3), the deviation from the linear pressure profile mirrors the cell distribution shown in Figure 8. Similarly, the axial velocity profile shown in Figure 11 is greatly affected by the cells' presence. Upstream, the cells do not influence the flow and u remains x-invariant; however, where n 0 the axial flow speed increases with x in an approximately linear fashion, the gradient increasing with n as required by (2.18). Again, this behaviour was indicated by the one-dimensional analysis; see equations (3.9)-(3.11) and Figure 5(b). We note that the choices P u = 1, P d = 0 set the velocity scale, U * w , used in the non-dimensionalisation (2.9) to be U * w = P * u L * /μ * w . In the absence of cells, the maximum axial velocity is given by u * = P * u h * 2 /(8μ * w L * ); for the parameter choice given, we therefore expect the maximum upstream dimensionless velocity to be u = 1/8. Figure 11 confirms this.
The transverse component of velocity (v) is initially small due to the one-dimensional initial cell population. However, as n increases, transverse variation is introduced within the tissue construct and v increases, achieving maxima on the periphery of the construct (see Figures 8 and 12). We note that the magnitude of the axial velocity, u, is much greater than the transverse velocity, v. This is because we impose an axial pressure gradient to drive the flow and, due to the no-penetration boundary conditions, the resistance to transverse motion in the channel is greater.

Mechanotransduction
The dominant mechanical stimuli relevant to specific tissue engineering applications have not all been elucidated. By extending our model to consider the effect of coupling the growth of the cell population to the local environment, we can determine the characteristic growth pattern associated with specific mechanical stimuli; in tandem with experimental data, this will allow optimisation of culture conditions to enhance yield. For illustrative purposes, we consider the response of the cells to the following stimuli: contact inhibition caused by cell-cell interactions, the effect of stress caused by increases in local cell density and the influence of the external fluid dynamics. The coupling is achieved by modifying the net mass production term, S n . We remark that since the "cell" phase comprises cells and ECM, modifying the growth and death rates (k m , k d ) in response to local environmental factors enables crude modelling of a phenotypic switch due to mechanical stimuli from, for instance, a proliferative phase to an ECM-producing phase.
The relevance of the following work hinges upon the appropriate choice of S n . We restrict attention to the following cell density and pressure-dependent responses (1) Cell density-dependent response: S n = [k m (n) − k d ] n = κ(n)n; (2) pressure-dependent response: in which k m , k d are the rates of growth and death, respectively. For ease, the death rate k d is kept constant; the coupling of cell behaviour and the mechanical environment is captured entirely through the growth rate, k m . We now motivate these choices. Although not explicitly modelled, the choice k m (n) enables us to capture the effect of contact inhibition [10] and tissue growth-induced stress [24,57] on cell behaviour. Furthermore, since the pressure, p, represents a "lumped pressure" (p = p w0 + n 0 Σ n0 ; see §2.1), the choice k m (p) incorporates these considerations as well as the cells' response to the culture medium pressure. The response of cells to culture medium pressure is well documented, especially with respect to bone tissue growth; for example, many authors have shown that bone cells respond to intermittent hydrostatic compression with diminished bone resorption and enhanced bone formation ( [34] and references therein), increased adhesion [27] and increased osteopontin (a protein implicated in the bone remodelling process) expression [48]. Excessively high hydrostatic pressure (> 200 kPa) has been shown to exert an inhibitory effect on bone-specific gene expression [56].
Considering first the density-dependent behaviour described above, we identify three distinct stages in the behaviour of the cell population: (i) a proliferative phase, (ii) an ECM-producing phase and (iii) an apoptotic phase. At low density, the cells proliferate at a rate, k m (n) = k 1n ; at intermediate density, due to the additional production of ECM, the growth rate of the cell phase is modified to a new value, k m (n) = k 2n (we assume k 2n > k 1n ); finally, when the local density is too high, the cells enter an apoptotic phase (k m (n) = 0) with death rate k d (we note that the "death rate", k d includes ECM degradation as well as cell death). The threshold cell densities that separate these three types of behaviour are denoted n 1 and n 2 .
Similarly, in the pressure-dependent regime, we assume that at intermediate levels of pressure, the cells exhibit enhanced proliferation and ECM deposition (k m (p) = k 2p ); at low pressure, the cells enter a quiescent state in which proliferation and ECM deposition is greatly reduced (k m (p) = k 1p < k 2p ); at excessively high levels, ECM deposition ceases and the cells become apoptotic (k m (p) = 0). The corresponding thresholds are denoted p 1 , p 2 .
We now present numerical solutions of the two-dimensional equations (2.20)-(2.22) subject to the boundary conditions (2.23)-(2.25) in which we employ the above choices for S n ; initial conditions are given by (4.1). We assume that the rates of growth and death of the cell phase (k 1n , k 1p , k 2n , k 2p , k d ) are constant and represent the proliferative responses described above with a smoothed step function, as defined below: In (4.2), ϕ represents the stimulus in question, with corresponding threshold values ϕ 1 , ϕ 2 and the parameter g dictates the level of smoothing; we note that smoothing is necessary in order to obtain numerical solution. Figure 13 shows a sketch of the function k(ϕ), highlighting the progression from one phase to the next. The effect of these choices of mass production term on the morphology of the resulting tissue is shown in Figures 14 and 15. Figure 14 illustrates that when the cells' response is density-dependent, the growth of the cell phase is arrested at n = n 2 due to the smoothed progression from proliferation and ECM-production (κ(n) = k 2n − k d ) to apoptosis (κ(n) = −k d ). We note that despite the presence of apoptosis in this formulation, regression back to the proliferative phase ensures that the cell density does not fall below n = n 2 .   Figure 15 shows the response of the cell phase to pressure-dependent growth. Rather than being arrested at a threshold density the cells become apoptotic where the pressure is high (near the upstream diffuse interface and in regions of high cell density; see Figure 10) and proliferation is reduced near x = 1 (where the pressure is low); between these regions, growth is enhanced. The result of this spatial variation in proliferative rate is a tissue construct which grows preferentially downstream in the regions of intermediate pressure.
Comparison of the cell phase distribution in each of the above growth regimes with that obtained in the case of constant growth and death rates (Figures 8, 14 and 15) shows that the composition of the tissue construct is dramatically affected by coupling the growth response of the cells to their environment. When cell proliferation and ECM deposition are density-dependent, a uniform tissue construct is obtained; in the pressure-dependent case, the predicted tissue construct composition is far less uniform. It is interesting to note that in the absence of an imposed flow, the pressure field is directly proportional to the cell phase distribution (see equation (3.3)) and the cell density-and pressure-dependent responses are identical. Inspection of the morphology of tissue constructs produced in static and dynamic culture, together with the characteristic growth patterns predicted by this model, therefore provides a simple means to identify the dominant regulatory growth stimulus in a cell population.

Discussion
In this paper, we have studied tissue growth in a perfusion bioreactor which is represented by a two-dimensional channel containing a two-phase mixture of interacting viscous fluids. The cells and ECM were modelled as a single phase; the second phase represented the culture medium. Guided by parameter estimation, we employed the limit of large interphase drag, in which we may describe each phase as being subject to a common velocity and pressure field.
We have considered the effect of a dynamic flow environment on tissue growth processes. Analytic predictions were obtained in the limit of a one-dimensional growing tissue defined by two sharp interfaces, between which the cell density remains spatially invariant. The cell phase displayed logistic growth, the interfaces were advected by the ambient flow and the tissue width increased exponentially. The stability of this tissue to periodic transverse perturbations was investigated, and the interfaces were found to be unstable; the long-time stability being regulated only by the net growth rate. In the presence of a corresponding perturbation to the cell volume fraction, the perturbations to the upstream interface reversed sign due to the variation in construct density. This effect diminished with increasing perturbation wave number and decayed to zero for large time.
Using numerical simulations in one and two-dimensions, the behaviour of a diffuse tissue construct under an ambient flow was calculated for constant growth and death rates and the advection behaviour predicted by the sharp interface analysis was observed, indicating that this asymptotic limit captures much of the qualitative behaviour of the full system. In the two-dimensional model, transverse variation in the tissue construct density was induced by the parabolic flow of culture medium; small transverse flows were induced at the construct periphery.
Our analysis indicates that cells and ECM are advected through the bioreactor at the same speed as the imposed flow, implying that a long bioreactor and/or a low rate of perfusion is required in order to prevent the tissue from being flushed from the bioreactor before tissue growth can occur. This is a consequence of the simplifying limit of large interphase drag employed in this paper which demands that the cell and culture medium phases are subject to a common velocity field.
We further extended this model formulation to account for complex coupling between the cells' proliferative response and their local environment. This was achieved by replacing the constant growth and death rates (k m , k d ) with appropriate functional forms. Specifically, motivated by a range of studies, we considered the response of a cell population to the local density and pressure. Simulations were presented showing that the growth of the cell population is profoundly altered by these effects, dramatically changing the composition of the construct. These simulations clearly demonstrate the importance of considering the effect of mechanotransduction mechanisms within tissue growth models. Furthermore, our model suggests that in static culture, regulation of proliferative behaviour by cell density and culture medium pressure results in indistinguishable tissue constructs. In principle, on provision of appropriate experimental data, this conclusion provides a simple mechanism for the identification of the dominant regulatory mechanism in a given cell population simply by observing the tissue construct morphology resulting from culture in static and dynamic conditions. However, we note that we have not considered nutrient-limited growth which is expected to become significant in the absence of perfusion (indeed, after many days in culture, delivery of nutrients to downstream sections of the scaffold may be problematic even under perfusion, especially in scaffolds of relevant clinical thickness) and may affect the robustness of our predictions [38]. The influence of nutrient-limited growth will be a focus of future work.
The proliferation/ECM deposition functions were chosen to reflect qualitatively a simple mechanotransduction mechanism. Our model admits more complex functional forms and dependence on combinations of the field variables; physiologically, it is expected that these effects are integrated in a complex way to produce the cells' overall response. The simple forms employed here allow clearer illustration of the importance of mechanotransductionaffected growth within a tissue growth modelling framework; however, the mathematical formulation and numerical scheme developed is highly versatile, permitting the study of more complex functional forms and an investigation of the interplay between many competing growth stimuli as appropriate experimental data become available.
The model employed here is applicable to tissue growth processes in which the solid characteristics of the system are unimportant. In a subsequent paper, we will present a three-phase formulation in which we distinguish between the ECM and cell phases and consider mechanotransduction mechanisms in more detail. This more complex model will allow further predictions regarding the effect of perfusion the composition and mechanical integrity of the construct. Furthermore, the effect of cell-cell and cell-scaffold interactions on cell movement will be studied.