Multiphase modelling of the influence of fluid flow and chemical concentration on tissue growth in a hollow fibre membrane bioreactor

A 2D model is developed for fluid flow, mass transport and cell distribution in a hollow fibre membrane bioreactor. The geometry of the modelling region is simplified by excluding the exit ports at either end and focusing on the upper half of the central section of the bioreactor. Cells are seeded on a porous scaffold throughout the extracapillary space (ECS), and fluid pumped through the bioreactor via the lumen inlet and/or exit ports. In the fibre lumen and porous fibre wall, flow is described using Stokes and Darcy governing equations, respectively, while in the ECS porous mixture theory is used to model the cells, culture medium and scaffold. Reaction–advection–diffusion equations govern the concentration of a solute of interest in each region. The governing equations are reduced by exploiting the small aspect ratio of the bioreactor. This yields a coupled system for the cell volume fraction, solute concentration and ECS water pressure which is solved numerically for a variety of experimentally relevant case studies. The model is used to identify different regimes of cell behaviour, and results indicate how the flow rate can be controlled experimentally to generate a uniform cell distribution under regimes relevant to nutrientand/or chemotactic-driven behaviours.


Introduction
In our society, there is an increasing need for replacement tissue and organs as a result of damage due to trauma, disease or old age.In a 2001 review (Lysaght et al., 2001), it was reported that existing organ and tissue replacement therapies account for 8% of all medical expenditure worldwide, equivalent to an annual cost of $350 billion.Tissue engineering is a promising alternative, as long as it has the potential to reduce these costs and improve patient prognoses.In vitro, it has been described as the cultivation and growth of transplantable human tissue to restore, maintain or improve tissue function (Langer et al., 1993).One approach is to seed cells on a biodegradable or resorbable 3D scaffold and culture the resulting construct in a bioreactor.When the desired quantity and type of fully developed, functional tissue has been produced, it is transplanted into the patient and, over time, the scaffold is degraded or resorbed so that only tissue remains (Tabata, 2000;Stock et al., 2001;Cortesini, 2005;Dawson et al., 2008).Alternative in vivo methods are also used, where the initial cell-scaffold construct is transplanted directly to the appropriate site and left to develop, exploiting the body as a bioreactor (MacArthur et al., 2004;Stevens, 2008;Geris et al., 2010).
The advantage of tissue engineering over donor transplants is that replacement tissue can be grown from cells that are much more readily available than donor organs, therefore reducing waiting times for the transplant.Since the patient's own (autologous) cells are often used, there is a significantly reduced risk of rejection of the new tissue (Stock et al., 2001).This approach can also be used as an alternative to artificial implants such as hip or knee replacements, which have a limited lifespan and can cause allergic reactions (Pörtner et al., 2005).
Tissue engineering is not without its drawbacks, however.The immense cost of even small-scale experiments has limited progress in the field to date, and prevented widespread scale-up for clinical use.This is despite significant research being carried out: e.g.there are over 5600 PubMed articles from the last 15 years on cartilage tissue engineering alone, yet no routine clinical solutions.The wide range of cell types under investigation, each with its own specific properties, also means that there is no one method that will work for all tissues.Issues to consider include how to seed the cells in order to control tissue generation and to ensure effective mass transport throughout the system.Furthermore, the correct choice of scaffold material and bioreactor geometry is vital in enabling the affordable production of a large enough quantity of viable cells.Cells are extremely sensitive to their surroundings, and internal bioreactor conditions must allow cells to be appropriately stimulated to mimic the in vivo environment, while at the same time avoiding cell damage or death (Stock et al., 2001;Martin et al., 2004Martin et al., , 2005;;Pörtner et al., 2005).As a result of this, a vast array of bioreactor designs and materials exist, each needing a different set of operating conditions to successfully grow a particular tissue (Pörtner et al., 2005).

Hollow fibre membrane bioreactors
We focus on hollow fibre membrane bioreactors (HFMBs), which consist of a cylindrical, hollow glass module with an exit port at the up-and down-stream ends (see Fig. 1).A synthetic, porous hollow fibre (which mimics a capillary in vivo and is impermeable to cells) is inserted through the centre and fixed to the module at the ends.The choice of fibre materials depends on the desired characteristics, e.g.pore size, porosity and hydrophilicity.Scaffold materials vary, but required features include controlled degradation, retention of shape, effective nutrient delivery and biocompatibility (Murphy et al., 2000;Khademhosseini et al., 2006;Dawson et al., 2008).Cell seeding involves suspending the cells in a collagen gel solution, which is then injected through the exit ports, and the gel left to set (see Ye et al., 2007 for more details).This method results in a cell-seeded natural scaffold throughout the extracapillary space (ECS).Alternatively, the cells can be seeded on the outer surface of the fibre in a similar manner (Ellis, 2006).
Culture medium (containing multiple nutrients and growth factors) is then driven through the fibre lumen and/or the exit ports by a pump at a prescribed flow rate.In the case of lumen-driven flow, the medium reaches the cells via the porous fibre walls and flows out through the end of the fibre lumen (the retentate) or either of the exit ports (the permeate).The ratio of retentate to permeate is controlled by imposing a downstream pressure via a clamp at the lumen exit (the exit ports being left at atmospheric Fig. 1.Top: photograph of a single HFMB module (ruler scale in cm).The exit ports (one of which is closed off here) can be seen, along with the hollow fibre which runs through the centre of the module.Cells are seeded in a natural scaffold throughout the ECS.Middle: cross-section of the boxed region (not to scale).The lower half (below the x-axis) is excluded, and the upper half is the idealized 2D modelling region with the x-axis running along the lumen centreline.Bottom: cross-section detail of a typical hollow fibre membrane.Images courtesy of Dr Marianne Ellis, Centre for Regenerative Medicine, University of Bath.pressure).The porous fibre wall, or membrane, separates the cells from the main flow in the lumen so that relatively large flow rates can be used without exposing cells to high fluid shear stresses.In addition, the alternative flow route in via the exit ports allows greater control of the cell environment.The great advantage of HFMBs is that the surface area to which culture medium is supplied is large compared with the bioreactor volume, providing efficient delivery conditions.As a comparison with standard flask culture, a 1 m 3 flask would be required in order to grow the same number of cells as in a 0.5 l HFMB.This bioreactor, therefore, has great potential as a method to generate replacement tissues on a clinical scale (Shipley et al., 2010).However, this is still very much an open research area, and outstanding problems include how best to initially seed the cells to obtain the desired final distribution and to find optimal flow rates to ensure sufficient delivery of nutrients to, and removal of waste products from, the cells.Only when these issues have been resolved will such procedures present feasible future alternatives to current tissue replacement strategies.

Multiphase modelling
Given the considerable cost and complexity of tissue engineering experiments, it is clear that mathematical modelling of the systems could provide great benefit.In particular, since many bioreactor experiments are still at a relatively early stage of development, a rare opportunity exists to use mathematical modelling to inform and direct experiments.In the interest of brevity, we summarize the existing literature most relevant to the model presented here; for a more complete overview, we refer the reader to the recent review by O' Dea et al. (2013).
We adopt a mixture theory approach, which can be applied to a system of two or more interacting and interpenetrating phases; see Whitaker (1967Whitaker ( , 1973)), Gray (1975), Drew (1983) and Quintard et al. (1994) for more details.Microscale governing equations for each phase are averaged over a suitable representative volume, resulting in a continuum model that describes the macroscale dynamics of each phase.This technique allows the behaviour of inherently complex systems to be modelled relatively simply, with the possibility of solution of the governing equations by analytical and/or numerical methods.Individual phases can be tracked independently on the macroscale without needing to resolve microscale dynamics, but still allowing the incorporation of important microscale effects.These effects, along with interactions between phases, are included via constitutive laws.As well as tissue engineering, the approach has been used to model tumour growth (Byrne et al., 2003;Franks et al., 2003), biofilms (Cogan et al., 2004), cell motility (Alt et al., 1999;Byrne et al., 2004;King et al., 2005;Kimpton et al., 2013) and even soil dynamics (Theodorakopoulos, 2003).
In vitro tissue growth has been modelled using mixture theory by Lemon et al. (2006).The scaffold, cells and culture medium (modelled as water) are all treated as separate phases.Constitutive equations are chosen to describe the interactions between the rate of production of and the stress terms in each phase, and the linear stability of the equilibrium state of the system is analysed.This analytical approach is used to show that cells seeded into a scaffold in a static set-up can either disperse or aggregate depending on the relative strengths of cell-cell interactions and cell-scaffold tractions, and the results agree qualitatively with those from experiments.This work is built upon in O' Dea et al. (2010) and Osborne et al. (2010), where the mixture theory approach is applied to a perfusion bioreactor setup.This consists of a tissue-scaffold construct placed in a cylindrical chamber through which culture medium is driven using a pump.In O' Dea et al. (2010), the dependence of cell population behaviour on fluid shear stress, pressure and local cell density is analysed in a 2D Cartesian geometry.The governing equations are simplified in the limit where the aspect ratio of the bioreactor channel is small, and then solved using both analytical and numerical approaches.The results show that cell movement depends greatly upon cell density and the relative importance of cell aggregation and repulsion, while different mechanotransduction mechanisms produce qualitative variations in the composition of the cell-scaffold construct.Osborne et al. (2010) add to this by solving the full system from O'Dea et al. (2010) numerically using finite element methods.They find that the long-wavelength approach used in O' Dea et al. (2010) is good at predicting average behaviour; however, to determine detailed spatial distribution of cells, it is necessary to consider solutions of the full system.
The model presented here is 2D for simplicity to enable analytical progress.It is expected, however, that results from a 3D (axisymmetric) set-up would be qualitatively similar; we hope to verify this in future work.The flow in the lumen is modelled by the Navier-Stokes equations.In the membrane, the flow is governed by Darcy's law, which we put into the context of the multiphase formulation in order to give self-consistent boundary conditions.In the ECS the governing equations are derived using mixture theory, and we model the scaffold (comprised of a natural matrix such as collagen), culture medium and cells as three distinct phases.Interactions between these phases are incorporated via extra pressure, force and drag terms in the governing equations.In particular, we are concerned with the dynamics of the cell population and how they depend on experimentally controllable aspects of the bioreactor environment such as flow rate and nutrient concentration.

Paper outline
In Section 2, we describe the model setup, outlining the governing equations and boundary conditions in each region in Sections 3.1-3.4.The resulting system is non-dimensionalized in Section 4.1, using a lubrication approach to exploit the small aspect ratio of the bioreactor.In Sections 4.2-4.4,we use systematic asymptotic methods to reduce the system from 16 to 3 unknowns.This results in a coupled system of partial differential equations (PDEs) at leading order in the lumen aspect ratio for the cell volume fraction, global concentration of solute and reduced water pressure, along with analytical expressions for the remaining variables.Typical dimensionless parameter values for this model are given in Section 4.5.In Section 5, we present and compare a number of different case studies relating to possible experimental setups, including nutrient-limited cell growth, response to a delivered or produced chemoattractant and backflow through the ECS.We consider the effect of flow rate on cell yield and distribution and determine optimal flow rates for each case study in order to achieve a spatially uniform cell distribution.We conclude by discussing our findings in Section 6.

Model description
We consider an idealized problem, motivated by taking a cross-section of an HFMB module and excluding the end regions containing the exit ports.The central domain is modelled in 2D Cartesian coordinates (x, y), with the x-axis along the lumen centreline (which is assumed straight) and the y-axis pointing upwards (see Fig. 1).This makes the problem amenable to analytic techniques while still capturing key features of the bioreactor.For simplicity, we consider the upper half of the bioreactor only, so that the modelling region 0 < x < L, 0 < y < H consists of three rectangular compartments: the lumen, porous membrane and ECS.The exact up-and down-stream boundary conditions will be chosen to mimic the up-and down-stream exit port configurations for a specific experimental setup.We first consider the setup shown in Fig. 1, in which only the downstream exit port is open.In Section 5.6, we briefly look at the case where both exit ports are open and there is a flow driven through the downstream exit port in the opposing direction to the flow in the lumen.
In the first flow regime (where only the downstream exit port is open), the culture medium (modelled as water) is pumped into the lumen and, due to an imposed pressure at the downstream end, flows out either through the lumen or through the ECS.The up-and down-stream ends of the porous membrane are sealed so that no fluid or cells can leave the system there.The ECS is comprised of the cells (one cell type only), culture medium and a natural scaffold.We assume that there is no deposition or degradation of this scaffold (so that it is fixed and inert) and that initially the cells are seeded uniformly throughout the ECS.In addition, the porous membrane is impermeable to the cells so that they cannot leave the ECS.When considering interphase pressures in the ECS, we include tractions between the cells and the scaffold only, assuming that tractions between the cells and fluid, and between the fluid and scaffold, are negligible by comparison.We assume also that the drag between the water and cells is much smaller than that between the scaffold and water or between the scaffold and cells.These simplifying assumptions are motivated by the fact that we expect the cells to strongly adhere to the scaffold.The culture medium contains a solute of interest, which is either a nutrient (e.g.oxygen) or a chemoattractant (a chemical that induces a motile response from the cells; e.g.see Postlethwaite et al., 1987;Giannobile, 1996;Chen et al., 2003).The simpler geometry and small aspect ratio of the domain motivate the non-dimensionalization in Section 4.1 and allow us to derive a system of governing equations in the lumen, porous membrane and ECS which can be studied analytically.

Governing equations
Now we outline the governing equations for each of the three regions of the bioreactor: the lumen, the porous membrane and the ECS.We denote by h 1 , h 2 and h 3 the heights of the lumen, porous membrane and ECS, respectively, with the total height of the half-bioreactor given by H = (h 1 + h 2 + h 3 ).The total bioreactor length is denoted by L, and we note that h 1 , h 2 and h 3 are comparable, with h 1 /L 1, so that the aspect ratio of each region is small.We state relevant boundary conditions in Section 3.4.Parameter values are determined later in Section 4.5.

Lumen
The lumen is defined as the region 0 < x < L, 0 < y < h 1 .Here, there are no cells, so we consider the flow of water and transport of solute only.The former is modelled using the incompressible Stokes equations.Conservation of mass gives ∇ • u l = 0, (3.1) where u l = (u l , v l ) is the water velocity in the lumen, and the density of water, ρ w , is assumed constant.We exploit the fact that inertial effects in the lumen are negligible (this is confirmed in Section 4.5), so that conservation of momentum yields where σ l is the stress tensor for the water in the lumen, given by pl is the water pressure in the lumen, μ w is the water viscosity (assumed constant) and g = (0, g) is the acceleration due to gravity.In this work, we shall consider flow rates that are towards the lower end of those employed experimentally (see Table 1), and we must account for the effects of gravity in this instance.We introduce the reduced water pressure p l , defined by p l = pl + ρ w gy, and subsequently similarly define reduced pressures in the porous membrane and ECS.Conservation of mass for the solute concentration is given by the advection-diffusion equation where t is time, c l is the solute concentration per unit volume in the lumen and D is the diffusion coefficient for the solute in water (assumed constant).

Porous membrane
The membrane occupies the region 0 < x < L, h 1 < y < h 1 + h 2 .Here, we model the flow using Darcy's law for flow in porous media, which takes into account the microscale effects of viscosity on the macroscale: where u m = (u m , v m ) is the water velocity in the porous membrane, k is the porous membrane permeability (assumed constant) and p m = pm + ρ w gy is the reduced water pressure in the membrane.Furthermore, conservation of fluid mass yields where we have assumed that the porous membrane pore fraction φ is constant.Similar to the above, conservation of mass for the solute concentration is given by the advection-diffusion equation where c m is the solute concentration per unit volume of water in the membrane and we have assumed that the effective diffusion coefficient for the solute in the membrane is proportional to the pore fraction.

Extracapillary space
Finally, we consider the ECS, occupying 0 < x < L, h 1 + h 2 < y < H and comprised of cells, water and scaffold.In prescribing the governing equations for the multiphase flow in the ECS, we follow closely the works of Lemon et al. (2006) andO'Dea et al. (2010).We denote the cell, water and scaffold volume fractions by θ n , θ w and θ s , respectively, where θ s is assumed to be constant.We apply the no voids constraint that Conservation of mass for the cells and water is given by where u n = (u n , v n ) and u w = (u w , v w ) are the cell and water velocities, respectively, J n represents the net cell production rate and J w represents the net water production rate.Here, we have made the assumption that the cells and water are of equal and constant density, i.e. that ρ n = ρ w , which is a good approximation since cells are mainly comprised of water (Lemon et al., 2006).We impose J w = −J n so that mass is conserved.To keep the analysis as general as possible, we assume only that J n is a function of θ n and c w (the solute concentration per unit volume of water in the ECS), prescribing a functional form in Section 5 for each experimentally motivated case study that is considered.Conservation of momentum for the two phases (neglecting inertia) is given by where σ n and σ w are the macroscale stress tensors for the cell and water phases, respectively, and f ij (i, j = n, w, s, i | = j) represents the interphase force exerted on the ith phase by the jth phase.By Newton's third law, we have f ij = −f ji .
We pose constitutive forms for σ i and f ij , following the approach of Lemon et al. (2006).Firstly, we model the cell phase as a viscous liquid, and so let where pn is the cell intraphase pressure, I is the identity tensor, μ n is the effective viscosity of the cell population and τ n is the deviatoric stress tensor for the cell phase, given by where superscript T denotes transpose.We decompose the cell pressure pn into the water pressure pw plus an extra component Π , (3.13) and for now we assume only that Π is a function of θ n and c w .In Section 5, motivated by O'Dea et al. (2010) and Byrne et al. (2004), we split Π into two components Π n and Π c which, respectively, account for cell-cell interactions (such as osmotic stress and surface tension in the cell membranes) and chemotaxis, and take the following forms: where ν, δ a , χ and c are positive constants.The first term in Π n represents the fact that cells tend to aggregate at low densities, i.e. when θ n is small, while the second models the repulsive forces which arise due to contact inhibition (O'Dea et al., 2010).The function Π c represents chemotactic movement up gradients in the solute, with a greater effect in areas of lower concentration, a behaviour which has been observed experimentally by Jeon et al. (2002).
For the water phase, we assume the stress tensor so that the force balance (3.10b) and constitutive laws prescribed below for the interphase forces (see (3.16)) lead to the transport of water being governed by Darcy's law.
Next we impose constitutive laws for the interphase forces f ij in (3.10) of the form where p ij = p ji (i, j = n, w, s) represents the interphase pressure exerted on the ith phase by the jth phase.Each interphase force f ij is decomposed into three components: see Lemon et al. (2006) for details.
The first represents the interfacial force exerted on the ith phase by the jth phase.The second term of each force represents the reaction to the corresponding interfacial force exerted on the jth phase by the ith phase.Finally, the third term represents drag between the phases, where γ ij is the drag coefficient (assumed constant).We note that any forces involving the scaffold phase will be simplified by the fact that ∇θ s = u s = 0 due to our assumption that the scaffold is rigid with constant volume fraction.
The pressure p ij = p ji (i, j = n, w, s) is assumed to consist of the water pressure pw , which is 'contact independent', and an extra pressure ψ ij due to the tractions between the ith and jth phases: (3.17) Since we consider the tractions between the cells and scaffold only, we can simplify the system by taking ψ ij = ψ ji = 0 unless i = n and j = s.We take ψ ns to be a constant function: where η is a positive constant representing the cells' affinity for the scaffold.We note that this is a simplified version of the form used in O' Dea et al. (2010) which depends on θ n , and we will keep the more general form ψ ns = ψ ns (θ n ) until Section 5 in order to keep our derivation as general as possible.
We define reduced pressures in the ECS as in the lumen and porous membrane by setting Making use of the constitutive relations (3.11-3.13)and (3.15-3.17),neglecting drag between the cells and water as stated in Section 2 and substituting in for p w from (3.19), the resulting governing equations for the cell phase in the ECS are given by and for the water phase are We note that in (3.21b) we have indeed recovered Darcy's law, which models the microscale effects of the water viscosity on the macroscale.Now we return briefly to the membrane flow model to define the stress tensor in the context of the multiphase formulation.We note that within this framework Darcy's law corresponds to the force balance (3.23) this being the only choice for σ m which is consistent with the choice for σ w above.Finally, conservation of mass for the solute in the ECS is given by where we have assumed that mass transport only occurs in the water phase.We have assumed that the effective diffusion coefficient for the solute in the water phase is proportional to the water volume fraction.Here, R is a reaction term which may depend on θ n and/or c w .Specific forms will be given in Section 5, where it accounts for either the uptake of nutrient, or the production or uptake of chemoattractant.
Across the three regions, we have a system of 16 equations for the 16 unknowns u i , v i (i = l, m, n, w), p i , c i (i = l, m, w), θ n and θ w .Now we impose boundary conditions to close the problem.

Boundary conditions
In the following, n l = n m = n e = (0, 1) are the upward pointing normals to the lumen, porous membrane and ECS, respectively.
As we consider the upper half of the bioreactor only, we impose symmetry and no flux conditions on the lumen centreline so that On the lumen-membrane interface, we impose no slip, as well as continuity of water flux, of normal stress, of concentration and of solute flux so that

403
where we have assumed that all of the fluid stress from the lumen is taken up by the fluid in the membrane so that condition (3.26c) reduces to continuity of (reduced) pressure in the leading-order model (the gravitational contributions will cancel).The assumption of continuity of pressure at the boundary between flow through materials of different permeabilities is also taken in Beavers et al. (1967) and Joseph et al. (1964).The no slip condition (3.26a) is motivated by the work of Shipley et al. (2010) who found that slip is insignificant for these membranes.
The membrane-ECS interface is impermeable to cells as the mean pore radius (≈ 1 µm) is much smaller than that of a cell (≈ 5 µm), and hence we impose no flux of cells here.In addition, we impose no slip of cells, along with continuity of flux of water, of normal stress, of concentration and of flux of solute.We make the constitutive assumption that all the stress is taken up by the water in the ECS, not the cells, which again ensures continuity of (reduced) water pressure across the membrane-ECS interface at leading order in the reduced model.These give The bioreactor boundary y = H is a solid boundary on which we impose no flux and slip of cells, and no flux of water or solute, so that Boundary conditions at the up-and down-stream ends of the bioreactor are chosen to represent conditions in the excluded regions near the exit ports and at the lumen inlet and outlet, and will be prescribed as required in Section 4.3.Initial conditions are imposed as necessary in Section 5.

Model reduction
The choice of non-dimensionalization is motivated by lubrication theory.This scaling is appropriate in the lumen, due to its small aspect ratio and negligible inertial effects, and results in a reduced pressure which is independent of y, allowing us to explicitly solve for the flow velocity there.The ECS also has a small aspect ratio and negligible inertia, and hence this scaling is physically relevant in the multiphase region too.This leads to the reduced pressure term dominating in the y-force balance in the water phase, and both the reduced pressure and intra-/inter-phase pressure terms dominating in the y-force balance in the cell phase.As a result, both the reduced pressure and cell volume fraction are independent of y at leading order in the ECS.This allows significant simplification of the governing equations and results in a coupled system of three second-order non-linear PDEs at leading order in the lumen aspect ratio to solve for the cell volume fraction, solute concentration and reduced water pressure in the ECS, along with corresponding expressions for the remaining leading-order unknowns.

Non-dimensionalization
We non-dimensionalize as follows, using hats to denote dimensionless variables: where ε = h 1 /L is the (small) aspect ratio of the bioreactor and C * , R * and U * are (respectively) typical concentration, reaction and porous membrane/ECS velocity scales.We have multiplied the velocity scales in the lumen by a factor of μ n /μ w so that the same pressure scales can be applied throughout the system.As previously stated, we expect inertial effects to be negligible compared with viscosity, and hence have employed a viscous pressure scaling throughout.The resulting horizontal velocity scale in the lumen, μ n U * /μ w , represents the typical horizontal velocity of the fluid flowing through the lumen as determined by the imposed lumen inlet flow rate.We choose the timescale L/U * and take the inverse of this as the scale for J n so that the timescales for advection in the ECS/porous membrane regions and proliferation are assumed comparable.We also choose the reaction rate timescale so that R * = DC * /L 2 without loss of generality.For now we consider the biologically plausible regime where the viscosity of the cell phase (on the macroscale) is much greater than that of water due to the presence of the cytoskeletal network, the cells' adherence to the scaffold and the effect of cell-cell interactions (on the microscale).In particular, we consider the limit in which μ w /μ n = λε for some λ = O(1) as ε → 0. We note that other limits may also be physically relevant, but would involve solving a much more complex system.Relevant dimensionless parameter values will be determined in Section 4.5.
In the following, we move into component form, substituting for the reduced pressures throughout and eliminating v m in favour of ∂p m /∂y in the boundary conditions (3.26b) and (3.27b) using the y-component of (3.5).Dropping hats, the dimensionless system of equations in the lumen (for 0 < y < 1) is therefore where Pe = LU * /(λεD) is the Péclet number for axial flow in the lumen.Similarly, the porous membrane equations become (for 1 < y < 1 + h 2 ) and  where κ = k/(λε 5 L 2 ) is the dimensionless permeability.We will assume a dominant balance between the radial terms in (4.4b), and so take κ = O(1) (see Table 2).
In the ECS (1 + h 2 < y < 1 + h 2 + h 3 ), the no voids condition is unchanged: The dimensionless conservation of mass equations are conservation of momentum for the cell phase is given by and for the water phase is where Finally, the boundary conditions on the lumen centreline become ∂u l ∂y = 0, v l = 0, ∂c l ∂y = 0 on y = 0; (4.12) on the lumen-membrane interface on the membrane-ECS interface and on the bioreactor boundary Next we use asymptotic analysis to derive a reduced system of governing equations at leading order in ε, before introducing up-and down-stream boundary conditions to close the problem.We begin by expanding each dependent variable as an asymptotic series in powers of ε.For example, we let u l ∼ u l 0 + εu l 1 + ε 2 u l 2 + • • • and similarly expand the remaining velocities v l , u i , v i (i = m, w, n); the reduced pressures p i (i = l, m, w); the concentrations c i (i = l, m, w) and the volume fractions θ n and θ w .We will consider the regime where Pe = O(1) (see Table 2), corresponding to a dominant balance between axial diffusion and advection in the lumen.In the following, we omit subscript 0 from the leading-order variables except where needed for clarity.
Substituting the expansions into lumen equations (4.2) and ( 4.3) and boundary conditions (4.12) and (4.13a,b) and equating coefficients of ε 0 gives with ∂u l ∂y = 0, v l = 0, ∂c l ∂y = 0 on y = 0 (4.17) and which yields In the porous membrane, governing equations (4.4) and (4.5) become for some function P(x, t) which is determined as part of the solution to the reduced model.
In the ECS (1 + h 2 < y < H), equating coefficients of ε 0 in (4.6) yields and From (4.27c) and boundary conditions (4.28e,f) and (4.29d), we find that the solute concentration is uniform across the bioreactor at leading order, with c l = c m = c w = c(x, t), for some global concentration c(x, t) which is also determined as part of the solution to the reduced model.
Equation (4.27b) tells us that the reduced water pressure p w is independent of y, and thus it follows from (4.26) and (4.23), respectively, that the leading-order cell and water volume fractions θ n and θ w are also independent of y.In addition, from (4.25) and boundary conditions (4.28a) and (4.29a) we can deduce that (4.30) where (4.31) where the cell flux Q n and water flux Q w are given by Now we determine an equation for c before introducing up-and down-stream boundary conditions to close the resulting coupled system.Equation (4.11) implies that ∂c w /∂y is of O(ε 2 ) and therefore ∂c w 1 /∂y = 0. Hence, equating coefficients of ε 2 in the solute equations in each region (briefly returning to subscript 0 notation for leading-order variables for clarity), we obtain The corresponding boundary conditions are given by ∂c l 2 ∂y = 0 on y = 0, (4.39) and We integrate (4.36-4.38)with respect to y across the lumen, porous membrane and ECS, respectively, and add them together (multiplying the integrated form of (4.37) by the pore fraction φ first).Applying boundary conditions (4.39), (4.40a), (4.41a) and (4.42) then yields where the solute flux Q c is given by (4.44)

Boundary conditions for the reduced model
To close the coupled system given by (4.33a,b) and (4.43), we prescribe up-and down-stream boundary conditions appropriate to this flow regime.In order to determine these boundary conditions exactly, it would be necessary to determine the local flow dynamics just outside the edges of our modelling region by solving the 2D problem for the cross-section of the whole bioreactor.This would be extremely computationally expensive, and for our simplified setup, we instead impose conditions that mimic the expected behaviour in these up-and down-stream regions.We note that this approach allows a rapid exploration of parameter space.Firstly, we prescribe a volume flux into the lumen, setting for a non-dimensionalized volume flux per unit length in the z-direction (perpendicular to x and y) Q l,in , which is assumed constant.Additionally, we anticipate a small flux into the ECS (from fluid which has already permeated through the fibre from the lumen), and set where Q e,in (assumed constant) is the non-dimensionalized volume flux per unit length in the z-direction.
Here, we fix a value that is much less than the flux into the lumen, motivated by the expectation that the permeate flow rate in the ECS will be much lower than that in the lumen.We also carry out a sensitivity analysis of this and other key parameters in Section 5.5.Next, we prescribe a downstream lumen pressure (controlled experimentally by a clamp as discussed previously), setting p l = P d at x = 1, (4.47) where P d is dimensionless, fixed and constant.At the downstream ECS end, we mimic the atmospheric pressure conditions at the bioreactor exit ports, and so take (given the choice of non-dimensionalization in (4.1)) We impose no flux of cells out of the reduced region by setting This mimics the fact that filters are used at the exit ports (which lie just outside our modelled region) to avoid losing cells out of the bioreactor.
Depending on whether the solute is being delivered at the lumen inlet or not, we impose either an experimentally fixed concentration c in upstream (in which case C * = c in ) or no total flux of solute in.These give, respectively, c = 1 or Q c = 0 at x = 0. (4.50) We also require a concentration boundary condition at the lumen exit.Here, there is no constraint on the concentration in an experimental setup.Our choice of boundary condition is motivated by the fact that we would expect the concentration to be constant in space after it leaves the lumen due to the effect of diffusion.We therefore impose no diffusive flux out: We note that this is also the boundary condition which has the least influence on the concentration field (prescribing a downstream concentration, for instance, would have a large effect on the concentration elsewhere in the bioreactor) and that this still allows an advective flux out of the bioreactor.These boundary conditions allow us to solve explicitly for more of the leading-order unknowns in terms of θ n , c and p w , so that in the lumen, in the membrane, and in the ECS, where P is given in terms of p w by and α 1 , α 2 and M (x, t) are as given in (4.31).We note that (4.24a,b) still hold for v n and v w , respectively.

Summary of reduced model
The reduced model for θ n (x, t), c(x, t) and p w (x, t) is given by where the fluxes are defined by the boundary conditions are and the initial condition on θ n is taken to be for some initial profile θ * n (x).Thus, the original system ((3.1-3.8) and (3.20-3.24))has been reduced from 16 to 3 unknowns, and from 2D to 1D, having eliminated the dependence on y.We have an advection-reaction-diffusion PDE for θ n , and quasi-steady second-order non-linear ordinary differential equations (ODEs) for c and p w .The main effects contributing to the cell flux Q n are cell advection (dependent on the water pressure gradient and chemotaxis), cell diffusion (dependent on the cell intraphase pressure and cell-scaffold interphase pressure) and the cell-scaffold drag (through F).The water flux Q w contains an advection term that also depends on the water pressure gradient and waterscaffold drag.Finally, the solute flux Q c has an advective term dependent on the lumen flow rate, and a diffusive term that is affected by the widths and (effective) porosities of the membrane and ECS.This coupled system can now be solved numerically, given appropriate constitutive forms for J n and R.

Parameter values
Now we determine the sizes of the dimensionless parameters, guided by our experimental collaborators and dimensional values obtained from the literature where possible.This justifies the distinguished limit taken to obtain our leading-order system in the previous section.Typical dimensional parameter values are given in Table 1, with dimensionless parameters summarized in Table 2.The column entitled 'Restriction' in Table 2 states any asymptotic bounds on the parameter in question required for the validity of the analysis (compared with our small parameter ε = 0.002).We also discuss parameter choices for which either no data are available in the literature and therefore assumptions must be made, or there exist a variety of possible values.
Firstly, we set the membrane/ECS velocity scaling U * by using an experimental lumen flow rate of 3.86 × 10 −4 ml min −1 , which is equivalent in volume to that used in flask culture (in which 50 ml of culture medium is supplied every 3 days to cells covering a 75 cm 2 cross-section).Achieving such a small flow rate in an HFMB setup is one of the aims of current experiments being undertaken by our collaborators in Dr Marianne Ellis' group in the Centre for Regenerative Medicine, University of Bath.The corresponding lumen velocity is 5.12 × 10 −6 m s −1 and hence our membrane/ECS velocity scaling U * = 1.02 × 10 −8 m s −1 .It is not known whether such a low flow rate will guarantee sufficient supply of nutrient to the cell population, and this is something that we will verify as part of our results in Section 5. We set C * to be 0.22 mol m −3 , a typical inlet concentration for oxygen when culturing a variety of cell types (Shipley et al., 2012).When the type of solute is specified in Section 5, it is assumed that the typical concentration scalings for chemoattractant are similar to those for oxygen and so the same value for C * is used throughout.
Using the parameter values in Table 1, we note that the reduced Reynolds number in the lumen ε 2 Re = 2.05 × 10 −6 , which justifies the lubrication approximation taken in (3.2).We can also find the Reynolds number in the membrane, Re m = ρ w LU * /μ w = 0.001, which justifies the use of Darcy's law in this region.As mentioned in Section 4.2, we consider the regime in which Pe = O(1).This is appropriate, since, given the choice of flow velocity U * , the Péclet number Pe = 170.75.
In the membrane, we set κ = O(1) as discussed previously (in fact, κ ≈ 2.1 using the parameter values in Table 1).
In the ECS, although we do not have data indicating typical sizes of the drag terms, we take both ζ ns and ζ ws to be of O(1) so that their effects are retained at leading order in ε, but set ζ ws < ζ ns (as we would expect physically).In Π c , the coefficient χ is chosen so that the effects of chemotaxis appear at leading order, while the parameter c is taken to have a value comparable with a typical concentration.
A number of dimensionless parameters appear in the constitutive forms of J n and R in Section 5.In J n , the (dimensional) cell proliferation rate coefficient Γ nw is chosen 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 (values based on estimations by our experimental collaborators).In R, the dimensionless uptake coefficient ΓR1 (which is the same for nutrient and chemoattractant) is taken to be of O(1), as is the dimensionless chemoattractant-production-rate coefficient ΓR2 , so that we obtain a balance between uptake and production.Finally, the uptake constants K and K 1 which appear in J n and R, respectively, are taken to be comparable with a typical dimensionless concentration value.
Finally, to determine values for Q l,in , we have taken a plausible range of 3D flow rates (3.86 × 10 −6 − 3.86 × 10 −3 ml min −1 ).These values have been converted into 2D fluxes (in m 2 s −1 ) by dividing by the vertical length scaling εL.The dimensionless downstream pressure P d can take a range of values in an experimental setup and hence is set to be of order unity in our model.The exact value has been calculated using the relationship between flow rate, downstream pressure and permeate to feed ratio determined by Shipley et al. (2010) using a flow rate of 3.86 × 10 −4 ml min −1 and (dimensional) permeate to feed ratio of 0.005 (see Section 4.3).

Numerical results: case studies
Having developed a model for the leading-order cell volume fraction θ n , solute concentration c and reduced water pressure p w , we consider a number of experimentally relevant case studies which correspond to typical nutrient and chemoattractant scenarios.In Section 5.1, we consider a nutrient delivered at the lumen inlet which limits cell proliferation rate and neglect the effect of chemotaxis.In Sections 5.2 and 5.3, we consider a chemoattractant, respectively, delivered at the lumen inlet or produced by the cells, and in Section 5.4 we model both a delivered nutrient and a produced chemoattractant.Therefore, in each case, we choose different forms for the cell mass transfer term J n , reaction term R and inlet concentration boundary condition.The forms for these terms are chosen based on discussions with experimental collaborators, and results indicate the typical behaviour to be expected from the system in question, with the aim of informing bioreactor operating conditions (specifically, flow rates).We investigate the sensitivity of these results to various parameter values of interest in Section 5.5.Finally, in Section 5.6 we consider an alternative flow regime where there is a downstream flux into the lumen and a backflow in the ECS instead of a downstream flux into both.
Results are obtained using a finite-difference scheme implemented in MATLAB.The system (4.57-4.62)was discretized in x, resulting in a coupled system consisting of an ODE in time for θ n and algebraic equations for c and p w at each grid point (denoted by θ n,i , c i and p i ).This was then solved in MATLAB using the inbuilt ODE solver ode15s given an initial profile for θ n (see (4.63)).We found that for arbitrary, smooth initial conditions, the steady state obtained is independent of the choice of θ * n .For the results presented in this paper, we take θ * n to be the constant value 0.3 without loss of generality.A convergence test was carried out to ensure that the second-order scheme converges as expected as the step size Δx → 0. Once the accuracy of the scheme had been verified, we applied our general model to specific experimental setups.
In each experimentally motivated case study in Sections 5.1-5.4,we first investigate the steady-state distributions of the cells, the relevant solute(s) and the reduced water pressure for a number of different flow rates in a representative range of Q l,in for which our asymptotic solution is valid (see Table 2).The system is deemed to have reached a steady state at time T if the maximum change in θ n , c and p w between time T − 0.062 and T (corresponding to a dimensional time of around 1 week) is < 10 −4 , and hence the end times are different for each case study.We note that results from each case study within an experimentally relevant time frame are very close to the steady states shown here.The maximum change in each of θ n , c and p w between an experimental end-time solution and the steady-state solution is <10% for each case study.Specifically, the results from Section 5.6 are within 10% of the steady state after a dimensional time of around 2 weeks, those from Sections 5.2 and 5.3 after 3 weeks, those from Section 5.5 after 4 weeks and results from Sections 5.1 and 5.4 after 6 weeks.
Next, we consider the effect of flow rate on cell yield and cell distribution in each case study by finding the mean and standard deviation of θ n for a range of Q l,in values.In this work, the mean μ is given by and the standard deviation σ by where T is the relevant end time.The MATLAB function trapz is used to estimate μ and σ for a given θ n .This analysis is motivated by maximizing either the cell yield or spatial uniformity of the cell population in a particular experimental scenario.Focusing on the latter, for each case study we find the flow rate Q l,opt which results in the most uniformly distributed cell population by minimizing σ with respect to Q l,in .For this, as well as plotting Q l,in versus σ as described above, we use the inbuilt MATLAB function fminbnd to find the minimum standard deviation of θ n within the flow rate range considered, i.e. 0.01 < Q l,in < 10. Results for this optimal flow rate are presented alongside those for the full range of Q l,in values.By also considering the value of μ as the flow rate is adjusted, we are able to determine whether or not cell yield is being compromised by requiring a uniform cell distribution.
In Section 5.5, we vary a number of key parameters and investigate their effect on μ and σ .Finally, in Section 5.6 where there is a backflow in the ECS, we consider steady-state results for a range of Q e,in values, including the optimal flow rate Q e,opt for which the standard deviation σ of θ n is at a minimum.
In the case studies in Sections 5.1 and 5.4 which explicitly include nutrient, we check that the optimal flow rates result in a nutrient concentration that is above the minimum required for cell viability throughout the bioreactor.This minimum is taken to be a dimensional concentration of 8 × 10 −2 mol m −3 , which is at the upper end of the range of minimum oxygen concentrations required for functional cells (see Shipley et al., 2011) and here corresponds to a dimensionless concentration of 0.36.For the case studies in which nutrient concentration is not explicitly taken into account (Sections 5.2, 5.3), we assume that the cells are supplied with sufficient nutrient for all flow rates in the relevant range, and that there is enough nutrient to ensure that there is no limiting effect on cell proliferation rate.The first assumption is confirmed by the fact that the nutrient concentration in Section 5.4 (where both nutrient and chemoattractant are explicitly modelled) is above 0.36 for all flow rates considered.However, in Section 5.1 where only nutrient is considered, the nutrient concentration does fall below 0.36 at the downstream end of the bioreactor for the lowest flow rate Q l,in = 0.01.Hence, if such low flow rates are to be used, the cells may not have a sufficient supply of nutrient throughout the bioreactor.The second assumption is shown to not have a significant effect on trends in results.This is seen by comparing the cell distributions from Sections 5.3 (where nutrient is not included and assumed to be in abundance) and 5.4 (where nutrient concentration is explicitly modelled).Unless stated otherwise, all parameter values are taken from Table 2 and the (dimensionless) flux ratio Q e,in /Q l,in is kept constant and equal to 2.5.We note that although the kinetic parameter values used here are fairly generic, specific cell types could be explored as required.

Nutrient-driven proliferation
In the first case study, we consider a nutrient of concentration c(x, t) that is supplied at the lumen inlet and does not induce a chemotactic response from the cells.This could represent, for instance, hepatocytes in a liver sinusoid (which can be modelled using an HFMB for tissue engineering applications).Gradients in many chemical components, including oxygen, are set up along the length of a sinusoid.These gradients correlate with changes in metabolic cell function along the sinusoid, but do not induce cellular migration (Williams et al., 2013).The corresponding inlet boundary condition is and χ = 0 in Π c in (4.33a).The cell mass transfer term J n takes the form where Γ nw , Γ wn and K are all constants.The first term is motivated by the form of the corresponding mass transfer terms in Lemon et al. (2006) andO'Dea et al. (2010) with an additional dependence on the cell volume fraction so that there is no proliferation in the absence of either water (and hence solute) or cells.We have also included a dependence on the nutrient concentration via Michaelis-Menten-type kinetics.The second term represents the rate of cell death which is assumed proportional to the cell volume fraction.The reaction term R is taken to be where Γ R1 and K 1 are constant.This represents the rate of nutrient uptake by the cells, which is assumed to be of Michaelis-Menten form and proportional to the volume fractions of cells and water (for the same reasons as discussed above).
The results for this case study can be seen in Figs 2(a) and 4(a).We see in the plot of θ n that for the lowest flow rate there is slight bunching upstream of the cell population due to the higher nutrient concentration (and hence cell proliferation rate) in that region.As the flow rate is increased, the distribution of cells becomes relatively uniform along the bioreactor, due to the concentration becoming more uniform (see c plot, Fig. 2(a)).For the highest flow rate, there is significant downstream bunching of cells as the cells are also advected along the bioreactor.We also note that c is above the required minimum of 0.36 for cell viability except for the lowest flow rate Q l,in = 0.01, when the cells towards the downstream end of the bioreactor may not have a sufficient nutrient supply.We can see that the optimal flow rate Q l,opt = 7.16 × 10 −2 results in a relatively uniform distribution of cells along the bioreactor, with σ taking a minimum value of 6.7 × 10 −3 .The plot of p w shows that the reduced water pressure is monotonically decreasing along the bioreactor for all flow rates.Figure 4(a) shows that there is a slight trade-off between cell yield and distribution, with the mean μ at a maximum around Q l,in = 1, which corresponds to a greater value of σ than the minimum one.

Chemoattractant: endocrine signalling
A second situation is that of a delivered chemoattractant to which the cells respond.This is similar to endocrine signalling, where the stimulus is delivered via the bloodstream (represented by the lumen in an HFMB).This type of signalling is involved in, e.g.regulating the metabolic activity of liver and muscle cells (Berridge, 2012).In this case, we impose an inlet concentration so that c = 1 at x = 0. (5.6) We assume that cell proliferation is not limited by the chemoattractant concentration and that there is also sufficient nutrient in the system (the concentration of which we do not track) so that its concentration does not affect the growth rate.Therefore, we take where Γ nw and Γ wn are constants.The reaction term takes the form where Γ R1 is constant.This represents the rate of chemoattractant uptake by the cells, which is taken to be proportional to the concentration of chemoattractant and the volume fractions of cells and water.
The results for this case study can be seen in Figs 2(b) and 4(b).Here, we see that for low flow rates there is upstream bunching of θ n due to the negative gradient in c, which causes an increase in cells upstream due to chemotaxis.However, as the flow rate increases further and this gradient decreases, advection dominates and the cell volume fraction increases downstream once more.In this case, the (a) Fig. 2. Steady-state distributions of the cell volume fraction θ n , nutrient concentration c and reduced ECS water pressure p w versus x for a range of flow rates Q l,in including the optimal value Q l,opt .(a) Nutrient-driven proliferation case study results (from Section 5.1), end time T = 1.2, optimal flow rate Q l,opt = 7.16 × 10 −2 and corresponding standard deviation σ = 6.7 × 10 −3 .(b) Endocrine signalling case study results (from Section 5.2), end time T = 0.6, optimal flow rate Q l,opt = 0.1033 and corresponding standard deviation σ = 1.6 × 10 −3 .(c) Paracrine signalling case study results (from Section 5.3), end time T = 0.6, optimal flow rate Q l,opt = 1.0532 and corresponding standard deviation σ = 1.49× 10 −2 .Parameter values are as in Table 2. most uniform distribution of cells occurs for Q l,opt = 0.1033, with σ = 1.6 × 10 −3 .Once again, the reduced water pressure is monotonically decreasing in x for all flow rates.In Fig. 4(b), we see that this time the value of μ is fairly constant across the range of flow rates considered, but drops off slightly for higher flow rates.

Chemoattractant: paracrine signalling
The third case study models the concentration c(x, t) of a chemoattractant that is produced by the cells and also induces a chemotactic response.This local signalling between cells is known as paracrine signalling, and controls processes such as neutrophil chemotaxis during inflammatory responses, and the migration of haematopoietic stem cells in bone marrow (Berridge, 2012).Bioreactor expansion of haematopoietic stem cells in particular is needed for many therapeutic applications (see Liu et al., 2006).In this situation, we assume that no chemoattractant is supplied at the inlet and hence impose the boundary condition (5.9) We take the same mass transfer term as in Section 5.2: (5.10) Since we must also consider production of the chemoattractant here, R takes the form where Γ R2 is constant.The additional first term in (5.11) signifies production of chemoattractant by the cells at a rate proportional to the cell volume fraction and represents situations where the chemoattractant is always produced by the cells (hence higher concentrations where there are more cells, a situation which could arise, for instance, with stem cells).
Results for this case are shown in Figs 2(c) and 4(c).Overall, we can see that as the flow rate increases the extent of the downstream bunching of cells increases.However, between Q l,in = 0.1 and Q l,in = 1, the opposite trend is seen.For the lower flow rate of Q l,in = 0.1, the (positive) gradient in c is sufficiently large for chemotaxis to have a noticeable effect, resulting in an increased cell volume fraction downstream.However, once Q l,in is increased to 1, the concentration is much more uniform, and therefore, so too is the distribution of cells.When the flow rate is further increased to 10, advection dominates and the cells are swept downstream (results not shown).As in Section 5.1, the reduced water pressure is again monotonically decreasing in x for all flow rates.
Here, the optimal flow rate within the range considered is Q l,opt = 1.0532,where the standard deviation of θ n is at a minimum value of 1.49 × 10 −2 .We also see in Fig. 4(c) that the mean μ does not vary much with Q l,in in this case, and hence we would expect a similar yield of cells no matter which flow rate was chosen.However, the optimal flow rate Q l,opt does correspond to a maximum cell yield, as well as a minimum in σ .

Nutrient and chemoattractant concentrations
In the next case study, we consider two solutes: a nutrient that is supplied at the lumen inlet and a chemoattractant that is produced by the cells and induces a chemotactic response.The coupled system to solve can be derived as before to yield (where subscript 1 refers to nutrient and subscript 2 to chemoattractant concentration) (5.12) where θ n , c 1 , c 2 and p w satisfy (5.13) and the fluxes are defined by (5.14) (5.16) The corresponding boundary conditions are (5.17) (5.18) The mass transfer term is taken to be that from Section 5.1, while the reaction terms take the appropriate forms from Sections 5.1 and 5.3: (5.19) with the uptake constants Γ R1,1 and Γ R1,2 and production constant Γ R2,2 set to the same values as the corresponding parameters in the single-solute case (where Γ R1 = Γ R2 = 50).
Results for this extended system can be seen in Figs 3 and 4(d) and correspond to the case study in Section 5.3 when nutrient is explicitly modelled.Figures 2(c) and 3 show the same overall trend in cell distribution as the flow rate increases.However, in Fig. 3 there is also an overall increase in cell yield as the flow rate increases.We also note that here the cell volume fraction is lower in general.This is due to the decreased proliferation rate since we are assuming here that nutrient can have a limiting effect on cell proliferation rate, while in Section 5.3 it was assumed to be in abundance.We also see similar trends in chemoattractant concentration and reduced water pressure, albeit with higher overall values when nutrient is not explicitly modelled.This is due to the higher cell volume fraction overall which results in more chemoattractant being produced.We note that the distribution of nutrient is similar to that seen in Fig. 2(a) and that here the nutrient concentration is always above that required for cell viability.In this case, the optimal flow rate is Q l,opt = 1.341, with a corresponding standard deviation Fig. 3. Steady-state distributions of the cell volume fraction θ n , nutrient concentration c 1 , chemoattractant concentration c 2 and reduced ECS water pressure p w versus x for a range of flow rates Q l,in including the optimal value Q l,opt .Results are from the nutrient and chemoattractant case study presented in Section 5.4, end time T = 1.5 and optimal flow rate Q l,opt = 1.341 (corresponding standard deviation σ = 1.31 × 10 −2 ).Parameter values are as in Table 2.
of σ = 1.31 × 10 −2 .Figure 4 (d) shows that this flow rate produces a relatively high cell yield, and confirms that σ is at a minimum.

Parameter sensitivity
We next investigate how sensitive our model results are to the key dimensionless parameters Q e,in , η, ζ ns , ζ ws , κ, φ and λ.We note that in each case the parameter value in question is varied within a range that retains the validity of our asymptotic analysis.5.5.1 Effect of Q e,in .To examine how sensitive our model results are to the dimensionless ECS flux, we firstly calculate the mean and standard deviation of θ n for a range of values of Q e,in .Results for all case studies are summarized in Table 3 and indicate that overall, increasing Q e,in greatly increases the standard deviation σ of θ n but only slightly reduces the mean μ of θ n .Plots of the final cell, nutrient concentration and pressure distributions for the nutrient-driven proliferation case study from Section 5.1 can be seen in Fig. 5 and show that increasing Q e,in results in increased downstream bunching of cells.6(a-c)).As expected, increasing either η or ζ ns results in a lesser degree of downstream cell bunching (and therefore smaller standard deviation σ of θ n ) due to the increased cell-scaffold affinity or scaffold drag, respectively.Increasing ζ ws gives the opposite effect, due to the higher drag between the water and cells and therefore increased advection.
5.5.3Effect of membrane properties κ and φ.Next we consider the effect of the membrane porosity φ and dimensionless permeability κ on the cell population distribution.Increasing both parameters results in more downstream bunching of cells and therefore greater standard deviation σ of θ n (see Fig. 6(d) and (e)).This is due to the increased flow into the ECS from the membrane (from either a greater porosity or permeability) and therefore greater effect of advection on the cells.2.
5.5.4Effect of viscosity parameter λ.We lastly look at the effect of λ, which enters the viscosity ratio through μ w /μ n = λε so that decreasing λ increases the effective viscosity of the cells μ n .Results show that when the cells are more viscous, the distribution is less uniform overall (i.e. the standard deviation σ of θ n is greater) and also have more of a tendency to bunch downstream (see Fig. 6(f)).

Backflow
Finally, with the aim of counteracting downstream bunching effects without using a chemoattractant, we consider the influence of an opposing flow in the ECS by imposing an inlet flux at the downstream end:  2.
This corresponds to the following downstream boundary condition for the water flux We also need to impose an upstream pressure boundary condition.Unlike the standard flow regime (where atmospheric pressure is imposed), here we impose a prescribed (dimensionless) upstream pressure P u so that p w = P u at x = 0. (5.22) Experimentally, this corresponds to an additional clamp controlling the pressure at the upstream exit port.This is necessary in order to avoid a dominating upstream flow in this setup, and so to obtain a uniform cell population in the absence of chemotaxis.If a zero pressure was imposed upstream, we would always have a positive water pressure gradient in the ECS and therefore net upstream advection of the cells.The value of P u could be controlled experimentally, and here we choose illustrative values to demonstrate possible behaviours in this flow regime.2. Streamlines were calculated using the streamslice function in MATLAB.
We investigate the effect of this flow regime using the nutrient-driven proliferation case study from Section 5.1 as an example.Our aim is to see whether a uniform distribution of cells can be obtained without the effect of chemotaxis but for higher flow rates than those in Section 5.1.Hence, once again we look for flow rates that minimize the standard deviation σ of θ n .Given that the fluxes Q l,in and Q e,in should now be imposed independently of one another, we fix a value of Q l,in and allow Q e,in to vary.Plots of the cell, nutrient concentration and reduced water pressure distributions and the proliferation function J n for a range of values of Q e,in can be seen in Fig. 7. From this, we can see that as Q e,in is increased, the cell distribution switches from downstream to upstream bunching as expected, and we note that the pressure gradient changes from negative to positive.It is also still possible to obtain a fairly uniform cell population in the case where Q l,in = 2 and P u = 6 (corresponding σ = 0.0015), whereas for Q l,in = 5 and P u = 15, the profile is peaked in shape (σ = 0.0103).The plot of J n confirms that regions of lower cell volume fraction correspond to a higher net proliferation rate.
To further illustrate the behaviour of the system in the case where Q l,in = 5 and P u = 15, Fig. 8 shows plots of the cell and water streamlines, and the advective and diffusive cell flux contributions for each Q e,in value.In Fig. 8(a), we observe a net flow of water downstream as the lumen flow is dominating, resulting in the downstream bunching of cells observed in Fig. 7(b).In contrast, Fig. 8(c) shows results for Q e,in = 10, when the ECS flow dominates and the net flow (and cell bunching) is upstream.In the flux figures, we note the switch in the respective directions of the advective and diffusive cell fluxes between Fig 8(a) and (c).Finally, in Fig. 8(b), for the optimal value of Q e,in we can see that there is a balance between forward and backward water flow in the ECS, with a stagnation point at around x = 0.6 which corresponds to the peak in θ n observed in Fig. 7(b).This stagnation point is also present in the cell velocity streamlines, which in general show a flow in the direction opposite to that of the water.Furthermore, we note that, apart from near x = 0, there is a flow of water into the membrane from the ECS.

Conclusions
In conclusion, we have developed a multiphase model of fluid flow, cell population evolution and solute transport in a simplified HFMB setup.The governing equations have been reduced through careful asymptotic analysis to exploit the separation of length scales inherent in the experimental setup.The resulting leading-order system of equations has been solved numerically.The model developed thus far is generic, and results speculative due to the number of assumptions necessary at this stage.However, we believe that the work presented here is instructive in demonstrating the range of possible behaviours in such a setup and in predicting trends arising from varying flow rates and certain other key parameters.Given the early stage of experiments in this field, there are currently no data available to validate our model.However, the results presented here provide a suggested range of flow rates to be experimentally tested to see whether the cell behaviours predicted by our model are observed.We have also identified a number of experiments which we hope will be carried out in the near future (for instance, determining the fluid flux into the ECS from the lumen which is currently unknown) and which will greatly improve the parametrization of our model.Furthermore, the model in its current form could easily be applied to specific cell types and experimental setups through the choice of boundary conditions and specific kinetic parameter values.
To demonstrate the variety of possibly steady-state cell distributions, we have investigated and compared a number of different case studies which represent different experimental setups.In each case, there is a complex interaction between the effects of advection, diffusion and chemotaxis on the cell population.This has been shown to result in a variety of steady states, with cells mainly aggregated up-or down-stream, or more uniformly distributed.We have shown that in almost all cases the requirement for a sufficient nutrient supply to maintain cell function is satisfied, but for the nutrient-only case study from Section 5.1, this may not be true throughout the bioreactor for the lowest flow rate considered.When chemoattractant is delivered into the system, we observe competition between advection and chemotaxis: initial upstream bunching of cells due to the chemoattractant concentration gradient is overcome by downstream advection of cells for large enough flow rates.This competition also results in a non-monotonic relationship between cell distribution and increasing flow rate which is observed in Section 5.3, where chemoattractant is produced by the cells.By comparing the results from Section 5.3 with those in Section 5.4, we have also shown that explicitly modelling nutrient gives a similar trend in cell distribution as the flow rate is increased.In the alternative flow regime where there is a backflow in the ECS, we observe that the competition between the up-and down-stream flow rates can also result in bunching at either end of the bioreactor, or towards the centre.
We have additionally found 'optimal' flow rates by minimizing the standard deviation of the cell volume fraction and considering the corresponding cell yield.The extent to which a given setup is able to achieve a uniform cell population, and the shape of the resulting cell distribution, varies from case to case.For instance, when cells are undergoing paracrine signalling and/or nutrient-limited proliferation, the most uniform distribution possible has slight downstream bunching of the cell population.In contrast, in the case where endocrine signalling occurs, the most uniform distribution involves both up-and down-stream cell aggregation, with a dip towards the middle of the modelling domain.However, our results all show that provided an appropriate flow rate is chosen, each case study is capable of providing a cell population with relatively little variation (the largest σ being 1.49 × 10 −2 in Section 5.3) and that despite the low experimental flow rates considered here (Q l,in = 0.01 − 10 approximately corresponding to 3.86 × 10 −6 − 3.86 × 10 −3 ml min −1 ), these flow rates ensure sufficient nutrient delivery to cells in all but one case.We also note that the minimum values of σ and optimal cell distributions found in the matching case studies with and without nutrient (from Sections 5.3 and 5.4) are very similar, and hence can conclude that the trends in uniform cell distributions found here do not heavily depend upon the inclusion (or exclusion) of nutrient.
Furthermore, in most cases we either see relatively little variation in the mean value μ of θ n , or find that minimizing the standard deviation σ of θ n gives a value of μ that is very close to its maximum, and hence by choosing optimal flow rates based on cell distribution, we are not compromising on cell yield.In the nutrient-only case (Section 5.1), however, the optimal flow rate gives a lower cell yield than could be achieved with higher flow rates, and hence in this situation a trade-off may have to be made between desired cell yield and distribution.We note that the framework outlined here could also be used to find optimal conditions in the case where a non-uniform distribution of cells is required.
Finally, we have used the nutrient-limited growth case study from Section 5.1 to investigate the sensitivity of our results to a number of key dimensionless parameters and the alternative flow regime where there is an upstream flux into the ECS.For instance, we have found that increasing the flux into the ECS does not significantly affect cell yield, but does decrease the uniformity of the cell distribution by causing more downstream bunching.In the case of backflow, we have shown that by choosing an appropriate upstream lumen flow rate and downstream ECS flow rate, it is possible to generate a uniform distribution of cells in this alternative flow regime, without the need for chemotaxis and with a higher flow rate than in the standard flow regime.
There are a number of assumptions or simplifications in our current model that should be explored further.Firstly, we greatly simplified the bioreactor geometry in order to obtain our 2D rectangular modelling region.The corresponding axisymmetric problem and comparisons with results presented here will be a focus of future work.In addition, solving the full problem numerically (in the whole bioreactor region including the exit ports) would give great insight into the validity of the model reduction presented here.However, this would be very difficult numerically and is beyond the scope of this paper.Secondly, we assume that the ECS scaffold is rigid and inert, neglecting any degradation or ECM deposition by the cells.This simplification allows us to make analytical progress which would not be possible if all the three phases were active.However, we hope to include the production and decay of scaffold in future work.In addition, we have assumed that the effective viscosity of the bulk cell phase is much greater than that of the water.Although this seems sensible, the validity of this assumption should be investigated if possible.
As mentioned previously, the constitutive forms of the kinetic functions J n and R have been chosen to demonstrate the range of behaviour that could be expected from this system, and would need validation in order for this model to be applied to a specific experimental setup.However, given that our analysis was carried out for general J n and R, it would be straightforward to obtain new results for different kinetic functions.In addition, the kinetic parameter values currently used are not specific to any one cell type and, again, the model could easily be altered for specific cell types, given appropriate data.It is noted that these alterations could quantitatively alter our findings (for instance, proliferation could become the dominating effect over advection and therefore lead to more uniform cell distributions for higher flow rates than observed here), but we would expect results to correlate qualitatively to one of the case studies presented here.
.20a-d)    and so our assumption on the size of κ induces a pressure gradient at leading order across the porous membrane.Boundary conditions (4.13c-e) yield p l = p m , c l = c m ,

Fig. 4 .
Fig. 4. Plots of the mean μ of θ n (solid line) and the standard deviation σ of θ n (dash-dotted line) versus Q l,in for case study results from (a) nutrient-driven proliferation (Section 5.1), (b) endocrine signalling (Section 5.2), (c) paracrine signalling (Section 5.3) and (d) nutrient and chemoattractant concentrations (Section 5.4).The dashed vertical line in each plot corresponds to the optimal flow rate Q l,opt for which σ is at its minimum value.Parameter values are as in Table2.

Fig. 5 .
Fig.5.Steady-state distributions of the cell volume fraction θ n , nutrient concentration c and reduced ECS water pressure p w versus x from the nutrient-driven proliferation case study (Section 5.1) for a range of values of the ECS flux Q e,in and fixed lumen flux Q l,in = 1 (end time T = 0.9).All other parameter values are as in Table2.

θFig. 6 .
Fig.6.Steady-state distributions of the cell volume fraction θ n versus x from the nutrient-driven proliferation case study in Section 5.1 (end time T = 0.9) for a range of values of (a) the cell-scaffold affinity parameter η, (b) the cell-scaffold drag ζ ns , (c) the cell-water drag ζ ws , (d) the dimensionless membrane permeability κ, (e) the membrane porosity φ and (f) the viscosity parameter λ.The (dimensional) lumen flux Q l,in is fixed at 1.02 × 10 −9 m 2 s −1 .All other parameter values are as in Table2.

Fig. 8 .
Fig. 8. Plots of streamlines for the cell (left) and water (centre) velocities, and advective and diffusive cell fluxes (right), for the nutrient-driven proliferation case study (Section 5.1), with fixed Q l,in = 5 and P u = 15.In (a) Q e,in = 0.01, in (b) Q e,in = Q e,opt = 5.8443 and in (c) Q e,in = 10.All other parameter values are as in Table2.Streamlines were calculated using the streamslice function in MATLAB.

Table 1
Dimensional parameters.The different units for the cell death rate coefficient Γ wn (which appears in J n in Section 5) depend on which case study is under consideration a Values based on estimations by our experimental collaborators.

Table 3
Values of the mean μ of θ n and the standard deviation σ of θ n for different values of the ECS flux Q e,in , and percentage change overall for each case study.The lumen flux Q l,in is fixed at 1 and all other parameter values are as in Table2