Optimal metabolite transport in hollow fibre membrane bioreactors

Abstract

Hollow fibre membrane bioreactors provide a fast and efficient method for engineering functional tissue for use in medical treatments. Flow is utilised to overcome mass transport limitations by perfusing a nutrient-rich culture medium through the fibre lumen, which can then transport along the fibre lumen or across the porous membrane wall. Cells seeded at the outer membrane wall consume the nutrient and subsequently produce waste metabolites, which are transported away through an external extra-capillary space (ECS) along with excess nutrient. We present and investigate a two-dimensional axisymmetric model for fluid flow and solute transport through a single-fibre bioreactor configuration, with cells seeded to the external fibre wall. Fluid flow is modelled by steady lubrication and Darcy equations, which are coupled to the solute transport problem modelled by a system of advection–diffusion equations, supplemented with a reaction term to model the cell layer. Our model analysis reveals how spatially varying wall permeability distributions can be utilised to provide uniform nutrient delivery to a spatially uniform, homogeneous cell population. We also reveal how maximising the transmural pressure drop across the membrane wall is the dominant mechanism for waste removal rather than traditional experimental methods of flushing the ECS.

Impact Statement

Mathematical models are used to understand and improve nutrient transport to cells and waste product removal in a bioreactor, for applications in growing organs and tissue for transplantation. Typically, the nutrient concentration that cells experience depends on the location along the porous membrane wall on which they are seeded. This can result in cells growing at different rates, depending on how much nutrient they receive. Furthermore, under certain operating conditions, some cells may receive too little nutrient and die as a result. We show that by spatially varying the permeability of the membrane scaffold, advective nutrient delivery can be controlled to ensure that all cells receive the same amount of nutrient and grow at a constant rate, regardless of where they are located. Our model also provides mechanistic insight into how transport mechanisms change throughout the system, and guides experimental strategies to maximise nutrient delivery to, and removal of waste products from, the cells.

1. Introduction

Tissue engineering aims to circumvent the challenges associated with organ donation by synthetically repairing and replacing damaged organs and tissue (Vacanti & Vacanti 1991). One approach to tissue engineering is to use autologous cells taken from individual patients and grow tissue that avoids immune rejection and bio-compatibility issues. These issues have historically been significant bottlenecks in translating tissue engineering approaches into the clinic (Sarkar et al. 2023). To meet the clinical demand for the treatment of damaged or diseased organs, we must develop quick and efficient methods of engineering functional three-dimensional tissue for transplantation.

Cell growth is very sensitive to the local chemical and mechanical environment, which affects cell signalling pathways that promote cell proliferation and differentiation (Vance et al. 2005). Hollow fibre membrane bioreactors (HFMBs) offer the ability to carefully control the local cell environment via fluid flow, and allow careful tuning of factors such as the pH, metabolite (e.g. cell nutrients, growth factors and waste products) concentration and mechanical stress. Experimentally, cells are seeded on the outer surface of one or several porous hollow fibre membrane scaffolds, which are then cultured in a bioreactor, as shown in figure 1. Nutrient and growth factors are delivered to the system via flow into the fibre lumen(s), where they transport to the cells by advection and diffusion through the porous membrane wall. Flow of fresh culture media can also be introduced into the extra-capillary space to wash away unused nutrient, alongside waste metabolites produced. The porous hollow fibre scaffolds protect cells from excess shear stress from the fluid flow through the lumen, while also offering a large surface area–volume ratio for efficient nutrient delivery to the cell population (Wung et al. 2014). These HFMBs have shown promise across a wide range of applications, such as in the development of bio-artificial livers and in bypass surgeries (Sachlos & Czernuszka 2003; Seifu et al. 2013).

Figure 1. Schematic of HFMB adapted from Burova, Wall & Shipley (2019). The nutrient-rich medium enters the lumen inlets, which can transport along the lumen or across the porous membrane walls (grey) to cells seeded externally (pink). An external flow of culture medium (blue) in the extra-capillary space washes away excess nutrient alongside any waste metabolites produced by the cells.

Complementary to experimental work, mathematical models of HFMBs are used to improve experimental design by efficiently identifying parameter regimes that maximise tissue growth and bioreactor efficacy (O'Dea, Byrne & Waters 2012). Such models typically use steady Stokes and Darcy equations in a two-dimensional axisymmetric geometry to describe the fluid flow, coupled to advection–diffusion–reaction equations to describe metabolite transport and uptake on time scales where time-dependent cell layer growth can be neglected (Piret & Cooney 1991; Ye et al. 2006; Shipley et al. 2011). In order to generate a mechanistic understanding of the underlying system, significant simplifications or model reductions are often employed that limit the applicability of the models. For example, radial advection across the membrane and extra-capillary space (ECS) is often neglected in the analysis of metabolite transport (Piret & Cooney 1991; Ye et al. 2006; Pearson et al. 2014), despite promotion of advective mass transport via fluid flow being widely recognised as necessary for maximising nutrient delivery in bioreactor systems (Schonberg & Belfort 1987; O'Dea et al. 2012).

Additionally, assumptions of constant metabolite uptake in cells assuming excess nutrient availability have been used to simplify nutrient uptake kinetics (Shipley et al. 2011; Chapman et al. 2017). For example, nutrients such as oxygen are a key requirement for cellular respiration and oxidative phosphorylation, and growth factors such as fibroblast growth factors and vascular endothelial growth factors are essential for promoting cell signalling pathways for proliferation and differentiation into functional tissue. In practice, however, it is not always feasible to provide these nutrients in excess due to their significant cost and limited availability. Models that can account for a full range of inlet nutrient concentrations by utilising e.g. Michaelis–Menten uptake kinetics usually rely on computationally expensive numerical methods to solve the resulting transport problem, prohibiting detailed resolution of the metabolite field and the ability to comprehensively explore parameter space (Davidson, Ellis & Chaudhuri 2010; Mohebbi-Kalhori et al. 2012).

Given the long thin geometry inherent to HFMBs, a key challenge is ensuring all cells receive enough nutrient to grow and proliferate irrespective of their position along the bioreactor. Therefore, significant effort has been placed on determining optimal structural and operating parameter values such as inlet concentration, flow rate and membrane permeability, that ensure all cells receive sufficient nutrient to avoid the formation of hypoxic cell regions (Schonberg & Belfort 1987; Shipley et al. 2011; Chapman et al. 2014). These hypoxic regions are typically found at the furthest point from the lumen inlet i.e. near the ECS outlet, where nutrient has the furthest distance to travel. While mass transport limitations can be overcome by promoting advective mass transport via an increase in flow rate into the lumen or decrease in fluid flow resistance across the membrane wall, this comes at a cost of increased mechanical shear stress to the local cell environment. Furthermore, such solutions also do not eliminate the heterogeneity in nutrient delivery to cells, with cells receiving different amounts of nutrient depending on where they are located. This is a particular issue when seeking to engineer homogeneous tissue, as spatial differences in nutrient delivery will lead to non-uniform uptake and thus non-uniform tissue growth (Chapman et al. 2017).

In this paper, we systematically derive and analyse a formal reduced framework for modelling advection, diffusion and reaction of cellular metabolites in a single-fibre HFMB with nonlinear nutrient uptake and waste metabolite production. Specifically, we utilise physiologically relevant Michaelis–Menten reaction kinetics which are parameterised by cell-specific parameters and depend on the local metabolite concentration. We exploit the small thickness of the cell layer and the long thin bioreactor geometry to systematically reduce the model complexity, and introduce a non-trivial coordinate transformation to simplify the governing equations further. This approach generates significant savings in computational time compared with full numerical simulations, and facilitates quick and efficient exploration of parameter space. Importantly, our asymptotic reduction reveals mechanistic insights into nutrient transport mechanisms through HFMB systems, and identifies the existence of interfacial concentration boundary layers in the lumen and ECS which capture key transitions in dominant metabolite transport mechanisms. Similar to Dalwadi, Griffiths & Bruna (2015) and Köry et al. (2020), who considered heterogeneous membrane permeabilities for different membrane filtration applications, we consider heterogeneous membrane permeability distributions as a method to ensure uniform nutrient flux to cells via advection, which is the dominant transport mechanism across the membrane wall.

We begin in § 2 with a description of our model set-up for nutrient and waste metabolite transport in a single-fibre bioreactor system. In § 3 we non-dimensionalise, and systematically reduce and simplify our model. Utilising typical parameter values for HFMB operation from literature, we consider the transport of oxygen and carbon dioxide in a bioreactor culturing human endothelial cells. In § 4 we discuss our model results by exploring the impact of varying key operating and structural parameters, such as inlet/outlet pressures (relating to inlet fluid flux), membrane permeability and thickness, on metabolite transport. We then compare nutrient delivery in bioreactors with homogeneous and heterogeneous membrane permeabilities. Finally in § 5 we discuss the impact of this work on experimental design, and suggest future work.

2. Model description

2.1 Model set-up

We consider solute transport in a viscous Newtonian fluid through a long thin bioreactor of length $L^*$, with geometry illustrated in figure 2 and typical parameter values given in table 1. In cylindrical coordinates ($r^*, \theta ^*,z^*$) with the longitudinal axis coinciding with the cylinder axis, we define $\hat {\boldsymbol {e}}_r, \hat {\boldsymbol {e}}_\theta, \hat {\boldsymbol {e}}_z$ as unit vectors in the radial, azimuthal and axial directions, respectively. The bioreactor geometry is axisymmetric around the cylinder axis. Motivated by this, we neglect azimuthal components and derivatives, and simplify our problem to a two-dimensional axisymmetric set-up. We define a lumen of radius $R^*_{l0}$ surrounded by a porous membrane of uniform annular cross-section that extends the entire length of the bioreactor, with internal and external radii $R^*_{l0}$ and $R^*_{m0}$, respectively. A cell layer is seeded to the external membrane wall at $r^*=R^*_{m0}$ with an outer radius at $r^* = R^*_{c0}$. A rigid wall surrounds the bioreactor at $r^*=R^*_{e0}$, thus bounding our four key domains: the lumen ($l$), membrane ($m$), cell layer ($c$) and ECS ($e$). The nutrient and waste metabolite concentration in region $i$ are denoted by $c^*_i$ and $w^*_i$, respectively, and the fluid velocity vector is denoted by $\boldsymbol {U}^*_i = U^*_i\hat {\boldsymbol {e}}_r + W^*_i\hat {\boldsymbol {e}}_z$, where $i=l,m,c,e$. Although we present a model for a general flow field $\boldsymbol {U}^*_i$, for definiteness, when presenting results in § 4 we consider the flow field due to pressure-driven lubrication flow in the lumen and ECS coupled to Darcy flow in the membrane, the details of which are discussed in Appendix A.

Figure 2. Axial (a) and radial (b) cross-section of the bioreactor set-up. The fibre lumen ($l$) is enclosed by a porous fibre membrane wall ($m$) denoted by the grey region, with a cell layer ($c$) seeded to its outer surface represented by the pink region. The ECS ($e$) surrounds this, enclosed by an outer wall. Nutrient is introduced through the lumen inlet as indicated by the direction of the blue lines, denoting the different pathways available.

Table 1. Typical model parameters.

$^{{\dagger}}$ Parameters determined from internal unpublished experimental data for growth of human endothelial cells.

A nutrient-rich culture medium enters the bioreactor at the lumen inlet, driven by an axial pressure drop across the length of the lumen. Nutrient can transport along the fibre lumen and out through the lumen outlet, or across the membrane wall via advective and diffusive mass transport. Cells seeded at the outer membrane wall consume the nutrient and produce waste metabolite, which are transported away through the ECS along with excess nutrient, again via pressure-driven advection and diffusion. We consider the transport of both nutrient and waste products with typical parameter values outlined in table 1. Motivated by the longer time scales for cell proliferation relative to advection, diffusion, uptake and production, we neglect cell proliferation and focus on determining the quasi-steady nutrient distribution.

2.2 Governing equations

We model the concentration of nutrient $c^*_i$ in the lumen, membrane and ECS using the steady advection–diffusion equations. In the cell layer, we also introduce Michaelis–Menten reaction kinetics to account for nutrient uptake and waste metabolite production. We non-dimensionalise by setting

(2.1)\begin{equation} ( r^ * , z^ * ) = L^ * ( r, z), \quad (U^ * _i, W^ * _i ) = \bar{W}^ * (U_i, W_i), \quad (c^ * _i,w^ * _i) = C^ * _{in}(c_i, w_i). \end{equation}

Here, spatial coordinates are scaled with the length of the bioreactor $L^*$, and velocity components are scaled relative to the typical axial lumen fluid velocity $\bar {W}^*$, which is set by the typical operating conditions defined in table 1 and remains constant in our analysis. Concentrations are scaled with the inlet nutrient concentration in the lumen, $C^*_{in}$. The scaled nutrient transport governing equations in the lumen, membrane and ECS are given by

(2.2)\begin{equation} \mathcal{P}_i\boldsymbol{U}_i\boldsymbol{\cdot}\boldsymbol{\nabla} c_i = \nabla^{2}c_i, \quad i = l,m,e, \end{equation}

where $\mathcal {P}_l = \mathcal {P}_e = \bar {W}^*L^*/\mathcal {D}^*$ is the Péclet number representing the ratio between advective and diffusive nutrient transport effects in the lumen and ECS bulk flow regions, and $\mathcal {D}^*$ is the dimensional diffusivity of the nutrient in the culture medium. We assume that the diffusivity of nutrient in the porous membrane and cell layer regions is proportional to the volume fraction, $\phi _i$, of the culture medium (i.e. $\mathcal {D}^*_i = \phi _i \mathcal {D}^*$ for $i=m,c$). The effective Péclet number in each porous region is then $\mathcal {P}_i = \mathcal {P}_l/{\phi _i}$.

The governing equations in the cell layer are

(2.3)\begin{equation} \mathcal{P}_c\boldsymbol{U}_c\boldsymbol{\cdot}\boldsymbol{\nabla} c_c = \nabla^{2}c_c - \mathcal{R}(c_c), \end{equation}

where

(2.4a–c)\begin{equation} V_{max} = \frac{V^ * _{max}{L^ * }^{2}}{\phi_c C^ * _{in}\mathcal{D}^ * }, \quad K_M = \frac{K^ * _M}{C^ * _{in}}, \quad \mathcal{R} = \frac{V_{max} c_c}{K_M + c_c}. \end{equation}

Here, $V^*_{max}$ is the maximum uptake rate and $K^*_M$ is the Michaelis constant, which denotes the nutrient concentration at which the uptake rate is half of $V^*_{max}$. Therefore, $V_{max}$ denotes the ratio between reaction and diffusion time scales, and $K_M$ denotes the dimensionless Michaelis constant.

2.3 Boundary conditions

To determine the nutrient concentration in each domain, we require interfacial conditions and boundary conditions. At the lumen–membrane, membrane–cell and cell–ECS interfaces, respectively, we prescribe continuity of concentration and concentration flux

(2.5a)\begin{gather} c_l = c_m, (\mathcal{P}_l c_l\boldsymbol{U}_l - \boldsymbol{\nabla} c_l) \boldsymbol{\cdot} \hat{\boldsymbol{e}}_r = \phi_m(\mathcal{P}_l c_m\boldsymbol{U}_m - \phi_m\boldsymbol{\nabla} c_m) \boldsymbol{\cdot} \hat{\boldsymbol{e}}_r \quad \mbox{at} \ r = \varepsilon, \end{gather}

(2.5b)\begin{gather}c_m = c_c, \phi_m(\mathcal{P}_l c_m\boldsymbol{U}_m - \phi_m\boldsymbol{\nabla} c_m) \boldsymbol{\cdot} \hat{\boldsymbol{e}}_r = \phi_c(\mathcal{P}_l c_c\boldsymbol{U}_c - \phi_c\boldsymbol{\nabla} c_c) \boldsymbol{\cdot} \hat{\boldsymbol{e}}_r \quad \mbox{at} \ r = \varepsilon R_{m0}, \end{gather}

(2.5c)\begin{gather}c_c = c_e, \phi_c(\mathcal{P}_l c_c\boldsymbol{U}_c - \phi_c\boldsymbol{\nabla} c_c) \boldsymbol{\cdot} \hat{\boldsymbol{e}}_r = (\mathcal{P}_l c_e\boldsymbol{U}_e - \boldsymbol{\nabla} c_e) \boldsymbol{\cdot} \hat{\boldsymbol{e}}_r \quad \mbox{at}\ r = \varepsilon R_{c0}, \end{gather}

where $R_{i0} = R^*_{i0}/R^*_{l0}$ denotes the non-dimensional radius of domain $i$, and $\varepsilon = R^*_{l0}/L^*$ denotes the lumen aspect ratio. Exploiting axisymmetry, we prescribe no diffusive flux at the lumen centreline

(2.6)\begin{equation} \boldsymbol{\nabla} c_l\boldsymbol{\cdot} \hat{\boldsymbol{e}}_r = 0 \quad \mbox{at } r = 0, \end{equation}

and at the solid outer bioreactor wall we prescribe a no-flux condition

(2.7)\begin{equation} (\mathcal{P}_e c_e\boldsymbol{U}_e - \boldsymbol{\nabla} c_e)\boldsymbol{\cdot} \hat{\boldsymbol{e}}_r = 0 \quad \mbox{at} \ r = \varepsilon R_{e0}. \end{equation}

At each of the axial ends we prescribe no flux in the membrane and cell layer regions

(2.8)\begin{equation} (\mathcal{P}_l c_i\boldsymbol{U}_i - \phi_i\boldsymbol{\nabla} c_i)\boldsymbol{\cdot} \hat{\boldsymbol{e}}_z = 0 \quad \mbox{at} \ z = 0,1, \quad i = m,c. \end{equation}

Given the bioreactor operating conditions, we introduce fresh nutrient through the lumen inlet but not through the ECS inlet, which we capture via the following constant inlet concentration conditions:

(2.9)\begin{equation} c_l = 1, \quad c_e = 0 \quad \mbox{at} \ z = 0. \end{equation}

Finally, as solute exits the bioreactor via the lumen or ECS, we prescribe no diffusive flux

(2.10)\begin{equation} \boldsymbol{\nabla} c_i\boldsymbol{\cdot} \hat{\boldsymbol{e}}_z = 0 \quad \mbox{at} \ z = 1, \quad i = l,e. \end{equation}

We model the transport of waste metabolite $w_i$ in the lumen, membrane and ECS using the advection–diffusion equation (2.2), replacing $c_i$ with $w_i$, and $\mathcal {P}_i$ with $\mathcal {P}_{w,i}$. We define $\mathcal {P}_{w,l} = \mathcal {P}_{w,e} = \bar {W}^*L^*/\mathcal {D}^*_w$ as the Péclet numbers for waste metabolite transport in the lumen and ECS, where $\mathcal {D}^*_w$ denotes the diffusivity of the waste metabolite in the culture medium. We let $\mathcal {P}_{w,i} = \mathcal {P}_{w,l}/\phi _i$ denote the Péclet number for waste metabolite transport in the porous membrane and cell layer regions, where $i = m,c$. We assume that waste product is generated in the cell layer at a rate proportional to the nutrient uptake. The governing equation for waste transport in the cell layer is then

(2.11)\begin{equation} \mathcal{P}_{w,c}\boldsymbol{U}_c\boldsymbol{\cdot}\boldsymbol{\nabla} w_c = \nabla^{2}w_c + \nu\mathcal{R}(c_c), \end{equation}

where $\nu$ denotes the constant of proportionality between nutrient consumption and waste production. We also assume that no waste product is introduced to the system from the inlet ports, and prescribe a zero inlet concentration condition at the lumen and ECS inlet

(2.12)\begin{equation} w_i = 0 \quad \mbox{at} \ z = 0, \quad i = l,e. \end{equation}

The remaining boundary conditions for $w_i$ are equivalent to those in the nutrient transport problem, i.e. boundary conditions (2.5)–(2.8), (2.10) remain the same, replacing $c_i$ with $w_i$.

3. Model reduction and solution

We provide typical model parameter values in table 1. As referenced, we use parameter values from literature wherever possible.

3.1 Cell layer reduction

We investigate nutrient transport and uptake in the cell layer of thickness $\delta = \varepsilon (R_{c0} - R_{m0})$. As justified by the parameter values in table 1, we consider the physically relevant asymptotic limit $0 < \delta \ll \varepsilon \ll 1$. By scaling in to the $O (\delta )$ cell layer region and integrating the cell layer governing equations (2.3) across the thickness of the cell layer, we derive the following effective membrane–ECS interfacial condition

(3.1)\begin{equation} c_m = c_e, \quad \phi_m^2\frac{\partial c_m}{\partial r} - \frac{\partial c_e}{\partial r} = {-}\phi_c^2\bar{\mathcal{R}}({c_m}) \quad \mbox{at} \ r = \varepsilon R_{m0}, \end{equation}

where $\bar {\mathcal {R}} = \delta \mathcal {R}$, which we take to be an $O (1)$ quantity. Equation (3.1) states that the jump in radial nutrient flux across the cell layer is proportional to the uptake of nutrients in the cell layer, formally replacing (2.3), (2.5b) and (2.5c), and reducing our problem to consideration of nutrient transport through the lumen, membrane and ECS only. Full details of the cell layer reduction of the fluid flow and nutrient transport problem are found in Appendices A and B. Similarly, the effective membrane–ECS interfacial condition for waste metabolite transport becomes

(3.2)\begin{equation} w_m = w_e, \quad \phi_m^2\frac{\partial w_m}{\partial r} = \frac{\partial w_e}{\partial r} + \nu\phi_c^2\bar{\mathcal{R}}(c_e) \quad \mbox{at} \ r = \varepsilon R_{m0}, \end{equation}

formally replacing (2.5b), (2.5c) and (2.11).

3.2 Streamfunction–arclength reformulation

The long, thin bioreactor geometry ($0 < \varepsilon \ll 1$) will motivate the application of lubrication scalings in § 3.3. We note, however, that although the small radial distances (compared with axial distances) promotes radial diffusion compared with axial diffusion, the flow is sufficiently fast for advective mass transport to remain important and balance radial diffusion in the lumen and ECS in (2.2). The spatially two-dimensional nature of the advective transport terms means that numerical solution of (2.2) remains computationally expensive.

To circumvent this issue, we transform into a streamfunction–arclength ($\psi, s$) coordinate system as shown in figure 3. Instead of specifying the nutrient concentration at a radial $r$ and axial $z$ position in space, we parameterise space by which fluid streamline ($\psi$) a nutrient particle is travelling along, and how far along the streamline it has travelled ($s$). Without loss of generality, we impose $0 \leq s \leq 1$ in the lumen, $1 \leq s \leq 2$ in the membrane and $2 \leq s \leq 3$ in the ECS. We note this is similar to the approach in Cummings & Waters (2007) (for nutrient transport in a rotating bioreactor system). In our system, nutrients or waste metabolites are advected along streamlines (lines of constant streamfunction value), and therefore advective terms become spatially one-dimensional with respect to arclength. We consider the full transformation from a two-dimensional axisymmetric set-up into streamfunction–arclength coordinates for an arbitrary incompressible flow field in Appendix C. In these transformed coordinates, (2.2) becomes

(3.3)\begin{equation} \mathcal{P}_i\frac{\partial c_i}{\partial s} = \frac{\partial}{\partial\psi}\left(r^2 l_{\psi,i}\tilde{U}_i\frac{\partial c_i}{\partial \psi}\right) + \frac{\partial}{\partial s}\left(\frac{1}{l_{\psi,i}\tilde{U}_i}\frac{\partial c_i}{\partial s}\right) \quad \mbox{where}\ \tilde{U}_i = \sqrt{U_i^2 + W_i^2}, \end{equation}

and $l_{\psi,i}$ is the total length of the streamline corresponding to $\psi$ in region $i$.

Figure 3. Transformation from an axisymmetric coordinate system ($r,z$) on the left, to a streamfunction–arclength coordinate system ($\psi,s$) on the right. Purple lines denote equally spaced streamlines of constant $\psi$, with $\psi =0$ corresponding to the streamline along the symmetry line $r = 0$. Orange lines denote contours of constant arclength $s$ along streamlines, and black solid lines denote the boundary between each physical domain, where the lumen corresponds to $r \in (0, 1)$, $s \in (0,1)$, the membrane corresponds to $r \in (1,2)$, $s \in (1,2)$ and the ECS corresponds to $r \in (2,5)$, $s \in (2,3)$. Dotted black lines denote the critical streamlines $\psi _{i,crit}$ defined in Appendix C which divide where fluid enters or leaves. Note that $\psi _{l,crit}=\psi _{e,min}$ and $\psi _{l,max}=\psi _{e,crit}$. Appropriate boundary conditions are annotated by their corresponding equation number.

3.3 Long thin reduction

From (3.3), we see that the governing equations in streamfunction–arclength formulation represent a balance between advective mass transport along streamlines and diffusive mass transport in all directions. Given the similarity to the original transport (2.2), it may not yet be clear what we have gained by performing this transformation. However, we are now in a position to exploit the long thin bioreactor geometry to simplify our nutrient transport governing equations further.

3.3.1 Lumen and ECS

In the lumen and ECS, we introduce the lubrication scalings $r = \varepsilon r^\star, U_i = \varepsilon U_i^\star$, where $r^\star, U_i^\star = {O}(1)$. Since we have an $O (1)$ axial flow entering the lumen through a circular area of $O (\varepsilon )$ radius, fluid flux through the bioreactor is $O (\varepsilon ^2)$, and hence the streamfunction $\psi$ scales with $\varepsilon ^2$. This motivates scaling the streamfunction via $\psi = \varepsilon ^2\psi ^\star$. Utilising these scalings, dropping the superscript $\star$, and taking the limit $\varepsilon \to 0$, (3.3) at leading order is

(3.4)\begin{equation} \varepsilon^2\mathcal{P}_i\frac{\partial c_i}{\partial s} = \frac{\partial}{\partial\psi}\left(r^2l_{\psi,i} W_i\frac{\partial c_i}{\partial \psi}\right), \quad i = l,e, \end{equation}

where $W_i \geq 0$ throughout. We note that $\varepsilon ^{2}\mathcal {P}_i$ is the reduced Péclet number in region $i$ which we take to be an $O (1)$ quantity, motivated by the parameters in table 1, and we see a balance between advection along streamlines and diffusion between streamlines. This process has systematically reduced a two-dimensional elliptic partial differential equation (PDE) (3.3) to a nonlinear one-dimensional heat equation with non-constant coefficients (i.e. a parabolic PDE) which is much less computationally expensive to solve, where arclength is the time-like coordinate and streamfunction is the space-like coordinate.

3.3.2 Membrane

In the membrane, we again scale $r = \varepsilon r^\star, U_i = \varepsilon U_i^\star, \psi = \varepsilon ^2\psi ^\star$, where $r^\star, U_i^\star, \psi ^\star = {O}(1)$. However, the analysis in Appendix A reveals streamlines now align with the radial axis, which is an $O (\varepsilon )$ length scale compared with the bioreactor length. Hence, $l_{\psi,m}$ scales with $\varepsilon$ in the membrane, and so we write $l_{\psi,m} = \varepsilon l_{\psi,m}^\star$, where $l^\star _{\psi,m} = {O}(1)$. Utilising the velocity scalings in the membrane (Appendix A), we find

(3.5)\begin{equation} \tilde{U}_m = \sqrt{\varepsilon^2U^{{\star} 2}_m + \varepsilon^4W^{{\star} 2}_m} \sim \varepsilon\tilde{U}_m^\star. \end{equation}

Again dropping the superscript $\star$ and taking the limit $\varepsilon \to 0$, the leading-order transport governing equation in the membrane (3.3) is

(3.6)\begin{equation} \varepsilon^2\mathcal{P}_m\frac{\partial c_m}{\partial s} = \frac{\partial}{\partial s}\left(\frac{1}{l_{\psi,m}U_m}\frac{\partial c_m}{\partial s}\right), \end{equation}

where the reduced Péclet number in the membrane $\varepsilon ^{2}\mathcal {P}_m = \varepsilon ^{2}\mathcal {P}_l/\phi _m = {O}(1)$ (noting that $U_m > 0$ throughout). Equation (3.6) is now an ordinary differential equation (ODE), which can be solved to determine an analytical solution for concentration in the membrane.

3.3.3 Interfacial and boundary conditions

At leading order in $\varepsilon$, the transformed lumen–membrane and membrane–ECS interfacial conditions (2.5a), (3.1) in streamfunction–arclength coordinates become

(3.7a)\begin{gather} c_l = c_m, \quad \frac{\partial c_m}{\partial s} = 0 \quad \mbox{at} \ s = 1, \psi \in [\psi_{l,crit}, \psi_{l,max}), \end{gather}

(3.7b)\begin{gather}c_m = c_e, \quad \frac{\partial c_m}{\partial s} = {-}\frac{\varepsilon l_{\psi,m}\phi_c^2\bar{\mathcal{R}}(c_m)}{\phi_m^2} \quad \mbox{at} \ s = 2, \psi \in (\psi_{e,min}, \psi_{e,crit}], \end{gather}

where we take $\varepsilon \bar {\mathcal {R}}$ to be an $O (1)$ quantity. Here, $\psi _{i,min}, \psi _{i,max}$ denote the minimum and maximum streamline values in domain $i$, and $\psi _{i,crit}$ denotes the critical streamline that divides which inlet or outlet fluid enters and leaves from domain $i$ as shown in figure 3.

The lumen symmetry condition (2.6) and the no-flux condition at the outer wall (2.7) become

(3.8)\begin{equation} \frac{\partial c_i}{\partial\psi} = 0 \quad \mbox{at} \ \psi = 0, \psi_{e,max}, \end{equation}

and the inlet concentration conditions (2.9) remain the same

(3.9a)\begin{gather} c_l = 1 \quad \mbox{at} \ s = 0, \end{gather}

(3.9b)\begin{gather}c_e = 0 \quad \mbox{at} \ s = 2, \psi \in (\psi_{e,crit}, \psi_{e,max}]. \end{gather}

Finally, the outlet no-diffusive-flux condition (2.10) becomes

(3.10)\begin{equation} \frac{\partial c_i}{\partial s} = 0 \quad \mbox{at} \ s = 1, \psi\in (0, \psi_{l,crit}] \quad \mbox{and} \quad s = 3. \end{equation}

In Appendices C and D, we consider the transformation and reduction of the full governing equations and boundary conditions for nutrient and waste metabolite transport.

3.4 Model solution

The resulting leading-order governing equations (3.4) and (3.6) are coupled via the appropriately transformed continuity of concentration and flux conditions at the lumen–membrane and membrane–ECS interfaces (3.7). Up until now, analysis has been presented for a general flow field. However, to solve (3.4), (3.6) for the solute concentration, we require a solution to the fluid velocities in each domain, $U_i$ and $W_i$.

3.4.1 Solution of fluid flow problem

The particular flow problem we solve for is detailed in Appendix A, and shown in figure 4. In summary, we consider slow, viscous lubrication flow of an incompressible Newtonian fluid in the lumen and ECS, driven by an axial pressure drop between the inlet fluid pressure, $P_{i,in}$, and the outlet fluid pressure, $P_{i,out}$, where $i=l,e$. Flow in the lumen and ECS is coupled to incompressible Newtonian Darcy flow across the porous membrane wall of permeability $\kappa$. Utilising the analysis outlined in Appendix A, we determine the fluid velocities in each region of the bioreactor ($U_i, W_i$) in terms of fluid pressures ($P_i$) (see (A12) and (A14a,b)), and substitute into the transport governing equations (3.4) and (3.6).

Figure 4. Canonical flow set-up considered in our model analysis. Slow, viscous, Newtonian, pressure-driven lubrication flow is prescribed in the lumen and ECS, coupled to Darcy flow across the membrane wall. Flow enters the lumen and ECS inlets at $z=0$ with fluid pressures $P_{l,in}, P_{e,in}$, respectively. Flow exits the lumen and ECS at $z=1$ with fluid pressures $P_{l,out}, P_{e,out},$ respectively.

3.4.2 Solution of transport problem

To solve (3.4), we use the method of lines with a central-difference discretisation for $\psi$, essentially solving a system of coupled ODEs with $s$ as a time-like variable using MATLAB's in-built ode15s solver. Initial conditions for the ODE solvers are provided by the inlet concentration condition (3.9a) in the lumen, and the continuity of concentration conditions (3.7b) and (3.9b) in the ECS.

We can solve the reduced equation (3.6) exactly by substituting in our analytic solution for membrane fluid flux (A14a,b) and imposing continuity of concentration at the lumen–membrane interface (3.7a) and the effective uptake condition at the membrane–ECS interface (3.7b). Therefore, nutrient concentration in the membrane is given by

(3.11)\begin{equation} c_m = 1 - \mathcal{A}((s-1)(R_{m0}-1) + 1)^{-\xi}, \quad \xi = \frac{\varepsilon^2\mathcal{P}_l\kappa(P_e - P_l)}{\phi_m^2\ln R_{m0}}, \quad \mathcal{A} = \frac{b + \sqrt{b^2 + 4ac}}{2a}, \end{equation}

where $P_l, P_e$ are the fluid pressures in the lumen and ECS, respectively, $s\in (1,2)$ in the membrane parameterises the position along streamlines, and

(3.12a–c)\begin{equation} a = \phi_m^2\xi R_{m0}^{-(2\xi + 1)}, \quad b = \phi_m^2\xi R_{m0}^{-(\xi + 1)}(K_M + 1) - \varepsilon\delta\phi_c^2V_{max} R_{m0}^{-\xi}, \quad c = \varepsilon\delta\phi_c^2V_{max}. \end{equation}

From (3.11), we see that nutrient concentration in the membrane is directly proportional to the transmural pressure drop across the wall, $P_l - P_e$. This transmural flow is a key driving force for nutrient transport in the membrane. Specifically, axial variations in nutrient concentration in the membrane arise from the axial dependency coming from the fluid pressure variations. We note that, while we solve our system of equations in streamfunction–arclength formulation, solutions presented in § 4 have been transformed back into two-dimensional axisymmetric coordinates.

We present full details of the metabolite transport solution method in Appendix E. We note that the same procedure outlined in § 3 can also be followed to mathematically reduce and solve the governing equations for waste metabolite transport in each domain, the details of which we discuss in Appendices C–E.

4. Results

4.1 Typical operating conditions

We first consider steady fluid flow and nutrient transport under typical bioreactor operating conditions (summarised in table 1) as shown in figure 5(a). In the lumen, the transformed governing equation for metabolite transport (3.4) is a nonlinear diffusion equation. Solving subject to constant inlet concentration (3.9a), the lumen concentration is constant along streamlines and therefore equal to the prescribed inlet concentration. Advective transport of nutrient across the membrane wall to the cells seeded at $r=R_{m0}$ is driven by a transmural pressure drop, which varies axially due to the axial variation in fluid pressure along the lumen and ECS, as shown in (3.11). Therefore, we expect to see a higher nutrient concentration delivered to the cells where the radial pressure drop across the membrane wall is highest. Any leftover nutrient is then washed away by flow in the ECS.

Figure 5. Nutrient concentration profiles across a cross-section of the bioreactor set-up. (a) Nutrient concentration profile under typical experimental conditions as outlined in table 1. Black lines denote the interfaces between the lumen ($r \in (0,1)$), the membrane ($r \in (1,2)$) and the ECS ($r \in (2,5)$). Fresh nutrient enters through the lumen inlet at $z=0$ and transports either along the length of the lumen or across the membrane wall to be consumed by cells seeded on the outer membrane surface at $r = 2$. Excess nutrient is washed away by flow in the ECS. Solid white lines denote streamlines beginning at the lumen inlet and (b) nutrient concentration profile at the outer membrane wall under varying axial pressure drops holding $P_{l,in}, P_{l,out}, P_{e,in}$ constant. Axial heterogeneity in nutrient concentration increases as axial variations in the transmural pressure drop increase. In particular, at a given axial position, the nutrient concentration delivered to the cells seeded at the outer membrane wall decreases when the transmural pressure drop decreases.

Focusing on figure 5(b) (left), we see that nutrient transport to the outer membrane interface under typical operating conditions is relatively uniform due to a pressure drop that remains relatively constant across the membrane. If we fix the lumen and ECS inlet pressures and vary the transmural pressure drop (e.g. by varying the ECS outlet pressure, $P_{e,out}$), we observe a strong dependence of nutrient delivery on the transmural pressure drop in figure 5(b). Towards the bioreactor outlet where the transmural pressure drop decreases, advective mass transport across the membrane wall decreases, and hence nutrient delivery also decreases. This spatial heterogeneity of nutrient delivery means that cells are receiving unequal amounts of nutrient depending on where they are seeded along the length of the membrane. This is of particular concern in engineering homogeneous tissue constructs, which have a uniform nutrient requirement, since heterogeneous nutrient delivery will lead to heterogeneous tissue growth. We see that to achieve uniform nutrient delivery across the membrane wall, we require a uniform transmural pressure drop irrespective of the given operating conditions.

4.2 Homogeneous vs. heterogeneous permeability distributions

One way to generate uniform nutrient delivery via advection is to use a membrane with a heterogeneous permeability. Advective mass transport is a function of the local permeability, $\kappa$, and transmural pressure drop, $P_l - P_e$, given by (3.11). At positions where advective nutrient flux across the membrane wall is smaller due to a smaller transmural pressure, introducing a larger permeability will reduce the resistance for advective transport across these regions, thus raising the flux (and vice versa for a smaller permeability). We can derive an analytic representation of this favourable permeability distribution from (3.11), where we seek a membrane permeability that exactly counteracts transmural pressure

(4.1)\begin{equation} \kappa(z) = \frac{1}{P_l(z) - P_e(z)}. \end{equation}

Figure 6(a) shows the concentration of nutrient delivered to cells seeded at the outer membrane wall as a function of the operating conditions, comparing a homogeneous membrane permeability to the heterogeneous permeability given in (4.1) under typical bioreactor operating conditions. Utilising the heterogeneous permeability distribution exactly cancels out the axial variations in $c_m$ from the transmural pressure drop in figure 6(a), resulting in a membrane concentration profile that is axially uniform. For the considered model parameters, the heterogeneous permeability distribution also leads to a larger average nutrient concentration. As the desired heterogeneous permeability distribution is a function of the lumen and ECS fluid pressures, any variations in $P_{l,in}$ and $P_{e,in}$ will alter the appropriate permeability distribution. However, because this will cancel out the fluid pressure terms in our solution to $c_m$, our concentration profile remains unchanged by variations in these parameters.

Figure 6. Comparing the nutrient concentration delivered to the outer membrane wall with a homogeneous vs a heterogeneous membrane permeability under varying structural and operating conditions. (a) Varying permeability, $\kappa$ (left), lumen inlet pressure, $P_{l,in}$ (middle), and ECS inlet pressure $P_{e,in}$ (right). Coloured lines denote the concentration profile at $R_{m0}$ with a homogeneous membrane permeability, and the black dotted line is the concentration profile with a heterogeneous membrane permeability. Note that the membrane concentration with a heterogeneous permeability is invariant to changes in the selected parameters and (b) Comparing the average nutrient concentration at the outer membrane wall, $\bar {c}|_{R_{m0}}$ with a homogeneous vs ideal heterogeneous permeability under varying membrane thicknesses $R_{m0}$ and reduced Péclet numbers $\varepsilon ^2\mathcal {P}_l$. The solid white lines denote contours of constant nutrient concentration. The average concentration of nutrient delivered to cells at the outer membrane wall is greater in the heterogeneous membrane at each corresponding membrane thickness and reduced Péclet number.

Our theoretically favourable solution does, however, depend on both the membrane thickness and the lumen Péclet number. Figure 6(b) shows how increasing the membrane thickness decreases the nutrient concentration at the outer wall as the radial membrane fluid flux decreases, thus decreasing advective mass transport. As might be expected, an increase in the lumen Péclet number increases advective mass transport, therefore increasing nutrient delivery across the membrane wall. In both cases, the optimal nutrient delivery problem delivers a much higher nutrient concentration to the outer wall compared with a membrane with constant membrane permeability, suggesting that in all scenarios, a heterogeneous permeability distribution is of significant benefit for both efficient and uniform nutrient transport.

4.3 Waste metabolite transport

Using the same methodology as § 4.2, we can model waste metabolite transport through our bioreactor system. This is of particular interest in understanding the intricate relationship between maximising nutrient delivery and waste removal, and how both factors can be affected by changes in structural and operating parameters.

Figure 7 shows how the waste metabolite concentration at the outer wall varies under different operating conditions. Similarly to nutrient concentration in figure 6, we see that any change in parameter value that increases advective mass transport (either by increasing the pressure drop or decreasing resistance to flow across the membrane) increases flux across and out of the membrane. For the same reasons that this increases nutrient concentration delivery to the outer membrane wall, this increase in advective mass transport also washes away any waste metabolite produced more efficiently, thus decreasing the waste metabolite concentration at the outer wall. Interestingly, we see that increasing the inlet ECS pressure, corresponding to an increase in flow rate, actually increases waste concentration at the outer membrane wall compared with lower ECS pressures. As the inlet ECS pressure increases towards the inlet lumen pressure, the pressure drop across the membrane decreases and thus fluid flux across the membrane wall also decreases. Therefore, waste product is removed less easily. These results suggest that, in the parameter regimes considered, the most important mechanism for waste product removal in HFMBs is maximising fluid flow across the membrane wall rather than increasing the axial flow rate through the ECS.

Figure 7. Waste metabolite concentration at the outer membrane wall under varying lumen inlet pressure (a), ECS inlet pressure (b) and membrane wall thickness (c). The red line denotes the average concentration with a homogeneous permeability given in table 1 with the shaded red region denoting the maximum and minimum concentrations. This shaded region is very thin in the right-most panel. The blue line denotes the concentration under the heterogeneous permeability distribution.

Similarly to before, we note that the specific heterogeneous permeability distribution we derive results in both a uniform production of waste metabolite, and an overall reduction in local waste concentration compared with a homogeneous permeability for the majority of parameter values. We conclude that the heterogeneous membrane allows for the delivery of a uniform nutrient concentration at higher concentrations compared with the homogeneous case under the same operating conditions, and washes away waste metabolites more efficiently.

5. Discussion

In this paper, we present a two-dimensional axisymmetric model for nutrient and waste transport through a hollow fibre membrane bioreactor. Utilising a coordinate transformation into a streamfunction–arclength formulation and exploiting the long thin geometry, we use an asymptotic analysis to systematically reduce the dimensionality of our governing equations and derive a reduced-order system which we solve efficiently using the method of lines. We interpret our results to understand how the membrane permeability can be manipulated to achieve uniform nutrient transport to the cell layer, as well as to understand how to tune physical operating and geometric parameters to improve the balance between nutrient delivery and waste product removal.

Our systematic derivation of the reduced model means that we can solve for the nutrient transport problem with greater speed and efficiency compared with full multiphysics simulations in software packages such as COMSOL. A full simulation of nutrient transport through our bioreactor design under typical operating conditions, using a fine mesh necessary to resolve variations in concentration in the radial direction, takes around seven minutes to compute in COMSOL on a standard desktop computer. Alternatively, a full simulation of our reduced model takes a fraction of a second to solve for the fluid flow problem, and around 1 minute to solve for the nutrient transport problem, for a comparative mesh size. Better still, we are able to solve for the nutrient concentration profile through the membrane and at the outer membrane wall using our analytic solution as per (3.11) in a fraction of a second. Our analysis also provides us with a full mechanistic understanding of the dominant nutrient transport mechanisms across the membrane and the transitions between different regions, which reveals the impact of system parameters on the model solution. This model reduction in turn also speeds up sensitivity analysis considerably.

For example, our model reduction allows us to understand the role of heterogeneous permeability distributions in achieving uniform nutrient delivery. From (E5), we derive an explicit relationship between the nutrient concentration in the membrane and the membrane permeability, as well as an explicit understanding of how axial variations in the membrane flux leads to axial variations in nutrient delivery. This means that we derive an explicit formula for the optimal membrane permeability that cancels out axial variations and leads to uniform nutrient delivery to the cell site. This is of particular importance for informing experimental design, whereby appropriate mathematical modelling can help inform the manufacture of membranes with the heterogeneous permeability distribution under desired operating conditions. For example, by changing the applied voltage or adjusting the distance between the spinneret and collector in the electrospinning of hollow fibre membranes, one can generate a gradient in fibre diameter and density, which in turn affects the membrane's permeability across different regions (Forgie et al. 2023). We note, however, that the spatial resolution is limited by the voltage, geometry and polymer properties in the electrospinning set-up, placing limits on the achievable permeability range.

This ‘optimal’ heterogeneous permeability distribution results in a uniform delivery of nutrient to a spatially uniform cell layer of uniform cell type and seeding density. In the case of a heterogeneous cell seeding density, or in a set-up containing several different cell types, the effect of a spatially varying nutrient demand will impact the optimal nutrient delivery, which in turn changes the corresponding heterogeneous permeability distribution. Given an axial dependency of nutrient uptake in the cell layer, it is straightforward to use our methodology to determine a heterogeneous permeability that depends on advection across the membrane wall and on the spatially varying reaction kinetics. For example, this could be theoretically achievable by setting the membrane concentration at the outer membrane wall in (3.11) to be constant, and rearranging to determine an explicit relationship between the permeability $\kappa$, and the axial dependency arising from the reaction kinetic parameters and the transmural pressure drop.

In our analysis of the nutrient and waste transport problem, it is clear that promoting advective mass transport not only improves nutrient delivery to the cell site, but also promotes the removal of waste cell products. However, in practice, an increase in advective mass transport comes at a cost of increased mechanical stress imparted on the cells. Cells are particularly sensitive to their local mechanical environment, which is of importance in regulating cell morphology and for mechanotransduction signalling pathways. Excessive shear stresses can damage and even detach cells seeded to the membrane wall in an experimental set-up, presenting a difficult trade-off between optimising the local chemical and mechanical cell environment. In our current model, the indirect impact of mechanical stress on nutrient uptake is not accounted for. However the derivation of the effective membrane–ECS interfacial condition holds for a general nutrient uptake reaction term, $\mathcal {R}$, and therefore can be used when introducing reaction kinetics that depend on the local mechanical stress. Furthermore, the semi-analytic nature of the framework presented in this study presents a quick and efficient method for implementing optimisation methods to determine the optimal system of structural and operating conditions that would lead to the desired trade-off between nutrient delivery and mechanical disruption.

The speed-ups offered by a reduced modelling approach are of particular importance in extending this work beyond traditional hollow fibre bioreactors to consider faster, more physiologically relevant fluid flow problems. The framework presented in this paper can be used alongside any underlying two-dimensional fluid velocity field provided that we can determine the accompanying streamfunction, either analytically or numerically. This will be of use in future work where we relax the assumption of Stokes flow, and introduce alternative flow profiles applicable in a range of experimental designs (Booth 2023).

The modelling approaches in this paper present a robust and flexible framework for the analysis of nutrient transport in hollow fibre membrane bioreactor systems. With these tools, we develop a mechanistic understanding of the system's response to varying structural and operating parameters, and the impact this has on the nutrient delivery/waste removal problem. This lays the groundwork for a plethora of mathematical extensions that more accurately describe the local chemical and mechanical environment of the cell layer, and provide us the tools to quickly and accurately traverse a large parameter space to determine the optimal bioreactor set-up for uniform nutrient delivery to the cells.

6. Conclusion

We have developed a mathematical model to describe steady fluid flow and cell metabolite transport in a single-fibre HFMB, and have considered the impact of varying structural and operating parameters on nutrient delivery and waste metabolite removal. We have developed the framework to accommodate any two-dimensional flow, and we showcase our methodology by considering the specific problem of lubrication flow in the lumen and ECS coupled to Darcy flow across the membrane. Having investigated the dependency of nutrient transport on the transmural pressure drop, we have determined a theoretically ideal heterogeneous membrane permeability distribution that facilitates uniform nutrient delivery to cells for engineering homogeneous cell populations of uniform size and density. This analysis provides a robust framework for understanding and optimising HFMB design and operation, valid across a wide range of possible operating Péclet numbers and limiting nutrient concentrations. Given the computational time savings relative to typical numerical simulations, the analysis in this paper can be utilised by experimentalists to quickly and efficiently understand and optimise experimental bioreactor design for applications in tissue engineering research.