Mathematical modelling of cell layer growth in a hollow ﬁ bre bioreactor ☆

Generating autologous tissue grafts of a clinically useful volume requires e ﬃ cient and controlled expansion of cell populations harvested from patients. Hollow ﬁ bre bioreactors show promise as cell expansion devices, owing to their potential for scale-up. However, further research is required to establish how to specify appropriate hollow ﬁ bre bioreactor operating conditions for expanding di ﬀ erent cell types. In this study we develop a simple model for the growth of a cell layer seeded on the outer surface of a single ﬁ bre in a perfused hollow ﬁ bre bioreactor. Nutrient-rich culture medium is pumped through the ﬁ bre lumen and leaves the bioreactor via the lumen outlet or passes through the porous ﬁ bre walls and cell layer, and out via ports on the outer wall of the extra-capillary space. Stokes and Darcy equations for ﬂ uid ﬂ ow in the ﬁ bre lumen, ﬁ bre wall, cell layer and extra-capillary space are coupled to reaction – advection – di ﬀ usion equations for oxygen and lactate transport through the bioreactor, and to a simple growth law for the evolution of the free boundary of the cell layer. Cells at the free boundary are assumed to proliferate at a rate that increases with the local oxygen concentration, and to die and detach from the layer if the local ﬂ uid shear stress or lactate concentration exceed critical thresholds. We use the model to predict operating conditions that maximise the cell layer growth for di ﬀ erent cell types. In particular, we predict the optimal ﬂ ow rate of culture medium into the ﬁ bre lumen and ﬂ uid pressure imposed at the lumen outlet for cell types with di ﬀ erent oxygen demands and ﬂ uid shear stress tolerances, and compare the growth of the cell layer when the exit ports on the outside of the bioreactor are open with that when they are closed. Model simulations reveal that increasing the inlet ﬂ ow rate and outlet ﬂ uid pressure increases oxygen delivery to the cell layer and, therefore, the growth rate of cells that are tolerant to high shear stresses, but may be detrimental for shear-sensitive cells. The cell layer growth rate is predicted to increase, and be less sensitive to the lactate tolerance of the cells, when the exit ports are opened, as the radial ﬂ ow through the bioreactor is enhanced and the lactate produced by the cells cleared more rapidly from the cell layer.


Introduction
The aim of in vitro tissue engineering is to produce cells and tissues in the laboratory that can be used to replace or repair damaged or lost tissues in a patient's body. Generating these cells and tissues from the patient's own cells (autologous cells) has several advantages, including decreased likelihood of immune rejection, but requires expansion of the original cell population taken from the patient. This can be achieved by seeding the cells onto a biomaterial scaffold and incubating the cell-scaffold construct in a bioreactor.
Hollow fibre bioreactors (HFBs) show great promise as cell expan-sion devices. They consist of a cylindrical glass module housing a single or multiple porous, hollow, biodegradable polymer fibres. Nutrient-rich culture medium is pumped through the fibre lumen(s) and forced through the fibre wall(s) (membranes) to cells seeded in the surrounding space (the extra-capillary space or ECS). There are ports at either end of the ECS, which may be opened to promote radial flow through the bioreactor, or left closed (see Fig. 1). With the ECS ports open, and hence at atmospheric pressure, the flow through the membrane is controlled by fixing the fluid pressure at the downstream lumen outlet, and this allows the nutrient delivery and fluid shear stress experienced by the cells in the ECS to be controlled. Cells can either be seeded directly onto the outer surface of the fibres, or in a hydrogel throughout the ECS (Wung et al., 2014). As well as enabling the chemical and mechanical environment of the cells to be controlled, HFBs provide a highly efficient means of expanding cells due to the large surface area of the fibres available for cell proliferation relative to the bioreactor volume-as many cells can be cultured in a 0.5 l HFB as in a 1000 l standard culture flask (Wung et al., 2014). Despite these advantages, there is still potential to improve operation and scale up of HFBs for clinical use. In particular, better understanding of the combined effects of cell seeding and operating conditions on cell population growth in HFBs is required. Mathematical modelling of cell growth in HFBs can help significantly in this regard by predicting parameter ranges for which growth is maximised, and thereby streamlining the choice of parameters and operating conditions to test experimentally, saving time and money. In this study we focus on the relatively simple single-fibre HFB shown in Fig. 1, to gain deeper insight into the factors that affect growth.
Experiments with the single-fibre HFB have shown that cells seeded onto the outer surface of the fibre will proliferate over the surface until they reach confluence given the right culture conditions and culture period (Ellis and Chaudhuri, 2007;Meneghello, 2010). Additionally, once cells have reached confluence, they can proliferate outwards into the ECS to establish a tissue layer multiple cells deep (Tharakan and Chau, 1986;Lu et al., 2005;Ye et al., 2007;De Napoli et al., 2011). Although oxygen concentration, lactate concentration and fluid shear stress are known to affect the rate of cell proliferation and death, how they interact and affect the rate and extent of cell layer growth into the ECS have not yet been investigated. This is the primary purpose of our study, which will be achieved through mathematical modelling.
Previous modelling studies of HFBs have focussed on describing fluid and solute transport through the bioreactor (Kelsey et al., 1990;Piret and Cooney, 1991;Sullivan et al., 2006;Shipley et al., 2010Shipley et al., , 2011Ye et al., 2006) and considered the timescales associated with these transport processes rather than the longer timescale associated with cell proliferation (see Shipley et al., 2010Shipley et al., , 2011Brotherton and Chau, 1996 for more detailed reviews). Hence the effects of cell proliferation on fluid flow and solute transport in the bioreactor have generally been ignored. Most transport models for HFBs are based on the Krogh cylinder approximation (Brotherton and Chau, 1996). In Krogh cylinder models, the flow and mass transport in the fibre lumen are modelled using the Stokes flow equations and advection-diffusion equations respectively, and convection effects in the membrane and ECS are usually ignored, which is representative of the HFB set-up in Fig. 2(a) with the ECS ports closed. When flow through the membrane and ECS has been considered (which is relevant when one or more ECS ports are open), it has been modelled using Darcy's law, and nutrient uptake and waste product synthesis have been described by adding appropriate reaction terms to the ECS mass transport equations. In terms of nutrient transport and consumption, oxygen is the most widely modelled solute (Schonberg and Belfort, 1987;Piret and Cooney, 1991;Patzer II, 2004;Sullivan et al., 2006;Davidson et al., 2010;Shipley et al., 2011;Pearson et al., 2014), but glucose (Ye et al., 2006;Abdullah et al., 2009;Das, 2007) and proteins (Labecki et al., 1996(Labecki et al., , 2004Shipley et al., 2009) have also been considered. By contrast, the transport of lactate, a waste product of cell metabolism toxic to cells in high concentrations, has only been modelled recently (Shipley and Waters, 2012). Shipley and Waters (2012) developed a model of fluid flow and oxygen and lactate transport in a single-fibre HFB with a cell-packed ECS in the absence of growth, and predicted that improving oxygen delivery to, and lactate removal from, the cells by opening the ECS ports would enable a larger cell population to be cultured in the bioreactor. Here we extend their model to allow for a non-uniform growing layer of cells attached to the outer surface of the fibre in the ECS. We assume that the rate of growth of the cell layer depends on the local oxygen and lactate concentrations and the fluid shear stress experienced by the cells. We consider both the flow configuration in which the ECS ports are closed (no flow in the membrane or ECS) and that in which they are open (flow throughout the bioreactor). The lumen outlet pressure controls the ratio of the flow through the membrane to that through the lumen. Although promoting radial flow through the bioreactor improves oxygen delivery and lactate removal to and from the cells, it also increases the shear stresses that the cells are exposed to, which can lead to cell death and detachment from the cell layer. Incorporating the effect of excess shear stress on the cell layer growth allows us to predict the optimal lumen inlet flow rate and lumen outlet pressure to maximise the growth of the cell layer for cells with different nutrient demands and sensitivity to shear stress. Modelling lactate transport and cell death due to excess lactate allows us to investigate the sensitivity of the cell layer growth to the lactate tolerance of the cells.
In addition to being the first model for a freely growing cell layer in a HFB, the work presented here differs from most previous theoretical studies of cell culture in HFBs in two key respects. First it accounts for the potentially negative effects on growth of high shear stress from high lumen inlet flow rates and lumen outlet pressures. Second it considers the feedback effect of tissue growth on fluid flow and solute transport in the bioreactor. The only other studies that have considered these effects are the multiphase models of cell culture in the single-fibre HFB of Pearson et al. (2014Pearson et al. ( , 2015a and our previous study on cell population expansion in the single-fibre HFB (Chapman et al., 2014). Multiphase models are continuum models that account for interactions between the flow, solute transport and cells by treating the cells and culture medium as separate phases with their own time-and spacedependent volume fractions. Mass transfer between the phases is described using constitutive source/sink terms in the mass conservation equations for each phase. Pearson et al. (2015b) used a 2D multiphase framework to describe shear-stress-dependent proliferation in a cell layer of constant depth in a HFB. They used the model to predict the effects of altering the flow rate into the ECS and the cell layer depth on the cell yield. We take a different approach here, by assuming that there is no movement of cells within the cell layer and that cell proliferation and death are localised to its outer surface, where they induce growth/recession of the cell layer, with the cell volume fraction in the cell layer remaining fixed. In our previous study (Chapman et al., 2014), we developed a 2D model of oxygen-and shear-stress-dependent cell aggregate growth along the outer surface of the fibre in a single-fibre HFB, applicable to the initial stages of cell culture, and used it to predict the lumen inlet flow rate, lumen outlet pressure and initial seeding distribution that minimised the time taken for the aggregates to reach confluence over the fibre surface. The model developed here is applicable to a later stage of cell culture, after the cells have reached confluence, when they are proliferating out into the ECS. The modelling framework presented here is similar to that used in Chapman et al. (2014), with equations for fluid and solute transport through the bioreactor coupled to equations representing cell population growth. However, here we consider an axisymmetric cylindrical geometry with a distinct cell layer rather than a simplified 2D geometry with cell aggregates, and employ additional equations to describe the fluid flow and solute transport through the cell layer as well as its growth. In addition to using the model to predict the optimal lumen inlet flow rate and lumen outlet pressure to maximise the growth of the cell layer, we also briefly explore the potential impact of non-uniformity in the initial cell layer depth on its subsequent growth.
The paper is organised as follows. The model set-up and governing equations for the system with the ECS ports open are presented in Section 2. Model parameterisation, reduction and solution are summarised in Section 3 (further details can be found in the Appendices). In Section 4 we compare the fluid flow and oxygen and lactate distributions for static cell layers of uniform depth and non-uniform depth, and describe how the flow and solute distributions change over time for a growing cell layer. We then present the results of simulations of cell layer growth for different cell types, including predictions of the optimal lumen inlet flow rate and lumen outlet pressure for growth. We summarise the key differences in the flow, solute transport and growth when the ECS ports are closed in Section 5. In Section 6 we discuss our key findings and suggest possible model extensions.

Model set-up
Fig. 2(a) shows a schematic of the single-fibre HFB set-up with a cell layer surrounding the fibre. We assume that the lumen is cylindrical with a circular cross-section and that the fibre, cell layer and ECS are coaxial and annular in cross-section. For simplicity, we view the flow and solute distributions as symmetric about the lumen axis, and treat the ECS ports as being distributed over the whole of the outer wall of the ECS following Shipley et al. (2010) and Shipley and Waters (2012) (Fig. 2(b)). The system can therefore be described using axisymmetric cylindrical polar coordinates (r, z), where r is the radial distance from the lumen axis and z is the distance along the lumen axis from the inlet at z=0 to the outlet at z=L. The corresponding unit vectors in the r-and z-directions are denoted by e r and e z . The radius of the lumen and outer radii of the fibre and ECS are denoted by R l , R m and R e .
Culture medium is pumped into the lumen inlet at a prescribed volumetric flow rate Q l in , and a pressure P l out , is imposed on the fluid at the lumen outlet. When the ECS ports are open, P l out , controls the ratio of the flow through the membrane to that down the lumen. Oxygen enters the system at constant concentration C in via the culture medium pumped into the lumen inlet, is transported through the bioreactor by a combination of advection and diffusion, and is consumed by the cells in the cell layer. Lactate is produced by the cells in the cell layer as a by-product of respiration and transported out of the system by advection and diffusion. While high oxygen levels promote proliferation, excess lactate causes cell death.
Cell proliferation and death are assumed to be localised to the outer surface of the cell layer, and to occur at a rate dependent on the local oxygen and lactate concentrations and interstitial fluid shear stress. As the cell layer grows or recedes, the flow, oxygen uptake and lactate production change, and this in turn affects the growth. Thus, the radial position of the outer surface of the cell layer, R c , varies with axial position and time, i.e. R R z t = ( , ) . The aspect ratio of the fibre lumen, ϵ, is therefore very small as are those of the membrane, cell layer and ECS. We exploit this fact to simplify the model of the fluid flow, solute transport and cell layer growth given in Sections 2.2.1-2.2.3, so that we can make analytical progress in solving the model.

Governing equations
The model for the cell layer growth is composed of equations describing the fluid flow through the bioreactor (Section 2.2.1), equations for the transport and consumption/production of oxygen and lactate overlaid on the fluid transport model (Section 2.2.2), and equations for the growth of the cell layer coupled to the fluid and solute transport equations (Section 2.2.3). All values for the parameters that appear in the following equations, the dimensional analysis of the model, and an explanation of how parameters unobtainable from the literature were estimated are given in Section 3, and Appendices A and B.

Fluid transport
The flow in the lumen and ECS can be modelled as steady, incompressible Stokes flow: where u i (i l e = , ) and p i (i l e = , ) are the fluid velocity and pressure in the lumen and ECS (denoted by subscripts l and e), and μ is the dynamic fluid viscosity.
Following several other authors (Brotherton and Chau, 1996;Ye et al., 2006;Abdullah and Das, 2007;Shipley and Waters, 2012), we model the membrane and cell layer as rigid porous media and describe the fluid flow through them using the incompressible Darcy flow equations where subscripts m and c denote the membrane and cell layer, ϕ i (i m c = , ) are the membrane and cell layer porosities (assumed constant since the cell volume fraction in the cell layer is taken to be fixed), u i and p i (i m c = , ) are the interstitial fluid velocities and pressures (averaged over the fluid volume), and k i (i m c = , ) are the fluid permeabilities of the membrane and cell layer.
We now give the boundary conditions for the fluid transport: first those on the interfaces between the lumen, membrane, cell layer and ECS and on the outer wall of the ECS, and then the axial boundary conditions at the ends of each region. For these boundary conditions, we require the definition of the fluid stress tensors in the different regions, σ i (i l m c e = , , , ), corresponding to Eqs.
(2) and (3) At the lumen-membrane and membrane-cell layer interfaces we prescribe continuity of normal fluid velocity and normal stress. Following previous studies (Shipley and Waters, 2012;Pearson et al., 2014Pearson et al., , 2015b, we assume that the stress is transmitted via the fluid phase, so that the continuity conditions at the lumen-membrane and membrane-cell layer interfaces are σ σ ϕ r R u e u e e e e e · = · , · · = · · on = , l r m m r r l r r m r l (6) σ σ ϕ ϕ r R u e u e e e e e · = · , · · = · · on = .
m m r c c r r m r r c r m By (5), the normal stress condition on the membrane-cell layer interface is equivalent to continuity of pressure On the moving outer boundary of the cell layer, r R z t = ( , ) c , we impose continuity of the normal fluid velocity relative to the moving boundary Flow of the culture medium past the permeable boundary of the membrane will cause a boundary layer to develop on the membrane surface in which the tangential fluid velocity is non-vanishing (Beavers and Joseph, 1967). This can be modelled using a Beavers-Joseph boundary condition (Shipley et al., 2010). However, following Shipley et al. (2010), who showed that the slip at the permeable boundary of the PLGA-PVA membrane has a negligible impact on the flow, we apply a no-slip boundary condition for the tangential component of the fluid velocity on r R = l ϕ r R u e u e · = · on = .
l z m m z l Similarly, we assume that the slip at the outer surface of the porous cell layer is negligible, and impose no slip there where R z R z t e e = (∂ /∂ + )/ 1 + (∂ /∂ ) c c r z c 2 is the tangent vector to the surface. To test the validity of this assumption the magnitude of the slip at the outer surface of the cell layer should be determined from measurements of the flow distribution with a cell layer attached to the hollow fibre. However, in the absence of such experimental data, we make the simplifying assumption that the cell layer permeability is low enough for (11) to provide a reasonable approximation to the Beavers-Joseph boundary condition.
Typical axial lumen flow velocity 0.013-0.13 m s −1 - Lumen outlet pressure 1.027 × 10 -2.068 × 10 Pa We treat the ECS ports as being distributed over the entire curved wall of the ECS (r R z L = , ∈ [0, ] e ), henceforth referred to as 'the ECS port', and set the fluid pressure to atmospheric pressure and the axial velocity to zero on this boundary p P r R u e = , · = 0 on = .
e atm e z e (12) Shipley et al. (2010) confirmed that this approach gives excellent agreement between model predictions and experimental measurements of the lumen and ECS flow rates when there is no cell layer, as the bulk of the pressure drop from the lumen to the ECS occurs across the membrane. Whilst the precise details of the flow in the cell layer and ECS will determine how closely the model predictions agree with experimental flow rate measurements when there is a growing cell layer, the membrane permeability will still be the key determinant of the flow in the ECS, since it is much lower than the cell layer permeability (see Table 1). Treating the ECS port as a distributed port should therefore still give a good approximation to explicitly modelling flow through the two ECS ports. The ends of the fibre are glued into place with epoxy resin and the walls of the ECS are solid, so there is no flow out of the ends of the membrane, cell layer and ECS. Hence, In the experimental set-up, the volumetric flow rate of fluid into the lumen and normal stress at the lumen outlet are fixed, so These conditions are sufficient to determine the fluid pressures and velocities in the reduced model for the fluid flow obtained by exploiting the small aspect ratio of the bioreactor (see Appendices A and C).

Mass transport: oxygen and lactate
As the timescales for advection, diffusion and consumption of oxygen and production of lactate are much shorter than the timescale for cell layer growth for the flow rates we consider (see Section 3), we assume that the solute transport is quasi-steady on the growth timescale, and thus governed by steady advection-diffusion equations in the lumen and ECS where c i (i l e = , ) is the solute (oxygen or lactate) concentration and D i (i l e = , ) are the assumed constant solute diffusivities in the lumen and ECS. In the membrane, the transport is governed by where c m is the solute concentration (per unit volume of fluid) and D m is the effective solute diffusivity (accounting for dispersion effects). In the cell layer, there is also oxygen uptake and lactate production, so the solute transport is described by the reaction-advection-diffusion equation where c c is the solute concentration and D c is the effective solute diffusivity, and the reaction term mϕ c ( ) c c describes the rate of oxygen uptake (for m = −1) or lactate production (for m=1) by the cells. Eqs. (17) and (18) are derived by volume averaging (reaction-)advectiondiffusion equations for the solute concentrations in the pore space over the membrane and cell layer, respectively (see Gray, 1975;Quintard and Whitaker, 1994 for further details).
Oxygen uptake by cell populations in HFBs is typically modelled by Michaelis-Menten kinetics (Abdullah and Das, 2007;Chen and Palmer, 2009;Das, 2007;Pillarella and Zydney, 1990;Shipley et al., 2011), for which where V max (in mol m s −3 −1 ) is the maximal uptake rate per unit volume of the cell layer and C 1/2 is the concentration at which the uptake rate is half-maximal. Here we restrict attention to cell types for which c C ≫ c 1/2 (see Section 3), so that (19) can be approximated as constant and (18) can be solved analytically for a given position of the outer boundary of the cell layer, r R z t = ( , ) c . We also assume that the lactate production rate is approximately constant, since the culture medium is sufficiently glucose-rich and the lactate concentration kept sufficiently low by lactate buffering for this to be a reasonable approximation (see Shipley and Waters, 2012).
Since the oxygen and lactate distributions are assumed to be axisymmetric, we impose no diffusive flux through r=0 At the lumen-membrane and membrane-cell layer interfaces we impose continuity of solute concentration and flux At the outer surface of the cell layer we impose continuity of the concentration and flux relative to the moving boundary With the normal velocity continuity conditions in (6), (7) and (9), the flux conditions above reduce to continuity of diffusive flux There is no concentration flux out of the ends of the membrane, celllayer and ECS, so Experimentally, the oxygen concentration at the lumen inlet is held constant and there is no lactate in the culture medium entering the lumen, so Following Shipley and Waters (2012), we assume that the culture medium leaving the lumen outlet and ECS port is well-mixed and impose zero-diffusive-flux at these outlets for both the oxygen and lactate transport

Cell layer growth
We assume that the cell density in the cell layer is constant and that changes in the cell number occur at its outer surface, which moves outwards if the cell proliferation rate exceeds the cell death rate and inwards if the death rate exceeds the proliferation rate (due to cells that die detaching from the surface and being carried away by the flow). Cells in the bulk of the cell layer are assumed to be stationary and quiescent due to contact inhibition (Engel et al., 2005;Machide et al., 2006;Lee et al., 2003;Guo et al., 1989).
Cell proliferation and death are assumed to depend on the oxygen and lactate concentrations, c and l, and interstitial fluid shear stress σ (the shear stress on the cells due to the flow through the spaces between the cells) at the surface of the cell layer. We assume that the growth depends on the interstitial shear stress rather than the shear stress on the cells from the ECS flow tangential to the surface, as the interstitial shear stress is an order of magnitude larger for the flow configuration considered here, both when the ECS ports are open and when they are closed. Thus the normal velocity of the outer surface of the cell layer is given by where the growth function G is to be prescribed. We note that this model is equivalent to the thin-rim (fast-consumption) limit of a twophase free-boundary model for the cells and culture medium in which fast nutrient consumption at the tissue surface prevents nutrient transport to the tissue interior so that cell proliferation is consigned to the free boundary (King and Franks, 2006), except that in the present model growth in the cell layer interior is prevented by contact inhibition rather than lack of nutrients. We estimate the interstitial shear stress from the interstitial fluid velocity by assuming that the flow through the interstitial spaces can be approximated as Poiseuille flow through a circular duct of diameter d with mean velocity u | | c . With r ∼ as the local radial coordinate, the interstitial velocity profile is therefore so the interstitial shear stress on the cells is Following previous studies (McElwain and Ponzo, 1977;Jones et al., 2000;Lewis et al., 2005), and based on experimental evidence that the proliferation rate of many cell types increases with oxygen availability Kim et al., 2001;Murrell et al., 1990), we assume that the cell proliferation rate increases linearly with the oxygen concentration provided it exceeds the minimum threshold required for proliferation, C min . As data on the relationships between cell death and lactate concentration, and cell death and shear stress, is more limited, we assume that cells die and detach from the cell layer surface at a constant rate when either the lactate concentration exceeds the threshold at which it is toxic to the cells, L max , and/or the shear stress exceeds a critical threshold, Σ d . The effects of excess lactate and shear stress are taken to be multiplicative, so that cells die at a faster rate when both l L > max and σ Σ > d . We use a smoothed Heaviside function rather than a normal Heaviside function to describe the lactate and shear stress effects as l and σ cross L max and Σ d , as a normal Heaviside function would lead to discontinuities in R t ∂ /∂ c and very sharp changes in R z t ( , ) c that would be difficult to resolve numerically and for which the reduced model derived in the small aspect ratio limit (Appendix A) would not be valid. Since the cell layer stops growing if it fills the entire ECS or recedes to the fibre surface (which is taken to represent all the cells dying), Thus we arrive at the following form for G where H(·) is the Heaviside function and x s sx F( , ) = (1 + tanh( )) 1 2 is the smoothed Heaviside function with smoothing factor s; the constant A p is the growth rate of the cell layer per unit concentration if c C > min and there is no cell death (l L < max and σ Σ < d ); and B d , B L and B s together determine the recession rate of the cell layer if l L > max and/or σ Σ > d and there is no cell proliferation (c C < min or σ Σ > d ). Finally, we prescribe the initial depth of the cell layer as a function of z The full model for the flow, solute transport and cell layer growth is thus given by Eqs.

Model parameterisation and reduction
Typical values of the model parameters are given in Tables 1-3. Where parameters are not obtainable from the literature, we determine them via physical arguments in Appendix B. We nondimensionalise Eqs.
(2)-(36) and exploit the small aspect ratio of the bioreactor to simplify the system (see Appendix A for details). For the inlet flow rates in Table 1, the reduced Reynolds number for the flow (the ratio of inertial to viscous forces in the fluid), Re ρUL μ ϵ =ϵ / 2 2 , is in the range 5.3 × 10 -5.3 × 10 −3 −2 (i.e. viscous forces in the fluid dominate inertial forces), so we are justified in neglecting inertia in the lumen and ECS flow equations (2). The reduced Péclet numbers in the different regions (the ratios of the rates of axial advection to those of radial diffusion), (i l m c e = , , , ), are all (1) , so advection and diffusion are equally important. Although the maximum oxygen uptake/lactate production rate, V max , is cell-type dependent (see Table 2), we assume that the dimensionless uptake/production rate is (1) to retain a balance between the rates of uptake/production and diffusion in the reduced model, thus keeping our analysis as general as possible. We estimate the cell layer growth rate per unit concentration, A p , for different cell types (Table 3) from their cell doubling times (see Table 2 Cell culture, oxygen uptake, and lactate production parameter values for different cell types.

Solute
Cell type Cell density Appendix B.2). Since the cell layer growth timescale is much longer than the timescales for advection and solute diffusion, we are justified in assuming that the flow and solute transport are quasi-steady on the growth timescale (see Appendix B.2). In our simulations we vary seven parameters to determine their effect on the cell layer growth: the lumen inlet flow rate Q l in , , the lumen outlet pressure P l out , , the maximum oxygen uptake/lactate production rate, V max , the growth rate A p , and the oxygen and lactate concentration thresholds and shear stress threshold for cell proliferation and death, C min , L max and Σ d . In the absence of experimental data for the concentration at which lactate becomes toxic to the cells, L max , we follow Whittaker et al. (2009)

Uniform vs non-uniform cell layer depth
Before simulating the cell layer growth, we compared the flow, and oxygen and lactate distributions in layers of rat cardiomyocytes of fixed uniform and non-uniform depths to assess the impact of non-uniformity in the depth on oxygen delivery to, and lactate removal from, the cells. Variations in the cell layer depth can arise experimentally from non-uniformity in the initial cell seeding on the surface of the fibre or in the subsequent cell proliferation, or both. Fig. E1 in Appendix E shows that the uniformity of the cell layer depth has little effect on the fluid pressure distribution, and therefore the flow profiles (Fig. E2), in the bioreactor. This is because the majority of the pressure drop from the lumen to the ECS occurs across the membrane and not the more permeable cell layer. However, the uniformity of the cell layer can have a significant effect on the oxygen and lactate distributions. For example, changing from a constant cell layer depth to one that decreases linearly along the fibre (but keeping the volume of the cell layer the same), switches the minimum oxygen and maximum lactate concentrations from being at the ECS wall around the lumen outlet  E4(b)). This is because more oxygen is taken up and more lactate produced where the cell layer is thicker, since the number of cells in the cross-section of the cell layer is greater (as the cell density is assumed constant). The sensitivity of the solute distributions to the uniformity of the cell layer depth suggests that fine control over the culture conditions may be required to maintain uniform solute distributions, and hence uniform growth, for cells with high nutrient demands.

Changes in flow profiles and solute distributions with growth
The cell layer grows outwards from the fibre provided that the shear stress is below the threshold at which cells die ) and the positive and negative effects of oxygen and lactate on the growth balance each other, then the cell layer attains a steady state depth. Since the radial flow velocity in the ECS decreases with r (as shown in Fig. E2b), the shear stress on the cells at the outer surface decreases as the layer grows. Hence, if σ| r R = c is initially below Σ d , it will remain so as the cell layer grows (see Fig. 4(c)). Conversely, if σ| r R = c is initially above Σ d , then the cell layer recedes and the shear stress on the surface cells increases. Hence, the layer will continue to recede until it reaches the fibre surface (see Figs. 5 and 6, which are discussed in detail in Section 4.2.1).

Cell layer growth simulations
For sustained growth of the cell layer it is required that . In this investigation we assume that the bioreactor dimensions and membrane properties (its depth and permeability) are fixed, and that our control parameters are Q l in , and P l out , and the initial depth of the cell layer , as these are the parameters that can be readily prescribed experimentally. As the governing equations for the flow, solute transport and cell layer growth must be solved numerically, it is not possible to define analytical operating equations for Q l in . We consider two cell types: rat cardiomyocytes, which have high oxygen requirements (max- . The simulations enable us to predict values of Q l in , and P l out , that maximise the average depth to which the cell layer grows in a set time or minimise the time it takes to grow to fill the ECS. We also briefly consider the impact of variation in the initial cell layer depth with z on the cell layer growth, and the sensitivity of the growth to the lactate toxicity threshold, L max . The simulations are run for a maximum period of T=60 days, or until a stopping criterion is reached (the cell layer grows to fill the entire ECS, recedes to the fibre surface or reaches a steady state). A maximum culture time of 60 days is chosen since significant degradation of the PLGA in the membrane occurs for longer times (Hoque et al., 2012;Azimi et al., 2014), increasing the membrane permeability and invalidating our assumption that it remains constant. Once 60 days have elapsed, the mean and standard deviation of the final outer radius (shown by the green circle in Fig. 3(a)), and decreases sharply when P > 1.965 × 10 Pa l out , 5 (to the right of the vertical green line in Fig. 3(a)). This change is due to the changes in the balance of the effects of increased oxygen concentration and increased shear stress at the surface of the cell layer as Q l in , and P l out , increase (l L ≪ max for both cell types over the range of flow rates and outlet pressures considered). For P < 1.965 × 10 Pa l out , 5 , the shear stress σ is below the cell death threshold, Σ = 0.03 Pa d , over the whole length of the fibre, so no cell death occurs (Fig. 4), and cell proliferation increases as Q l in , and P l out , increase due to improved oxygen delivery. For P > 1.965 × 10 Pa l out , 5 , as P l out , increases σ exceeds Σ d over an increasing length of the fibre, so cells die and the cell layer recedes towards the fibre surface at r = 400 μm (Fig. 5). For chondrocytes, R c increases monotonically as Q l in , and P l out , increase, since the shear stress is much lower than the threshold at which they die, Σ = 2 Pa d .
Our model predicts that the growth rate of the cardiomyocyte layer is maximised if we fix the flow rate and outlet pressure to be Q = 3.33 × 10 m s ; for the chondrocytes the growth rate is maximised when Q = 3.33 × 10 m s (green circle in Fig. 3(b)). For the cardiomyocytes, the standard deviation in the outer radius of the cell layer at t=60 days, SD, decreases as Q l in , decreases and as P l out , increases for P < 1.931 × 10 Pa l out , 5 , ranging from 0.2 μm (when Q = 2 × 10 m s ). For the chondrocytes, SD also decreases as Q l in , decreases and P l out , increases, but is much smaller than for the cardiomyocytes for each combination of Q l in , and P l out , , SD = 0-2 μm. The greater variation in the depth of the cell layer for the cardiomyocytes is due to their greater oxygen uptake rate. Their greater uptake rate means that at low outlet pressures increasing the flow rate only increases oxygen delivery to the upstream end of the cell layer, so it grows more quickly than the downstream end.

Impact of variation in initial cell layer depth on growth
In the previous simulations the initial depths of the cell layers were uniform along the fibre. However, since the oxygen and lactate distributions are affected by the degree of uniformity of the cell layer, variation in its initial depth along the fibre can significantly affect its subsequent growth. Such variation in the cell layer depth often arises as a result of spatial heterogeneity in the cell seeding distribution, so in this section we explore its potential impact on the growth of the cell layer.
We simulate the growth of the layer of rat cardiomyocytes described in Section 4.1.1, whose depth decreases linearly along the fibre (as a simple example of a non-uniform cell layer), at the flow rate and outlet pressure predicted to be optimal for the growth of the uniform layer in ). We obtain a very different result from Fig. 4, as shown in Fig. 6. Since the non-uniform cell layer is initially thinner at its downstream end, the shear stress on the cells at the surface is higher than for the uniform layer and exceeds Σ d . Hence the cell layer quickly recedes towards the fibre surface at its downstream end (for z > 7.38 cm), but grows outwards over the upstream part of the fibre (z < 7.36 cm) (Fig. 6(d)), where σ Σ < d (Fig. 6(c)). The shear stress on the cells at the downstream end increases as the cell layer recedes, so the layer continues to recede until it reaches the membrane at r = 400 μm after 5 days. Consequently, after 60 days the mean final outer radius of the cell layer, R = 750 μm c is much smaller (and the standard deviation in the final outer radius, SD = 211 μm, much larger) than for the initially uniform cell layer with the same flow rate and outlet pressure. Whilst the change in the depth of the cell layer between z = 7.36 cm and z = 7.38 cm in Fig. 6(d) appears sharp, it should be noted that the cell layer is less than 0.5 mm deep, so the gradient of the surface in this region is of order 1. Nevertheless, the change in depth may be sharper than would occur in practice, and reflects the simplicity of the form of shear stress dependence assumed in the cell layer growth law (35).
Reducing P l out , brings the shear stress on the cells near the outlet below Σ d , so that the cell layer grows outwards over its whole length to a mean outer radius of R = 870 μm c after 60 days. However, the cell layer remains thicker at its upstream end and the standard deviation in the final outer radius (SD = 16 μm) is larger than for the uniform cell layer. Similar results are observed when the initial depth of the cell layer increases with z (R z R z z ( , 0) = ( ) = (419 + 6. 1 × 10 ) μm , so that R = 705 μm c and SD = 14 μm, but grows over its whole length at lower outlet pressures.

Comparison with flow, solute transport and growth with ECS port closed
Finally we briefly describe the key differences in the flow, solute . Arrows indicate direction of increasing time, solid lines the dependent variables at regular intervals of 6 days, dash-dot line in (c) the shear stress cell death threshold Σ d , and dashed line in (d) the initial cell layer outer radius R z ( ) = 450 μm c init , transport and cell layer growth when the ECS port is closed.

Flow
When the ECS port is closed, there is no normal flow out of, or axial flow along, the ECS outer boundary, r R u e u e · = 0, · = 0 on = , e r e z e so the axial flow in the lumen dominates the radial flow and flow is negligible in the rest of the bioreactor (see Appendix D.1). For the cell types and range of flow rates considered here, the interstitial shear stress remains far below the cell death threshold, Σ d . Thus changes in the interstitial shear stress due to changes in Q l in , and P l out , do not affect the growth of the cell layer.

Solute transport
As before, we nondimensionalise the governing equations ((2)-(13), (14)-(36), (37)) and exploit the small aspect ratio of the bioreactor to simplify the resulting system as in Appendix A. The mass transport problem for the oxygen and lactate reduces to solving an advection-diffusion equation for the lumen concentration at each growth time step (see Appendix D.2). Since radial flow through the bioreactor is much weaker than when the ECS port is open, the rates at which oxygen is transported through the membrane to the cells and lactate is advected away from the cell layer are much lower. Consequently the minimum oxygen and maximum lactate concentrations are much lower and higher respectively. For example, for the uniform cell layer in Fig. 4(d) (R = 450 μm c ) the minimum oxygen concentration with the ECS port closed is 0.108 mol m −3 (compared to 0.216 mol m −3 with the port open) and the maximum lactate concentration is 0.119 mol m −3 (more than 60 times higher than with the port open).
With the ECS port closed, the leading order oxygen and lactate concentrations at the outer surface of the cell layer, and therefore the cell layer growth rate, depend on Q l in , (through its effect on the oxygen and lactate transport through the lumen) but not on P l out , . As when the ECS port is open, the growth rate of the cell layer decreases as it grows out into the ECS, since the surface oxygen concentration decreases with distance from the fibre. Since cell layer growth does not affect the dominant flow in the bioreactor, it only affects the oxygen and lactate distributions by altering uptake/production and diffusion within the . Lines and arrows as in Fig. 4. Chapman et al. Journal of Theoretical Biology 418 (2017) 36-56 ECS (it changes the latter as the solutes diffuse more slowly in the cell layer than in the free fluid in the ECS).

Growth
The reduced system for the flow, solute transport and growth with the ECS port closed is solved as described in Appendix D.3, and used to simulate the growth of initially uniform layers of rat cardiomyocytes and bovine chondrocytes as in Section 4.2.1. Fig. 7 shows how R c , the mean outer radius of the cell layer at t=60 days, varies with Q l in , for the cardiomyocytes and chondrocytes (R z ( ) = 450 μm c init , as in Fig. 3). The cardiomyocyte layer quickly reaches a steady state depth (R = 463-466 μm c ), since the surface oxygen concentration c| r R = c rapidly approaches C min . By contrast, the depth of the chondrocyte layer has not reached equilibrium after 60 days, when R = 656-657 μm c . In this case the oxygen uptake rate is lower, so c| r R = c remains well above C min as the cell layer grows. For both cardiomyocytes and chondrocytes lactate levels remain below the threshold for cell death, ), whereas the reduction in growth is very small for the chondrocyte layer (R = 656-657 μm c compared to 665-675 μm). The reduction in growth when the ECS port is closed is due to reduced radial flow and oxygen delivery to the cells. The effect is more pronounced for the cardiomyocytes as they consume oxygen at a higher rate, causing the oxygen concentration at the outer surface of the cell layer to decrease more rapidly as it grows. Fig. 7 shows that increasing Q l in , has only a modest effect on the final cell layer depth for both cardiomyocytes and chondrocytes. This is because increasing Q l in , increases only advective transport through the lumen at leading order when the ECS port is closed, leading to a marginal increase in the oxygen concentration at the inner surface of the membrane, which propagates through the membrane and cell layer by diffusion. The variation in the cell layer depth at t=60 days for the cardiomyocytes is much smaller when the ECS port is closed (SD = 0.5-1.4 μm) than when it is open (SD = 0.2-44 μm), as expected given that the cardiomyocyte layer grows very slowly. The reduction is ). Lines and arrows as in Fig. 4. Chapman et al. Journal of Theoretical Biology 418 (2017) 36-56 less pronounced for the chondrocytes (SD = 0.1-0.4 μm compared to SD = 0-2 μm). From these results it is clear that opening the ECS port and using higher flow rates and outlet pressures can enhance cell layer growth. However, it is important to know the level of shear stress that the cells can withstand.

Sensitivity of growth to lactate concentration
So far, we have assumed that the lactate concentration at which cell death occurs is the same for all cells (L = 0.4 mol m max −3 ), and have found that for the cell types and operating conditions considered the lactate concentration at the outer surface of the cell layer remains much lower than L max . In the absence of cell-type-specific data for L max , the sensitivity of the cell layer growth to L max warrants further investigation. Therefore we now simulate the growth of an initially uniform rat cardiomyocyte layer as L max varies in the range 0.001-0.4 mol m −3 . Fig. 8 shows how the growth of the cell layer (measured as the difference, R Δ c , between its mean final outer radius, R c , and initial outer radius, R = 450 μm c init , ) varies with L max when the ECS port is closed and when it is open. With the ECS port closed, the cell layer recedes if L < 0.14 mol m max −3 (i.e. if L max is to the left of the vertical line in Fig. 8(a)), the recession increasing as L max decreases, but grows if L > 0.14 mol m

Discussion
The simulation results show that opening the ECS ports, increasing Q l in , and increasing P l out , all improve growth, provided that the shear stress on the cells at the surface does not exceed that at which the cells die or detach from the layer (both processes being captured by recession of the cell layer in the model). Our results also show that the improvement in growth associated with opening the ECS ports is greater for cells with higher oxygen demands (Fig. 3).
Increasing Q l in , and/or P l out , increases the oxygen concentration and Fig. 7. Variation in the mean, R c, and standard deviation, SD, of the cell layer outer radius at t=60 days with Q l in , for (a) rat cardiomyocytes and (b) bovine chondrocytes, when the ECS port is closed. Parameter values: (a) as in Fig. 3(a); (b) as in Fig. 3(b). Solid line shows R c, dashed line SD. . Cell layer growth measured by R R R Δ = − decreases the lactate concentration at the cell layer surface, but also increases the shear stress on the cells. For shear-sensitive cells it is necessary to avoid imposing too high an outlet pressure and inlet flow rate. In view of this, we have used the model to predict the inlet flow rate and outlet pressure that maximise the growth of initially uniform layers of rat cardiomyocytes (assumed to be shear-sensitive) and bovine chondrocytes (shear-tolerant) when the ECS ports are open, and have predicted a lower optimal P l out , for rat cardiomyocytes (Fig. 3). This highlights the extent to which the optimal flow rate and outlet pressure will be cell-type-specific. Cells should be cultured at the highest possible Q l in , and P l out , for which the interstitial shear stress in the cell layer remains below Σ d , to maximise oxygen delivery and lactate clearance.
Although our model simulations predict that cells such as cardiomyocytes, with high oxygen demands (high V max O and C min ), should be grown with the ECS ports open, it may be better to keep the ECS ports closed for cells that are very sensitive to shear stress (have a very low Σ d ) or have sufficiently low oxygen requirements that the culture medium can be recycled directly from the outlet to the inlet (i.e. does not need to be replaced) without any adverse effect on the growth. When the ECS port is open, the predicted percentage increase in the cell number for an initially uniform layer of cardiomyocytes is 1300% after 60 days (for Q = 3.33 × 10 m s ), based on the volume change of the cell layer and assuming that the cells are all spherical with a fixed radius of 10 μm (Table B1) ). When the ECS port is closed, the percentage increases in cell number for cardiomyocytes and chondrocytes are 30% and 500% respectively (for Q = 3.33 × 10 m s l in , −8 3 −1 ). We anticipate that opening the ECS port and using a high flow rate and outlet pressure could lead to similar improvements in growth for rat hepatocytes and pancreatic beta cells to those for rat cardiomyocytes, given their similar oxygen requirements and higher shear stress tolerances (Tables 2 and 3).
For the cell types we have considered, the surface lactate concentration remains well below the estimated toxicity threshold, L max , for typical operating conditions (see , and oxygen is the growthrate limiting solute. However, in the absence of experimental estimates of L max , we have used the model to predict the sensitivity of the cell layer growth to L max . Our results suggest that the range of values of L max for which the lactate concentration is growth-rate limiting (or causes recession of the cell layer) is smaller when the ECS ports are open. Therefore it may be better for the ECS ports to be open if L max is not known for the cell type being cultured (provided the cells are not highly shear-sensitive).
Our simulations demonstrate that cell layer growth (and the optimal flow rate and outlet pressure) may be sensitive to variation in the initial cell layer depth along the fibre (Fig. 6). They predict that the cell layer will grow less and have a less uniform depth after 60 days if its initial depth is non-uniform. This suggests that cell seeding and initial culture conditions should be aimed at producing a uniform cell layer depth (by making the initial cell density on the fibre surface and the solute distributions as uniform as possible) to maximise growth and the uniformity of the final cell layer. However, in practice it can be difficult to achieve such uniformity in the initial cell density and therefore in the growth (Chapman et al., 2014). We have considered initial cell layers, , whose depth is uniform or decreases/ increases linearly with distance along the fibre as an initial exploration of the impact of non-uniformity on the cell layer growth. However, further investigation of choices for R z ( ) c init , is needed to confirm whether non-uniformity has only a negative impact on the growth.
across a range of experimental values of Q l in , and P l out , would also enable the prediction of optimal combinations of R z ( ) c init , , Q l in , and P l out , . There are several ways in which this study could be extended and the accuracy of its predictions improved. The first step would be to validate the model predictions against experiments for specific cell types. This would involve estimating parameters such as the lactate toxicity threshold L max and shear stress death threshold Σ d , and the death rates under excess lactate and shear stress, B B L d and B B s d . It may be necessary to use a more detailed growth law to capture the precise dependence of the proliferation and death rates on the oxygen and lactate concentrations.
Although contact inhibition is believed to be significant for the cell types considered here, it would be instructive to consider simulations in which growth occurs throughout the cell layer, rather than just at the free surface, to see how the growth predictions differ. An assessment of which assumption is more realistic, i.e. surface or volumetric growth, could then be made by comparing model predictions to experimental measurements of cell layer growth.
A simplifying assumption made in deriving the model was that the cell layer is initially axisymmetric and remains so as it grows. In practice, there will be some azimuthal variation in the cell layer depth, and this may affect cell layer growth. However, to include such variation in the model would require fully 3D simulations, which would be highly computationally expensive, and hence lie beyond the scope of this study.
Further possible extensions include using a multiphase description of the cell layer to account for variation in the cell layer porosity with time due to cell proliferation and extra-cellular matrix deposition (Pearson et al., 2015b), and using a more detailed model of cell metabolism incorporating glucose and protein transport (Shipley et al., 2009) and coupling between the oxygen, glucose and lactate concentrations (Casciari et al., 1992).

Conclusions
We have developed a model of cell layer growth in a single-fibre HFB by extending a previous HFB fluid and mass transport model (Shipley and Waters, 2012). Our model incorporates dependence of cell layer growth on local oxygen and lactate concentrations and on local fluid shear stress. This enables predictions to be made of the effect of varying the lumen inlet flow rate and lumen outlet pressure on the growth of different cell types. We have observed that opening the ECS ports and increasing the flow rate and outlet pressure enhance growth for shear-tolerant cells by increasing oxygen delivery, but that it may be necessary to use lower outlet pressures for shear-sensitive cells to maximise growth.
Once parameterised against experimental data, our model could be used to make quantitative predictions of the optimal operating conditions for different cell types.  The advection-diffusion equations in the lumen, membrane and ECS, (16) and (17), and the reaction-advection-diffusion equation in the cell layer (18)   are the ratios of the cell layer recession rates due to excess lactate and excess shear stress, and the two combined, to the cell proliferation rate.
The initial condition for the position of the outer surface of the cell layer becomes  The permeability of the cell layer, k c , is very difficult to measure experimentally and varies significantly with the initial seeding density and cell type. Assuming that the cell layer is densely packed with cells, we can estimate its permeability from its porosity ϕ c and the average cell radius d cell using the Kozeny-Carman equation (Kozeny, 1927;Carman, 1937)    reasonably well with the range of values 10 −12 -10 −8 used in other studies (Heath et al., 1990;Hay et al., 2001).

B.2. Cell layer growth rate
We estimate the growth rate of the cell layer per unit concentration, A p , for different cell types from their cell doubling times (see Table B1). As we have assumed that only the outermost layer of cells in the cell region proliferates, we can assume that in the absence of cell death the depth of the cell region will increase by approximately one cell diameter d cell in the cell doubling time T d , so that Estimates for A p for different cell types are given in Table B1. For simplicity we assume that the baseline recession rate of the cell layer due to cell death and detachment from excess lactate or shear stress, B d , is the same as the rate at which the cell layer grows due to cell proliferation, A C p in , and the dimensionless weights for the excess lactate and shear stress recession rates, B L and B s , are both 1, so that β For the cell types in Table B1 the cell layer growth timescale R A C /( ) l p i n is ∼10 s 5 , whereas the advection timescale L U / is in the range 0.75-7.5 s for the flow rates in Table 1, and the timescales for radial diffusion and uptake/production of oxygen and lactate are in the ranges 13-290 s and 30-3000 s respectively. Hence, we are justified in assuming that the fluid and mass transport equations (2) and (16)-(18) are quasi-steady on the growth timescale.

B.3. Shear stress parameters
To estimate the dimensionless shear stress constant, κ k R d = 8 /( ) s c l , we need to estimate the size d of the interstitial spaces in the cell layer. Using the empirical relationship for the permeability of a membrane and the average pore diameter derived by Bear (1988) k Cd = , where C = 6.54 × 10 −4 , and our estimate for k c from Appendix B.1 (7.5 × 10 m −13 2 ) gives d k C ≈ / =3.4×10 m c −5 , so κ ≈ 8.9 × 10 s −4 . The maximum shear stress that the cells can withstand without dying or detaching from the surface of the cell layer, Σ d , depends on the cell type. While bovine chondrocytes grow well under shear stresses of up to 2 Pa (Smith et al., 1995), rat cardimyocytes can be damaged or killed by shear stresses higher than 0.16 Pa (Radisic et al., 2005). The level of shear stress that cells can withstand generally decreases (from order 1 Pa to order 0.1 Pa and smaller) as the length of time that they are exposed to it increases (Chapman et al., 2014, Supporting Information, Table S6). Studies with human foreskin fibroblasts (HFFs) and mouse osteoblasts in micro-channel bioreactors have shown that cell detachment can occur at shear stresses as low as 0.03 Pa (Korin et al., 2007;Leclerc et al., 2006). To incorporate the possibility of thinning of the cell layer due to cell detachment, we use values of Σ d , appropriate to the cell types we consider, in the range 0.03-2 Pa. . All other parameter values as in Table 1. Lumen axis running from right to left. Black dashed lines indicate the membrane surfaces and the white dashed lines the outer surface of the cell layer.   Fig. E1, and oxygen transport parameter values for rat cardiomyocytes as in Table 2. Dashed lines as in Fig. E1.