Multiphase modelling of the effect of fluid shear stress on cell yield and distribution in a hollow fibre membrane bioreactor

We present a simplified two-dimensional model of fluid flow, nutrient transport and cell distribution in a hollow fibre membrane bioreactor, with the aim of exploring how fluid flow can be used to control the distribution and yield of a cell population which is sensitive to both fluid shear stress and nutrient concentration. The cells are seeded in a scaffold in a layer on top of the hollow fibre, only partially occupying the extracapillary space. Above this layer is a region of free-flowing fluid which we refer to as the upper fluid layer. The flow in the lumen and upper fluid layer is described by the Stokes equations, whilst the flow in the porous fibre membrane is assumed to follow Darcy’s law. Porous mixture theory is used to model the dynamics of and interactions between the cells, scaffold and fluid in the cell–scaffold construct. The concentration of a limiting nutrient (e.g. oxygen) is governed by an advection–reaction–diffusion equation in each region. Through exploitation of the small aspect ratio of each region and asymptotic analysis, we derive a coupled system of partial differential equations for the cell volume fraction and nutrient concentration. We use this model to investigate the effect of mechanotransduction on the distribution and yield of the cell population, by considering cases in which cell proliferation is either enhanced or limited by fluid shear stress and by varying experimentally controllable parameters such as flow rate and cell–scaffold construct thickness.


Introduction
The need for a reliable and sufficient source of replacement tissue and organs is constantly increasing, due to our ageing population and consistent lack of donors. The potential for tissue engineering to meet this demand is great, but there are still many barriers to be overcome before the field can provide alternative treatments on a clinical scale. For each tissue type, a different experimental technique must be developed to ensure that the engineered substitute is viable and can perform the same essential functions as the original tissue. This often means being able to mimic the in vivo environment of a particular tissue and has resulted in a vast array of protocols, bioreactors, culture conditions, and scaffold materials (see e.g. Martin et al. 2004;Martin and Vermette 2005;Pörtner et al. 2005;Stock and Vacanti 2001). For each set-up and cell type, optimal operating conditions must therefore be determined, a process which can be extremely time-consuming and expensive to resolve purely experimentally.
Mathematical modelling of bioreactor systems can be of great benefit for several reasons. It is reproducible and efficient, allowing the large numbers of parameters to be investigated relatively quickly and cheaply. Moreover, it can give insight into the combined effects of the physical processes involved in a particular set-up, or even focus on one process in isolation-this is hard to achieve in an experimental set-up. An example of one such process, and the focus of this work, is mechanotransduction: the mechanism by which forces are converted into biochemical signals and integrated into a cellular response. Below, we review a selection of mathematical models that specifically incorporate the effects of mechan-otransduction. For more general detail on mathematical models in tissue engineering, we refer to the recent review by O'Dea et al. (2013a).
Previous work on multiphase modelling of mechanotransduction in tissue engineering includes a series of papers by O'Dea et al. (2008O'Dea et al. ( , 2010O'Dea et al. ( , 2013b. In O'Dea et al. (2008), the authors develop a two-phase model of tissue growth in a perfusion bioreactor. They include dependence of cell proliferation, extracellular matrix (ECM) deposition and cell death on the local cell density and fluid pressure. Results indicate that these effects can dramatically alter the composition of the resulting cell construct. This model was extended to include an additional porous scaffold phase in O' Dea et al. (2010), in which the system was simplified analytically by taking the long-wavelength limit appropriate for a small channel aspect ratio. The effect of fluid shear stress, as well as cell density and fluid pressure, was considered. Findings supported the earlier conclusion that mechanotransduction can significantly affect cell distribution. The same set-up was modelled and solved numerically by Osborne et al. (2010), highlighting situations in which a numerical approach is necessary, and when the analytical limits taken in O' Dea et al. (2010) are relevant. A further extension of this work is given in O' Dea et al. (2013b) in which the degradation of a scaffold is included, and whose spatially non-constant porosity is informed by experimental data.
As well as cell proliferation and death, mechanotransduction can affect cell differentiation. This is the focus of a paper by Byrne et al. (2007), in which finite element analysis is used to model a poroelastic tissue-scaffold construct for bone regeneration, with the objective of determining optimal scaffold porosity under different mechanical loading conditions. The dependence of stem cell differentiation on fluid velocity and shear strain is considered by employing the mechanoregulation algorithm developed by Prendergast et al. (1997). Depending on the stimulus level, differentiation into fibroblasts, chondrocytes or osteoblasts can occur. Aiming to identify possible treatment options in fracture healing, Lacroix and Prendergast (2002) develop a three-dimensional model of a human tibia fracture, again using the mechanoregulation algorithm proposed in Prendergast et al. (1997). Two different compressive loading magnitudes were simulated, and it was shown that only the lower load led to successful healing.
A finite element approach was also taken by Driessen et al. (2003), who modelled the effect of the different mechanical loading regimes resulting from closed and open configurations on collagen fibre content and orientation in an aortic valve. The mechanical properties of the tissue construct were shown to depend on the type of loading, as the fibres aligned with the principal strain directions, and the results were validated against experimental data where possible. Cardiac tissue was also the focus of the work in Latimer et al. (2003), who determined the stress and strain distribution in two-and three-dimensional tissues. Analytical solutions were found and used to demonstrate the strain in and around an ischemic region, and have the potential to be further used to predict mechanotransduction effects on the tissue leading to arrhythmias, or following hypoxia. Finally, the importance of scaffold design in bone tissue engineering motivated the study by Sanz-Herrera et al. (2009). The need for an implantable construct that has the required macroscale properties, but which is affected by biophysical phenomena on the microscale, led to a multiscale model which incorporated a realistic scaffold microstructure. Bone regeneration was assumed to be dependent on mechanical stimuli, and properties for the macroscale model such as the scaffold stress and strain and cell diffusion were obtained through the analysis of the underlying microstructure, which in turn evolved as a result of macroscale cell migration.
In this paper, we consider various experimentally relevant case studies, motivated by specific cell types, to investigate the effect of fluid shear stress on the proliferation rate, and thus distribution, of cells in a hollow fibre membrane bioreactor (HFMB) (see Fig. 1). A HFMB consists of a cylindrical, glass module with a port at each end of the extracapillary space (ECS) and a porous hollow fibre inserted through the centre. A number of different cell seeding and flow regimes can be employed; for this work, we consider the set-up in which the cells are seeded in a scaffold in the ECS between the hollow fibre and the outer glass wall. Culture medium is pumped in via the lumen inlet and upstream ECS port. The pressure is set at the downstream lumen end via a clamp, and the downstream ECS port is left open to the atmosphere. Depending on the flow and pressure conditions, the fluid may pass through the hollow fibre walls and then flow out of either the lumen outlet or downstream ECS port.
We build on the model developed in Pearson et al. (2013) which describes fluid flow, solute concentration and cell distribution in a HFMB. In Pearson et al. (2013), the entire ECS was filled by the cell-scaffold construct, through which media could flow. Fluid velocities in the membrane and ECS were assumed to be an order of ε smaller than those in the lumen as a result of the resistance to flow through these porous regions, where ε 1 is the lumen aspect ratio. Here, we consider the alternative experimental set-up in which the cell-scaffold construct only fills part of the ECS (that directly above the porous membrane, defined as the cell layer), and media can flow freely throughout the remainder. The width of this cell layer can be controlled experimentally and will be varied as part of our investigation. Fluid is pumped into the system via both the lumen inlet and upstream ECS port; as a result, we expect the fluid velocity in this upper fluid layer to be of the same order of magnitude as the lumen. As well as enabling the relevant solutes to be delivered to, or removed  Pearson et al. (2013). Bottom cross section of the boxed region (not to scale), where the lower half (not shown) has been excluded based on symmetry. This depicts the idealised two-dimensional modelling region with the x axis running along the lumen centreline. The solid black arrows show the direction and location of the fluid fluxes into the system, and the star denotes the origin (x, y) = (0, 0) from, the cells more quickly, we expect the cells to be exposed to a higher and potentially significant level of fluid shear stress as the flow through the cell layer will be enhanced by the upper fluid flow above. The impact of this set-up remains an open question experimentally: depending on the cell type, fluid shear stress can either be beneficial or detrimental to cell growth and survival. In this work, we therefore consider mechanotransduction effects, which were not included in Pearson et al. (2013). Specifically, we explore the effect of fluid shear stress on cell yield and distribution, via changes in the fluid flux into the upper fluid layer. In addition, the fluid flux affects the concentration of a nutrient (taken to be oxygen) which can limit the cell proliferation rate. We determine the range of possible behaviours from this set-up and find optimal flow rates for which the fluid shear stress has an advantageous effect: that is, it enables a more spatially uniform cell population and/or a higher cell yield to be obtained.

Paper outline
We begin in Sect. 2 by describing the simplified modelling domain before introducing the governing equations and boundary conditions for each region in Sects. 2.1 and 2.2. In Sect. 3.1, we introduce relevant parameter values, motivating the non-dimensionalisation of the system (details of which are given in "Appendix"). In Sects. 3.2-3.4, we exploit the small aspect ratio of the lumen to reduce the system to two coupled partial differential equations for the cell volume fraction and solute concentration. Results for the cases in which cell proliferation is either enhanced or limited by fluid shear stress are presented in Sects. 4.1 and 4.2, and the sensitivity of our results to the cell layer width is investigated in Sect. 4.3. Finally, key findings and conclusions are discussed in Sect. 5.

Model description
We describe the system using two-dimensional Cartesian coordinates, for simplicity and to enable analytical progress. Although results would be quantitatively different in an axisymmetric or three-dimensional set-up, we would not expect them to change qualitatively and hope to verify this in future work. We define the axial dimension of the modelling domain by 0 < x < L. The lumen and membrane are then respectively given by 0 < y < h 1 and h 1 < y < h 1 + h 2 in the transverse direction. The cell layer is of comparable thickness to the membrane, occupying h 1 + h 2 < y < h 1 + h 2 + h 3 , and the upper fluid occupies H − h 4 < y < H , where H = h 1 + h 2 + h 3 + h 4 (see bottom half of Fig. 1). As in Pearson et al. (2013), we consider the cell layer to be a multiphase region consisting of the cells, culture medium (modelled as water) and a rigid, inert scaffold, closely following the formulation from Lemon et al. (2006). Water variables in the lumen, membrane, cell layer and upper fluid layer are denoted by subscripts l, m, w and f, respectively, and cell phase variables by subscript n, with the velocities given by u i = (u i , v i ), the water pressures by p i , and the solute concentrations per unit volume of water by c i (i = l, m, w, f). In the cell layer, we track the dynamics of the cell, water and scaffold phases through their respective volume fractions θ n , θ w and θ s , where θ s is constant in space and time due to the assumption of a rigid, inert scaffold.

Governing equations
The dimensional governing equations in each region take the form of conservation of mass and momentum for the water and cell phases, and conservation of mass for the solute. In the lumen, porous membrane and cell layer are the same as for the three-region system in Pearson et al. (2013), and the dynamics in the additional upper fluid layer are governed by the same equations as the lumen. We take the membrane porosity φ to be constant and work with reduced pressures throughout since gravitational effects are not negligible at the flow rates considered (see Pearson et al. 2013 for details). In the lumen (0 < y < h 1 ), the relevant water equations are those of Stokes flow, along with an advection-diffusion equation for the solute: where t represents time, μ w the water viscosity (assumed constant) and D the diffusion coefficient for the solute in water (also assumed constant). In the porous membrane (h 1 < y < h 1 + h 2 ), we use Darcy's law for flow in porous media, where k is the (constant) membrane permeability. In the cell layer (h 1 + h 2 < y < H − h 4 ), the no-voids condition is given by while conservation of mass and momentum for the cell phase is given by Here, J n represents the net cell production rate, p n is the (reduced) cell pressure, ψ ns is the interphase pressure due to tractions between the cells and scaffold, γ ns is the (constant) cell-scaffold drag coefficient, and τ n is the deviatoric stress tensor for the cell phase, given by where superscript T denotes transpose, I is the identity tensor and μ n is the effective viscosity of the cell phase (assumed to be much greater than μ w ). Here, we have neglected tractions and drag between the cells and water, assuming that they are much smaller than those between the cells and scaffold. This simplifying assumption is motivated by the expectation that the cells will strongly adhere to the scaffold. In the equations for the water phase, below we similarly neglect tractions between the cells and water. The corresponding equations for the water phase are where we have assumed that the net water production rate is −J n so that mass is conserved. Conservation of solute mass in this region yields the following: where R is a reaction term which accounts for solute uptake by the cells. Finally, in the upper fluid layer (H − h 4 < y < H ), the governing equations are Constitutive forms need to be prescribed for a number of the terms above, namely J n , R (which we present and discuss in Sect. 4) and p n and ψ ns . For now, we assume only that J n and R are functions of the cell volume fraction θ n , cell layer solute concentration c w and (in the case of J n ) the fluid shear stress in the cell layer. As will be seen in Sect. 3.2, the fluid flow in the cell layer is independent of y, and thus, we take the fluid shear stress to be proportional to the fluid pressure gradient ∂ p w /∂ x. The cell pressure is assumed to be equal to the water pressure plus an extra component Π which accounts for cell-cell interactions, for constants ν and δ. The form for the extra component Π is motivated by O'Dea et al. (2010), and the first term models aggregation of cells at low densities, whilst the second term represents contact inhibition. As in Pearson et al. (2013), ψ ns is assumed to be a negative constant representing the cells' affinity for the scaffold: (2.11)

Boundary conditions
On y = 0 symmetry requires that on the lumen/membrane interface y = h 1 we impose no slip of fluid, and continuity of fluid flux, of normal stress, of solute concentration and of solute flux, viz.
(2.13) and on the membrane/cell layer interface y = h 1 + h 2 we impose no flux and no slip of cells, and continuity of fluid flux, of normal stress, of solute concentration and of solute flux, viz. (2.14) On the cell layer/upper fluid interface y = H −h 4 , we impose no flux and no shear stress on the cell phase, continuity of flux and of normal stress on the water phase, no slip of fluid and continuity of concentration and of flux of solute: v n =0, n e · σ n · t e = 0, (2.15) We note that, in the above, the most appropriate choice of stress condition on the cell phase is unclear, and we have chosen to apply no shear stress as an example first case. Equally, a nonzero shear stress condition could be chosen, and the analysis follows through in a similar manner. However, this would require prescription of another constitutive term for the partition of the shear stress, and hence, we have chosen not to consider this here for simplicity. Finally, on the top of the bioreactor y = H , we impose no slip and no flux of fluid, and no flux of solute: We impose boundary conditions at the up-and down-stream ends of the bioreactor once the system has been reduced in dimension, as we shall now describe.

Model reduction
In this section, we discuss relevant dimensional parameter values, which motivate our choice of non-dimensionalisation (details are given in "Appendix"). The resulting dimensionless parameter values, and asymptotic expansion of the system variables in powers of the lumen aspect ratio, allow the system of governing equations and boundary conditions to be reduced significantly to four coupled PDEs. Imposing upand down-stream boundary conditions closes the problem and eliminates two of these equations, leaving a system of two coupled PDEs for the cell volume fraction and solute concentration at leading order in the lumen aspect ratio. The remaining leading order unknowns can be determined via analytical expressions.

Parameter values
Typical dimensional parameter values are taken from the literature or guided by our experimental collaborators where possible. We present these in Table 1 along with corresponding dimensionless parameters in Table 2. Dimensional parameter values for which data could not be found are omitted from Table 1; instead, their dimensionless equivalents are defined in Table 2 and values chosen so that our asymptotic analysis retains as many features as possible at leading order. A few choices are worth mentioning here: for a more indepth discussion see Pearson et al. (2013). Firstly, we note that the aspect ratio of each region is small; in particular, we define the lumen aspect ratio by ε = h 1 /L 1 and use this as the relevant parameter for our asymptotic expansion in Sect. 3.2. In addition, the reduced Reynolds number in the lumen and upper fluid layer is ε 2 Re = 2.05 × 10 −6 which justifies neglecting inertial effects in these regions in Typical porous membrane velocity We set the viscosity ratio μ w /μ n = λε for some constant λ = O(1). This is motivated by the fact that we would expect the effective (macroscale) viscosity of the cell phase to be much greater than that of water as a result of the cytoskeletal network, the cell-scaffold affinity and microscale cell-cell interactions. We choose the dimensionless cell-scaffold drag ζ ns and water-scaffold drag ζ ws to be of O(1), so that their effects are retained in the leading order model, but assume that ζ ws < ζ ns (as would be expected physically). We define the velocity scale U * to be a typical horizontal velocity of the water in the porous membrane. We then take the corresponding horizontal velocity scale in the cell layer to be U * , and in both the lumen and upper fluid layer to be μ n U * /μ w , so that the same pressure scales apply throughout the system in our non-dimensionalisation. These choices are also motivated by the fact that we would expect the porous structures to hinder the fluid flow in the membrane and cell layer. In addition, this choice enables substantial analytical progress to be made in the mathematical reduction in Sect. 3.2. We note that other lim-its could be considered, for instance taking the cell and water phases to be of comparable viscosity (as in O'Dea et al. 2010O'Dea et al. , 2013a or the cell phase to be inviscid (as in Lemon et al. 2006).
We take the reduced Péclet number in the lumen and upper fluid layer, ε 2 Pe, to be order ε 2 so that diffusion dominates advection throughout. The concentration scale C * and inlet concentration c in are both set to be a typical oxygen concentration used in cell culture (Shipley and Waters 2012), as oxygen is the solute of interest in Sect. 4. The two-dimensional lumen and upper fluid inlet fluxes Q l,in and Q f,in (which come into the up-and down-stream boundary conditions in Sect. 3.3) have been converted from three-dimensional experimental values by first calculating the corresponding velocity and then multiplying by the length scale in the y direction, εL. The downstream lumen outlet pressure P d and atmospheric pressure p atm are also introduced in Sect. 3.3. The cell proliferation/death rates Γ nw /Γ wn and solute uptake rate Γ R 1 are parameters in the constitutive laws for J n and R, both of which we define in Sect. 4. Based on estimations by our experimental collaborators, we take the (dimensional) cell proliferation rate coefficient Γ nw to correspond to one cell division every 48 h, and we assume that cells live, on average, for 28 days when fixing the (dimensional) cell death rate coefficient Γ wn .

Derivation of the reduced model
Having non-dimensionalised the governing equations and boundary conditions (see "Appendix"), we asymptotically expand each dependent variable in powers of ε, for instance setting u l ∼ u l 0 + εu l 1 + ε 2 u l 2 + · · · and similarly for the remaining velocities, reduced pressures, solute concentrations and volume fractions. In the following, we omit the subscript 0 from leading order variables except where needed for clarity. Equating coefficients of ε 0 in the lumen then yields

2b)
∂ 2 p m ∂ y 2 = 0, (3.2c) in the cell layer and in the upper fluid layer The leading order boundary conditions on the lumen centreline are (3.5c) on the lumen/membrane interface on the membrane/cell layer interface and finally on the bioreactor top Consideration of the above solute equations and boundary conditions shows that at leading order, there is a global concentration that is independent of y: c l = c m = c w = c f = c(x, t). To close the problem and find c(x, t), we must also consider the solute equations at O(ε 2 ). We note that the O(ε) components of the solute equations and boundary conditions in each layer indicate that ∂c i 1 /∂ y = 0(i = l, m, w, f), and so (in the lumen, membrane, cell layer and upper fluid layer, respectively) together with boundary conditions ∂c l,2 ∂ y = 0 on y = 0, (3.11) A substantial amount of analytical progress can be made in the distinguished limit considered above. Equations (3.1c), (3.3g) and (3.4c) respectively tell us that the leading order fluid pressures in the lumen, cell layer and upper fluid layer are independent of y. From (3.3e), and given θ s is constant, we can thus conclude that the cell volume fraction θ n is also a function of x and t only, which allows substantial simplifications to be made. We can obtain the following expressions for the leading order variables in terms of p l (x, t), p f (x, t), θ n (x, t) and c(x, t): where M(x, t) is given by Hence, at leading order, we have Poiseuille flow in the lumen and upper fluid layer and transverse flow only in the porous membrane. In the cell layer, the horizontal fluid flow is Darcylike. The variables v n and v w can be determined from the leading order conservation of mass Eqs. (3.3b, c), respectively.
In addition, we have the following equations for p l , θ n , the global leading order concentration c and p f : where the cell flux Q n is defined by and the solute flux Q c by (3.25) The equations for p l , θ n and p f come from integrating the respective conservation of mass Eqs. (3.1a), (3.3b) and (3.4a) across the lumen, cell layer and upper fluid layer. The remaining equation (3.23c) has been obtained by adding together the O(ε 2 ) solute equations (3.10a-d) after averaging each across the appropriate region. We now impose up-and down-stream boundary conditions to close the problem.

Boundary conditions for the reduced model
We mimic the experimental set-up where fluid is pumped into the bioreactor via both the lumen inlet and (upstream) ECS port, and so impose where Q l,in and Q f,in are prescribed (dimensionless) volume fluxes per unit length in the z direction (perpendicular to x and y) into the lumen and upper fluid layer, respectively, and are assumed constant. We note that in this case, the dimensional Q f,in is of the same order as the dimensional Q l,in given that the velocity scales in the lumen and upper fluid layer are the same. We also prescribe a down-stream lumen pressure and mimic the atmospheric pressure conditions at the downstream ECS ports by setting For the cells, we impose no flux out of the modelling domain (which would be achieved using filters in an experiment): Finally, we prescribe an inlet solute concentration and assume no diffusive flux at the downstream end of the bioreactor: As there is no downstream constraint on the solute concentration experimentally, condition (3.31) is motivated by the fact that we would expect the concentration to be constant in space as it leaves the bioreactor, due to the effects of diffusion. Applying boundary conditions (3.26)-(3.28) allows us to solve (3.23a, d) explicitly for p l and p f , giving and thereby specific expressions for u l , v m , p m , u w , p w and u f from (3.16)-(3.19). In particular, we note the form for the fluid pressure in the cell layer, which is equal to p f (see 3.18b), (3.33)

Summary of reduced model
In summary, we have derived the following coupled system to solve for θ n and c: with the boundary conditions Here M and Φ are as in (3.21) and (3.22), and the appropriate initial condition is the prescription of θ n at t = 0 for 0 < x < 1. We can compare this to the corresponding system from Pearson et al. (2013): it is interesting to see that, despite the apparently more complex set-up here (with four modelling regions instead of three), we can reduce the governing equations to a system of two instead of three coupled PDEs. This is due to the addition of the upper fluid layer region which was not present in Pearson et al. (2013), and in which we can explicitly solve for the reduced water pressure. The reduced water pressure in the cell layer is the same as that in the upper fluid layer, and hence only two unknowns remain, θ n and c.

Numerical results
The reduced system (3.34)-(3.40) was solved numerically using the method of lines, first discretising in x and then performing the time integration using the MATLAB function ode15s. We are interested in steady states of the system, as these arise on timescales comparable to the long culture time of experiments (see Pearson et al. 2013 for a more thorough discussion). These states were found to be independent of the choice of initial condition in our simulations, but for completeness, we note that in the simulations presented in this paper, we set θ n (x, 0) = 0.3, c(x, 0) = 1 for 0 < x < 1. We consider the effect of the fluid shear stress on the cell distribution and yield, in the case where the solute of interest is a nutrient (e.g. oxygen), and with the aim of determining optimal experimental conditions which result in both a spatially uniform cell population and high cell yield.
In each of the following sections, we consider the Michaelis-Menten type reaction term used in Pearson et al. (2013) for nutrient-limited proliferation: We take the cell mass transfer term J n to have the form where S is the (dimensionless) fluid shear stress in the cell layer. If F(S) were equal to unity, then we would recover the form used in the case of the nutrient-driven proliferation discussed in Pearson et al. (2013). In the following analysis, we take F(S) to have different forms depending on the cells' response to shear. As we effectively have Darcy flow of water through the cell layer, we use the following estimate for S motivated by a similar expression in Whittaker et al. (2009): based upon the assumption that there is Poiseuille flow through each of the 'pores' in the cell layer, where the pores are approximated as circular ducts. Given the assumed form for S, and the expression for ∂ p w /∂ x in (3.33), we can see that the fluid shear stress is dependent on both the upper fluid layer flux Q f,in and cell layer width h 3 (via h 4 = H −(1+h 2 +h 3 )), both of which can be controlled experimentally. In what follows we investigate the sensitivity of our results to variations in these two parameters. Moreover, for a given flow rate and cell layer width, we note that variation in S arises purely through variations in θ w . In each of the following case studies, we present results for the steady-state cell volume fraction θ * n (x) and use the mean μ and standard deviation σ of the cell volume fraction θ n to characterise the cell yield and distribution, where (4.5)

Shear stress-enhanced proliferation
First of all, we consider the case when cell proliferation is enhanced by increased levels of fluid shear stress. This has been observed in cells such as osteoblasts (Kapur et al. 2003), chondrocytes (Malaviya and Nerem 2002) and (more controversially) endothelial progenitor cells (Yamamoto et al. 2003), and is thought to be an important factor in ensuring the mechanical integrity of tissue such as bone and cartilage (Bancroft et al. 2002;Gemmiti and Guldberg 2006). The exact nature of this dependence is not known and will depend on the cell type in question. To capture the qualitative aspects of this dependence, we set so that the cell proliferation rate increases with shear stress up to a maximal level, but (as would be expected) high levels of shear stress result in a reduced proliferation rate (see Fig. 2). The coefficients in (4.6) have been chosen to capture the possible behaviours of the cell population within the range of computed shear stresses in our model.
Results from this case study can be seen in Fig. 3, where the steady-state cell population distribution obtained with F(S) as in Eq. (4.6) is compared to that when F(S) = 1 (denoted the shear-insensitive population, i.e. when the fluid shear has no effect on the cells) for a range of upper fluid layer flow rates Q f,in . From Fig. 3a, b, we see that for the lowest flow rate of Q f,in = 0.1, the cell volume fraction of the shear-sensitive (hereby called shear-enhanced) population is less than that of the shear-insensitive population for 0 < x < 1, but then increases above the shear-insensitive volume fraction once Q f,in is increased to 0.5. This corresponds to the shear stress exceeding the value at which F(S) is increased above 1 (around 2.5) and is to be expected given our chosen form of F(S). The increase in θ n is observed across the whole bioreactor domain, instead of at a specific point, as a result of the relatively small variation in θ n (and hence θ w ) for 0 < x < 1, which in turn means that there is little spatial variation in S. At the intermediate flow rate (Q f,in = 0.5), we also see that the shear-enhanced population becomes much more spatially uniform than the shear-insensitive cells. For the highest flow rate Q f,in = 1, the shear-enhanced θ n decreases at each point in x as the shear stress reaches levels which decrease the cell proliferation rate. In addition, greater up-and downstream variation is seen in θ n for both populations, but the changes are less extreme for the shear-enhanced cells. Hence, overall, the effect of the shear stress is to 'smooth out' the cell population: an increase in θ n corresponds to a decrease in θ w , and thus a higher value of S. Thus, regions where θ n is relatively high experience a fluid shear stress sufficient to decrease F(S) and hence the proliferation rate, whilst regions  Table 2 where θ n is relatively low experience a lower shear stress and therefore increased proliferation rate.
The effect of shear on the cell yield and population uniformity can be seen in Fig. 4, and the plots support the findings discussed above. The mean μ and standard deviation σ of θ n are plotted for a range of values of Q f,in , for the shearenhanced and shear-insensitive cases. This clearly shows that the shear-enhanced population has a lower yield and greater standard deviation than the shear-insensitive population for low Q f,in values, but that this relationship is reversed when Q f,in increases beyond a critical value just below 0.2. This switch corresponds to the critical value of S being reached at which F(S) > 1, and the cell proliferation rate is enhanced. As can be seen in the graphs for both μ and σ , the change occurs relatively abruptly; this is due to the fact that (as mentioned above) there is relatively little variation in θ n at this value of Q f,in , and so the fluid shear stress is relatively constant along the domain. The standard deviation of the shear-enhanced population remains lower than the shearinsensitive population for all subsequent values of Q f,in , but the yield drops below that of the shear-insensitive population once another critical value of Q f,in is reached, around 0.7, when the value of F(S) falls below 1. Thus, between  Table 2 these critical values (given by vertical dashed lines in Fig. 4), there is an optimum range of Q f,in within which shear stress helps to promote cell yield and also population uniformity. Within this range, the most uniform population is obtained at the (optimal) flow rate of Q f,opt = 0.2048. We also note that, for a flow rate above the first critical value, the shearenhanced population is always more spatially uniform than the shear-insensitive population for a given yield.

Shear stress-limited proliferation
In the second case study, we take F(S) to be of the following form so that cell proliferation is inhibited for sufficiently high values of fluid shear stress (see Fig. 2). Cells that have been  Table 2 shown to be sensitive to shear stress in this way include smooth muscle cells (Papadaki et al. 1996) and certain types of endothelial cells (e.g. human aortic endothelial cells Imberti et al. 2002 and human umbilical vein endothelial cells Akimoto et al. 2000). As in the previous section, the exact form of the cells' shear stress dependence is unknown, and our choice of F(S) and corresponding parameter values have been chosen so that our model captures the range of possible behaviours. Throughout this section, the cell population with proliferation rate limited by shear stress will be known as shear-limited. Figure 5 shows the results obtained from this case study. For low values of Q f,in , the shear stress has no effect on cell proliferation (as F(S) is equal to 1) and thus the shear-limited and shear-insensitive population distributions are the same. As Q f,in is increased, however, the downstream region, where θ n is greater, experiences a higher level of shear stress, and thus a decrease in proliferation rate of the shear-limited population. This results in the shear-limited population becoming more uniform than the shear-insensitive population. In contrast to the previous section, we note that regions of lower cell volume fraction are unchanged by the shear sensitivity, since in this case, proliferation is never enhanced by the fluid shear. Thus, the increased uniformity of the shearlimited population is at a cost of a lower overall volume fraction. The plots of μ and σ in Fig. 4 show that there is just one critical value of Q f,in in this case (around 0.4, given by the vertical dash-dotted line in Fig. 4), below which the shear-limited and shear-insensitive cases are indistinguishable. Above this, the shear-limited population is more uniform than the shear-insensitive population, but also has a lower yield. However, Table 3 demonstrates that the percentage decrease in σ is more significant than the percentage decrease in μ for 0.4 Q f,in 0.8, and thus, this provides a potentially optimal flow rate range when a uniform cell population is more important than a high cell yield.

Sensitivity to cell layer height
Finally, we investigate the sensitivity of our results to the cell layer height h 3 , which may be prescribed experimentally. In each case, we compare the cell volume fraction θ n and note that larger values of h 3 correspond to greater overall cell numbers, and vice versa for smaller values of h 3 . The cell layer height also affects the pressure gradient in the cell and upper fluid layers, since ∂ p w /∂ x depends linearly on Fig. 6). Figure 7a, b shows plots of the mean μ and standard deviation σ of θ n for a range of values of Q f,in for a thinner cell layer than Fig. 4, whilst Fig. 7c, d shows the corresponding plots for a thicker cell layer. In Fig. 7a, we see that the optimal ranges for Q f,in are now higher. The corresponding optimal flow rate is also higher than before, with Q f,opt = 1.048. This is a result of the form for |∂ p w /∂ x|, which increases as h 3 does as stated above. The shear stress therefore also decreases as h 3 decreases, and we thus need higher flow rates in order to obtain the same behaviour from the cell population. As expected, this trend is reversed in Fig. 7c when h 3 takes a higher value, with an optimal flow rate in this case of Q f,opt = 0.085. In addition, for Q f,in greater than around 1, we see the cell population decreases significantly, no matter what the shear stress dependence. We therefore conclude that a thinner cell layer makes the cell population less sensitive to the upper fluid layer flow rate, as higher flow rates can be used before compromising cell yield, and vice versa for a thicker cell layer.

Discussion
We have developed a multiphase model of fluid flow, solute transport and cell distribution in a simplified HFMB. This extends the work presented by Pearson et al. (2013), by considering the alternative set-up in which the cell-scaffold construct only partially occupies the ECS and there is an additional upper layer of free-flowing fluid. Fluid is pumped into the upper fluid layer via the upstream ECS port, and hence, this region experiences higher flow rates than when fluid is only supplied at the lumen inlet. It is therefore hypothesised that this could result in a greater fluid shear stress in the cell layer which could have an important, and possibly beneficial, effect on cell dynamics. This is confirmed from our model results which imply that the fluid shear stress in the cell layer is dependent on both the upper fluid inlet flux and cell layer width, both of which are varied as part of our investigation.
We presented results from two different case studies. In the first (relevant to osteoblasts, chondrocytes or endothelial progenitor cells), cell proliferation rate was enhanced by fluid shear stress, but only up to a certain critical level. Results demonstrated that the mechanosensitive cell population required a certain minimum ECS flow rate in order to reach yields comparable to the shear-insensitive population. For flow rates that resulted in fluid shear within the 'enhancing' range, however, the mechanosensitive population had a higher yield and was more uniform than the shearinsensitive cells. This uniformity persisted even once the shear levels became detrimental to the proliferation rate and the yield decreased. In the second case study, fluid shear had no effect on cell proliferation until a critical value was surpassed, after which it decreased the proliferation rate. This could be applicable to smooth muscle cells, or certain types of endothelial cells. Results again showed that the mechanosensitive cell population was more uniform than the shearinsensitive cells, but always had a lower cell yield. However, an ideal flow rate range could again be found, within which the mechanosensitive population is more uniform than the shear-insensitive population, and the yield has not yet been significantly compromised. Finally, we considered the effect of the cell layer width on these results and found that a thinner cell layer experiences a lower level of shear stress than a thicker layer, and hence can withstand higher flow rates.
The results presented here are qualitative in nature, demonstrating the range of expected behaviours in each sit- Fig. 7 Plots of the mean μ and standard deviation σ of θ n , for a range of Q f,in values in the shear-insensitive, shear-enhanced proliferation and shear-limited proliferation cases, for (a, b) h 3 = 0.5 and (c, d) h 3 = 1.5. The vertical dashed lines again indicate the optimal flux range for the shear-enhanced case, and the vertical dash-dotted line the critical flux value for the shear-limited case. The lumen flux Q l,in = 0.1 and other parameter values are as in Table 2 (a) shear insensitive shear−enhanced shear−limited uation. In particular, the functional forms for F(S) are constitutive and relatively simple, having been chosen to clearly demonstrate the possible effects of mechanotransduction on cell distribution in a HFMB. Any choice of F(S) would need experimental validation to be applied to a specific experimental setup. There are currently no relevant experimental results to which we can compare the results from our mathematical model, as the modelling is ahead of the experiments. However, the model framework developed here is extremely flexible and permits interrogation of different forms for any of the parameter values and constitutive laws presented. Hence, once relevant data become available, the model could be readily adapted and applied to a particular cell population, and proper validation would be possible. This could then be used in a predictive way to stimulate future experimental work.
have assumed that the timescales for advection in the porous c l = c m , ∂c l ∂ y = φ ∂c m ∂ y on y = 1; on the membrane/cell layer interface (again substituting in for p m in place of v m ) (6.11) on the cell/upper fluid layer interface v n = 0, ∂u n ∂ y + ε 2 ∂v n ∂ x = 0, u f = 0, (6.12) c w = c f , θ w ∂c w ∂ y = ∂c f ∂ y on y = 1 + h 2 + h 3 ; and finally on the bioreactor top u f = v f = 0, ∂c f ∂ y = 0 on y = H. (6.13)