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

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

Introduction

The need for a reliable and sufficient source of replacement tissue and organs is constantly increasing, due to our ageing population and consistent lack of donors. The potential for tissue engineering to meet this demand is great, but there are still many barriers to be overcome before the field can provide alternative treatments on a clinical scale. For each tissue type, a different experimental technique must be developed to ensure that the engineered substitute is viable and can perform the same essential functions as the original tissue. This often means being able to mimic the in vivo environment of a particular tissue and has resulted in a vast array of protocols, bioreactors, culture conditions, and scaffold materials (see e.g. Martin et al. 2004; Martin and Vermette 2005; Pörtner et al. 2005; Stock and Vacanti 2001). For each set-up and cell type, optimal operating conditions must therefore be determined, a process which can be extremely time-consuming and expensive to resolve purely experimentally.

Mathematical modelling of bioreactor systems can be of great benefit for several reasons. It is reproducible and efficient, allowing the large numbers of parameters to be investigated relatively quickly and cheaply. Moreover, it can give insight into the combined effects of the physical processes involved in a particular set-up, or even focus on one process in isolation—this is hard to achieve in an experimental set-up. An example of one such process, and the focus of this work, is mechanotransduction: the mechanism by which forces are converted into biochemical signals and integrated into a cellular response. Below, we review a selection of mathematical models that specifically incorporate the effects of mechanotransduction. For more general detail on mathematical models in tissue engineering, we refer to the recent review by O’Dea et al. (2013a).

Previous work on multiphase modelling of mechanotransduction in tissue engineering includes a series of papers by O’Dea et al. (2008, 2010, 2013b). In O’Dea et al. (2008), the authors develop a two-phase model of tissue growth in a perfusion bioreactor. They include dependence of cell proliferation, extracellular matrix (ECM) deposition and cell death on the local cell density and fluid pressure. Results indicate that these effects can dramatically alter the composition of the resulting cell construct. This model was extended to include an additional porous scaffold phase in O’Dea et al. (2010), in which the system was simplified analytically by taking the long-wavelength limit appropriate for a small channel aspect ratio. The effect of fluid shear stress, as well as cell density and fluid pressure, was considered. Findings supported the earlier conclusion that mechanotransduction can significantly affect cell distribution. The same set-up was modelled and solved numerically by Osborne et al. (2010), highlighting situations in which a numerical approach is necessary, and when the analytical limits taken in O’Dea et al. (2010) are relevant. A further extension of this work is given in O’Dea et al. (2013b) in which the degradation of a scaffold is included, and whose spatially non-constant porosity is informed by experimental data.

As well as cell proliferation and death, mechanotransduction can affect cell differentiation. This is the focus of a paper by Byrne et al. (2007), in which finite element analysis is used to model a poroelastic tissue–scaffold construct for bone regeneration, with the objective of determining optimal scaffold porosity under different mechanical loading conditions. The dependence of stem cell differentiation on fluid velocity and shear strain is considered by employing the mechanoregulation algorithm developed by Prendergast et al. (1997). Depending on the stimulus level, differentiation into fibroblasts, chondrocytes or osteoblasts can occur. Aiming to identify possible treatment options in fracture healing, Lacroix and Prendergast (2002) develop a three-dimensional model of a human tibia fracture, again using the mechanoregulation algorithm proposed in Prendergast et al. (1997). Two different compressive loading magnitudes were simulated, and it was shown that only the lower load led to successful healing.

A finite element approach was also taken by Driessen et al. (2003), who modelled the effect of the different mechanical loading regimes resulting from closed and open configurations on collagen fibre content and orientation in an aortic valve. The mechanical properties of the tissue construct were shown to depend on the type of loading, as the fibres aligned with the principal strain directions, and the results were validated against experimental data where possible. Cardiac tissue was also the focus of the work in Latimer et al. (2003), who determined the stress and strain distribution in two- and three-dimensional tissues. Analytical solutions were found and used to demonstrate the strain in and around an ischemic region, and have the potential to be further used to predict mechanotransduction effects on the tissue leading to arrhythmias, or following hypoxia. Finally, the importance of scaffold design in bone tissue engineering motivated the study by Sanz-Herrera et al. (2009). The need for an implantable construct that has the required macroscale properties, but which is affected by biophysical phenomena on the microscale, led to a multiscale model which incorporated a realistic scaffold microstructure. Bone regeneration was assumed to be dependent on mechanical stimuli, and properties for the macroscale model such as the scaffold stress and strain and cell diffusion were obtained through the analysis of the underlying microstructure, which in turn evolved as a result of macroscale cell migration.

In this paper, we consider various experimentally relevant case studies, motivated by specific cell types, to investigate the effect of fluid shear stress on the proliferation rate, and thus distribution, of cells in a hollow fibre membrane bioreactor (HFMB) (see Fig. 1). A HFMB consists of a cylindrical, glass module with a port at each end of the extracapillary space (ECS) and a porous hollow fibre inserted through the centre. A number of different cell seeding and flow regimes can be employed; for this work, we consider the set-up in which the cells are seeded in a scaffold in the ECS between the hollow fibre and the outer glass wall. Culture medium is pumped in via the lumen inlet and upstream ECS port. The pressure is set at the downstream lumen end via a clamp, and the downstream ECS port is left open to the atmosphere. Depending on the flow and pressure conditions, the fluid may pass through the hollow fibre walls and then flow out of either the lumen outlet or downstream ECS port.

Fig. 1

Top Photograph of a single HFMB module (ruler scale in cm) as shown in Pearson et al. (2013). Bottom cross section of the boxed region (not to scale), where the lower half (not shown) has been excluded based on symmetry. This depicts the idealised two-dimensional modelling region with the axis running along the lumen centreline. The solid black arrows show the direction and location of the fluid fluxes into the system, and the star denotes the origin

We build on the model developed in Pearson et al. (2013) which describes fluid flow, solute concentration and cell distribution in a HFMB. In Pearson et al. (2013), the entire ECS was filled by the cell–scaffold construct, through which media could flow. Fluid velocities in the membrane and ECS were assumed to be an order of smaller than those in the lumen as a result of the resistance to flow through these porous regions, where is the lumen aspect ratio. Here, we consider the alternative experimental set-up in which the cell–scaffold construct only fills part of the ECS (that directly above the porous membrane, defined as the cell layer), and media can flow freely throughout the remainder. The width of this cell layer can be controlled experimentally and will be varied as part of our investigation. Fluid is pumped into the system via both the lumen inlet and upstream ECS port; as a result, we expect the fluid velocity in this upper fluid layer to be of the same order of magnitude as the lumen. As well as enabling the relevant solutes to be delivered to, or removed from, the cells more quickly, we expect the cells to be exposed to a higher and potentially significant level of fluid shear stress as the flow through the cell layer will be enhanced by the upper fluid flow above. The impact of this set-up remains an open question experimentally: depending on the cell type, fluid shear stress can either be beneficial or detrimental to cell growth and survival. In this work, we therefore consider mechanotransduction effects, which were not included in Pearson et al. (2013). Specifically, we explore the effect of fluid shear stress on cell yield and distribution, via changes in the fluid flux into the upper fluid layer. In addition, the fluid flux affects the concentration of a nutrient (taken to be oxygen) which can limit the cell proliferation rate. We determine the range of possible behaviours from this set-up and find optimal flow rates for which the fluid shear stress has an advantageous effect: that is, it enables a more spatially uniform cell population and/or a higher cell yield to be obtained.

Paper outline

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

Model description

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

Governing equations

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

Boundary conditions

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

Model reduction

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

Parameter values

Typical dimensional parameter values are taken from the literature or guided by our experimental collaborators where possible. We present these in Table 1 along with corresponding dimensionless parameters in Table 2. Dimensional parameter values for which data could not be found are omitted from Table 1; instead, their dimensionless equivalents are defined in Table 2 and values chosen so that our asymptotic analysis retains as many features as possible at leading order. A few choices are worth mentioning here: for a more in-depth discussion see Pearson et al. (2013). Firstly, we note that the aspect ratio of each region is small; in particular, we define the lumen aspect ratio by and use this as the relevant parameter for our asymptotic expansion in Sect. 3.2. In addition, the reduced Reynolds number in the lumen and upper fluid layer is which justifies neglecting inertial effects in these regions in Sect. 2.1 and motivates our choice of non-dimensionalisation. We set the viscosity ratio for some constant . This is motivated by the fact that we would expect the effective (macroscale) viscosity of the cell phase to be much greater than that of water as a result of the cytoskeletal network, the cell–scaffold affinity and microscale cell–cell interactions. We choose the dimensionless cell–scaffold drag and water–scaffold drag to be of , so that their effects are retained in the leading order model, but assume that (as would be expected physically).

Table 1

Dimensional parameters

Parameter	Dimensional value and units	Definition
h_1	200μm^a	Lumen height
h_2	200μm^a	Porous membrane height
h_3+h_4	600μm^a	ECS height (cell layer + upper fluid layer height)
L	0.1m^a	Length of bioreactor modelling domain
ρ_w	1 g cm^-3a	Water density
μ_w	10^-3 Pa s^a	Water viscosity
k	6.73×10^-16 m^2b	Porous membrane permeability
p_atm	14.69 psia^a	Atmospheric pressure
Q_l,in,Q_f,in	1.02×10^-11-1.02×10^-8 m^2 s^-1c	Lumen/upper fluid layer inlet flux
Γ_nw	5.79×10^-6 s^-1c	cell proliferation rate coefficient
Γ_wn	4.13×10^-7 s^-1	Cell death rate coefficient
c_in	0.22 mol m^-3a	Inlet solute concentration
D	3×10^-9 m^2 s^-1a	Solute diffusion coefficient in water
U^∗	1.0239×10^-8 m s^-1c	Typical porous membrane velocity
C^∗	0.22 mol m^-3a	Typical solute concentration

Values taken from Shipley et al. (2010)

Experimentally obtained values

Values based on estimations by our experimental collaborators

Table 2

Dimensionless parameters

Parameter	Definition	Value	Restriction (if any)
h^_2	h_2/(εL)	1	h^_2>0
h^_3+h^_4	(h_3+h_4)/(εL)	3	h^_3,h^_4>0
θ_s	Scaffold volume fraction	0.4^a	0<θ_s<1
ε	h_1/L	2×10^-3	0<ε≪1
λ	μ_w/(εμ_n)	1	ε≪λ≪1/ε
ε^2Pe	εLU^∗/(λD)	6.83×10^-4	ε^3≪ε^2Pe≪ε
ε^2Re	ερ_wLU^∗/(λμ_w)	2.05×10^-6	ε^2Re≪1
ϕ	Porous membrane porosity	0.77^b	0<ϕ<1
κ	k/(λε^5L^2)	2.1	ε≪κ≪1/ε
Q^_i,in,i=l,f	λQ_i,in/(LU^∗)	0.01–10	ε≪Q_i,in^≪1/ε
ν^	λε^3Lν/(μ_wU^∗)	0.3^c	–
δ^	λε^3Lδ/(μ_wU^∗)	0.1^c	–
η^	λε^3Lη/(μ_wU^∗)	0.3^c	–
ζ_ns	γ_nsL^2λε^3/μ_w	1	ε≪ζ_ns≪1/ε
ζ_ws	γ_wsL^2λε^3/μ_w	0.1	ε≪ζ_ws≪ζ_ns
Γ^_nw	LΓ_nw/U^∗	56.52	ε≪Γ^_nw≪1/ε
Γ^_wn	LΓ_wn/U^∗	4.04	ε≪Γ^_wn≪1/ε
Γ^_R1	L^2Γ_R1/(C^∗D)	50	ε≪Γ^_R1≪1/ε
K^,K^_1	K/C^∗,K_1/C^∗	1	–
c^_in	c_in/C^∗	1	–
P_d	–	2.5	–

Values taken from Lemon et al. (2006), Meneghello et al. (2009), and O’Dea et al. (2010)

We define the velocity scale to be a typical horizontal velocity of the water in the porous membrane. We then take the corresponding horizontal velocity scale in the cell layer to be , and in both the lumen and upper fluid layer to be , so that the same pressure scales apply throughout the system in our non-dimensionalisation. These choices are also motivated by the fact that we would expect the porous structures to hinder the fluid flow in the membrane and cell layer. In addition, this choice enables substantial analytical progress to be made in the mathematical reduction in Sect. 3.2. We note that other limits could be considered, for instance taking the cell and water phases to be of comparable viscosity (as in O’Dea et al. 2010, 2013a) or the cell phase to be inviscid (as in Lemon et al. 2006).

We take the reduced Péclet number in the lumen and upper fluid layer, Pe, to be order so that diffusion dominates advection throughout. The concentration scale and inlet concentration are both set to be a typical oxygen concentration used in cell culture (Shipley and Waters 2012), as oxygen is the solute of interest in Sect. 4. The two-dimensional lumen and upper fluid inlet fluxes and (which come into the up- and down-stream boundary conditions in Sect. 3.3) have been converted from three-dimensional experimental values by first calculating the corresponding velocity and then multiplying by the length scale in the direction, . The downstream lumen outlet pressure and atmospheric pressure are also introduced in Sect. 3.3. The cell proliferation/death rates and solute uptake rate are parameters in the constitutive laws for and , both of which we define in Sect. 4. Based on estimations by our experimental collaborators, we take the (dimensional) cell proliferation rate coefficient 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 .

Derivation of the reduced model

Having non-dimensionalised the governing equations and boundary conditions (see “Appendix”), we asymptotically expand each dependent variable in powers of , for instance setting and similarly for the remaining velocities, reduced pressures, solute concentrations and volume fractions. In the following, we omit the subscript from leading order variables except where needed for clarity. Equating coefficients of in the lumen then yields in the membrane in the cell layer and in the upper fluid layer The leading order boundary conditions on the lumen centreline are on the lumen/membrane interface on the membrane/cell layer interface on the cell layer/upper fluid interface and finally on the bioreactor top

Consideration of the above solute equations and boundary conditions shows that at leading order, there is a global concentration that is independent of : . To close the problem and find , we must also consider the solute equations at . We note that the components of the solute equations and boundary conditions in each layer indicate that , and so (in the lumen, membrane, cell layer and upper fluid layer, respectively)

together with boundary conditions A substantial amount of analytical progress can be made in the distinguished limit considered above. Equations (3.1c), (3.3g) and (3.4c) respectively tell us that the leading order fluid pressures in the lumen, cell layer and upper fluid layer are independent of . From (3.3e), and given is constant, we can thus conclude that the cell volume fraction is also a function of and only, which allows substantial simplifications to be made. We can obtain the following expressions for the leading order variables in terms of and : where is given by Hence, at leading order, we have Poiseuille flow in the lumen and upper fluid layer and transverse flow only in the porous membrane. In the cell layer, the horizontal fluid flow is Darcy-like. The variables and can be determined from the leading order conservation of mass Eqs. (3.3b, c), respectively.

In addition, we have the following equations for , the global leading order concentration and : where the cell flux is defined byand the solute flux byThe equations for and come from integrating the respective conservation of mass Eqs. (3.1a), (3.3b) and (3.4a) across the lumen, cell layer and upper fluid layer. The remaining equation (3.23c) has been obtained by adding together the solute equations (3.10a–d) after averaging each across the appropriate region. We now impose up- and down-stream boundary conditions to close the problem.

Boundary conditions for the reduced model

We mimic the experimental set-up where fluid is pumped into the bioreactor via both the lumen inlet and (upstream) ECS port, and so impose where and are prescribed (dimensionless) volume fluxes per unit length in the direction (perpendicular to and ) into the lumen and upper fluid layer, respectively, and are assumed constant. We note that in this case, the dimensional is of the same order as the dimensional given that the velocity scales in the lumen and upper fluid layer are the same. We also prescribe a down-stream lumen pressure and mimic the atmospheric pressure conditions at the down-stream ECS ports by settingFor the cells, we impose no flux out of the modelling domain (which would be achieved using filters in an experiment):Finally, we prescribe an inlet solute concentration and assume no diffusive flux at the downstream end of the bioreactor: As there is no downstream constraint on the solute concentration experimentally, condition (3.31) is motivated by the fact that we would expect the concentration to be constant in space as it leaves the bioreactor, due to the effects of diffusion.

Applying boundary conditions (3.26)–(3.28) allows us to solve (3.23a, d) explicitly for and , givingand thereby specific expressions for and from (3.16)–(3.19). In particular, we note the form for the fluid pressure in the cell layer, which is equal to (see 3.18b),

Summary of reduced model

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

Numerical results

The reduced system (3.34)–(3.40) was solved numerically using the method of lines, first discretising in and then performing the time integration using the MATLAB function ode15s. We are interested in steady states of the system, as these arise on timescales comparable to the long culture time of experiments (see Pearson et al. 2013 for a more thorough discussion). These states were found to be independent of the choice of initial condition in our simulations, but for completeness, we note that in the simulations presented in this paper, we set for .

We consider the effect of the fluid shear stress on the cell distribution and yield, in the case where the solute of interest is a nutrient (e.g. oxygen), and with the aim of determining optimal experimental conditions which result in both a spatially uniform cell population and high cell yield.

In each of the following sections, we consider the Michaelis–Menten type reaction term used in Pearson et al. (2013) for nutrient-limited proliferation:We take the cell mass transfer term to have the formwhere is the (dimensionless) fluid shear stress in the cell layer. If were equal to unity, then we would recover the form used in the case of the nutrient-driven proliferation discussed in Pearson et al. (2013). In the following analysis, we take to have different forms depending on the cells’ response to shear. As we effectively have Darcy flow of water through the cell layer, we use the following estimate for motivated by a similar expression in Whittaker et al. (2009):based upon the assumption that there is Poiseuille flow through each of the ‘pores’ in the cell layer, where the pores are approximated as circular ducts. Given the assumed form for , and the expression for in (3.33), we can see that the fluid shear stress is dependent on both the upper fluid layer flux and cell layer width (via ), both of which can be controlled experimentally. In what follows we investigate the sensitivity of our results to variations in these two parameters. Moreover, for a given flow rate and cell layer width, we note that variation in arises purely through variations in .

In each of the following case studies, we present results for the steady-state cell volume fraction and use the mean and standard deviation of the cell volume fraction to characterise the cell yield and distribution, where

Shear stress-enhanced proliferation

First of all, we consider the case when cell proliferation is enhanced by increased levels of fluid shear stress. This has been observed in cells such as osteoblasts (Kapur et al. 2003), chondrocytes (Malaviya and Nerem 2002) and (more controversially) endothelial progenitor cells (Yamamoto et al. 2003), and is thought to be an important factor in ensuring the mechanical integrity of tissue such as bone and cartilage (Bancroft et al. 2002; Gemmiti and Guldberg 2006). The exact nature of this dependence is not known and will depend on the cell type in question. To capture the qualitative aspects of this dependence, we setso that the cell proliferation rate increases with shear stress up to a maximal level, but (as would be expected) high levels of shear stress result in a reduced proliferation rate (see Fig. 2). The coefficients in (4.6) have been chosen to capture the possible behaviours of the cell population within the range of computed shear stresses in our model.

Fig. 2

Plot of , showing the dependence of cell proliferation on the fluid shear stress in the case where increasing shear stress first enhances cell proliferation (solid line, Sect. 4.1), and where increasing shear stress inhibits cell proliferation (dashed line, Sect. 4.2)

Results from this case study can be seen in Fig. 3, where the steady-state cell population distribution obtained with as in Eq. (4.6) is compared to that when (denoted the shear-insensitive population, i.e. when the fluid shear has no effect on the cells) for a range of upper fluid layer flow rates . From Fig. 3a, b, we see that for the lowest flow rate of , the cell volume fraction of the shear-sensitive (hereby called shear-enhanced) population is less than that of the shear-insensitive population for , but then increases above the shear-insensitive volume fraction once is increased to 0.5. This corresponds to the shear stress exceeding the value at which is increased above 1 (around 2.5) and is to be expected given our chosen form of . The increase in is observed across the whole bioreactor domain, instead of at a specific point, as a result of the relatively small variation in (and hence ) for , which in turn means that there is little spatial variation in . At the intermediate flow rate (), we also see that the shear-enhanced population becomes much more spatially uniform than the shear-insensitive cells. For the highest flow rate , the shear-enhanced decreases at each point in as the shear stress reaches levels which decrease the cell proliferation rate. In addition, greater up- and down-stream variation is seen in for both populations, but the changes are less extreme for the shear-enhanced cells. Hence, overall, the effect of the shear stress is to ‘smooth out’ the cell population: an increase in corresponds to a decrease in , and thus a higher value of . Thus, regions where is relatively high experience a fluid shear stress sufficient to decrease and hence the proliferation rate, whilst regions where is relatively low experience a lower shear stress and therefore increased proliferation rate.

Fig. 3

Plots of for a range of values of the upper fluid layer flux , in the a shear-insensitive and b shear-sensitive (enhanced proliferation) cases. The lumen flux , the cell layer width and other parameter values are as in Table 2

The effect of shear on the cell yield and population uniformity can be seen in Fig. 4, and the plots support the findings discussed above. The mean and standard deviation of are plotted for a range of values of , for the shear-enhanced and shear-insensitive cases. This clearly shows that the shear-enhanced population has a lower yield and greater standard deviation than the shear-insensitive population for low values, but that this relationship is reversed when increases beyond a critical value just below . This switch corresponds to the critical value of being reached at which , and the cell proliferation rate is enhanced. As can be seen in the graphs for both and , the change occurs relatively abruptly; this is due to the fact that (as mentioned above) there is relatively little variation in at this value of , and so the fluid shear stress is relatively constant along the domain. The standard deviation of the shear-enhanced population remains lower than the shear-insensitive population for all subsequent values of , but the yield drops below that of the shear-insensitive population once another critical value of is reached, around 0.7, when the value of falls below . Thus, between these critical values (given by vertical dashed lines in Fig. 4), there is an optimum range of within which shear stress helps to promote cell yield and also population uniformity. Within this range, the most uniform population is obtained at the (optimal) flow rate of . We also note that, for a flow rate above the first critical value, the shear-enhanced population is always more spatially uniform than the shear-insensitive population for a given yield.

Fig. 4

Plots of (a) the mean and (b) the standard deviation of , for a range of values in the shear-insensitive, shear-enhanced proliferation and shear-limited proliferation cases. The vertical dashed lines indicate the optimal flux range for the shear-enhanced case, within which the cell yield is higher than the shear-insensitive population. The vertical dash-dotted line indicates the critical flux value for the shear-limited case, below which the shear-limited and shear-insensitive cases are indistinguishable. The lumen flux , the cell layer width and other parameter values are as in Table 2

Shear stress-limited proliferation

In the second case study, we take to be of the following formso that cell proliferation is inhibited for sufficiently high values of fluid shear stress (see Fig. 2). Cells that have been shown to be sensitive to shear stress in this way include smooth muscle cells (Papadaki et al. 1996) and certain types of endothelial cells (e.g. human aortic endothelial cells Imberti et al. 2002 and human umbilical vein endothelial cells Akimoto et al. 2000). As in the previous section, the exact form of the cells’ shear stress dependence is unknown, and our choice of and corresponding parameter values have been chosen so that our model captures the range of possible behaviours. Throughout this section, the cell population with proliferation rate limited by shear stress will be known as shear-limited.

Figure 5 shows the results obtained from this case study. For low values of , the shear stress has no effect on cell proliferation (as is equal to ) and thus the shear-limited and shear-insensitive population distributions are the same. As is increased, however, the downstream region, where is greater, experiences a higher level of shear stress, and thus a decrease in proliferation rate of the shear-limited population. This results in the shear-limited population becoming more uniform than the shear-insensitive population. In contrast to the previous section, we note that regions of lower cell volume fraction are unchanged by the shear sensitivity, since in this case, proliferation is never enhanced by the fluid shear. Thus, the increased uniformity of the shear-limited population is at a cost of a lower overall volume fraction. The plots of and in Fig. 4 show that there is just one critical value of in this case (around 0.4, given by the vertical dash-dotted line in Fig. 4), below which the shear-limited and shear-insensitive cases are indistinguishable. Above this, the shear-limited population is more uniform than the shear-insensitive population, but also has a lower yield. However, Table 3 demonstrates that the percentage decrease in is more significant than the percentage decrease in for , and thus, this provides a potentially optimal flow rate range when a uniform cell population is more important than a high cell yield.

Fig. 5

Plots of in the shear-sensitive (limited proliferation) and shear-insensitive cases, for a range of values of the upper fluid layer flux . The lumen flux , the cell layer width and other parameter values are as in Table 2

Table 3

Comparison of cell yield and standard deviation for the shear-limited and shear-insensitive cases, for a range of values

Q_f,in	Stress–insensitive	Stress–limited	% Change
μ
0.4	0.4557	0.4557	-5×10^-4
0.6	0.4548	0.4493	-1.209
0.8	0.4523	0.422	-6.695
1	0.4464	0.3827	-14.28
σ
0.4	0.0174	0.0174	-0.013
0.6	0.0279	0.0234	-16.06
0.8	0.0438	0.0343	-21.67
1	0.0671	0.053	-21.0

Sensitivity to cell layer height

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

Fig. 6

Plot of versus the cell layer width , for fixed

Fig. 7

Plots of the mean and standard deviation of , for a range of values in the shear-insensitive, shear-enhanced proliferation and shear-limited proliferation cases, for (a, b) and (c, d) . The vertical dashed lines again indicate the optimal flux range for the shear-enhanced case, and the vertical dash-dotted line the critical flux value for the shear-limited case. The lumen flux and other parameter values are as in Table 2

Discussion

We have developed a multiphase model of fluid flow, solute transport and cell distribution in a simplified HFMB. This extends the work presented by Pearson et al. (2013), by considering the alternative set-up in which the cell–scaffold construct only partially occupies the ECS and there is an additional upper layer of free-flowing fluid. Fluid is pumped into the upper fluid layer via the upstream ECS port, and hence, this region experiences higher flow rates than when fluid is only supplied at the lumen inlet. It is therefore hypothesised that this could result in a greater fluid shear stress in the cell layer which could have an important, and possibly beneficial, effect on cell dynamics. This is confirmed from our model results which imply that the fluid shear stress in the cell layer is dependent on both the upper fluid inlet flux and cell layer width, both of which are varied as part of our investigation.

We presented results from two different case studies. In the first (relevant to osteoblasts, chondrocytes or endothelial progenitor cells), cell proliferation rate was enhanced by fluid shear stress, but only up to a certain critical level. Results demonstrated that the mechanosensitive cell population required a certain minimum ECS flow rate in order to reach yields comparable to the shear-insensitive population. For flow rates that resulted in fluid shear within the ‘enhancing’ range, however, the mechanosensitive population had a higher yield and was more uniform than the shear-insensitive cells. This uniformity persisted even once the shear levels became detrimental to the proliferation rate and the yield decreased. In the second case study, fluid shear had no effect on cell proliferation until a critical value was surpassed, after which it decreased the proliferation rate. This could be applicable to smooth muscle cells, or certain types of endothelial cells. Results again showed that the mechanosensitive cell population was more uniform than the shear-insensitive cells, but always had a lower cell yield. However, an ideal flow rate range could again be found, within which the mechanosensitive population is more uniform than the shear-insensitive population, and the yield has not yet been significantly compromised. Finally, we considered the effect of the cell layer width on these results and found that a thinner cell layer experiences a lower level of shear stress than a thicker layer, and hence can withstand higher flow rates.

The results presented here are qualitative in nature, demonstrating the range of expected behaviours in each situation. In particular, the functional forms for are constitutive and relatively simple, having been chosen to clearly demonstrate the possible effects of mechanotransduction on cell distribution in a HFMB. Any choice of would need experimental validation to be applied to a specific experimental setup. There are currently no relevant experimental results to which we can compare the results from our mathematical model, as the modelling is ahead of the experiments. However, the model framework developed here is extremely flexible and permits interrogation of different forms for any of the parameter values and constitutive laws presented. Hence, once relevant data become available, the model could be readily adapted and applied to a particular cell population, and proper validation would be possible. This could then be used in a predictive way to stimulate future experimental work.