A strategy to determine operating parameters in tissue engineering hollow fiber bioreactors

The development of tissue engineering hollow fiber bioreactors (HFB) requires the optimal design of the geometry and operation parameters of the system. This article provides a strategy for specifying operating conditions for the system based on mathematical models of oxygen delivery to the cell population. Analytical and numerical solutions of these models are developed based on Michaelis–Menten kinetics. Depending on the minimum oxygen concentration required to culture a functional cell population, together with the oxygen uptake kinetics, the strategy dictates the model needed to describe mass transport so that the operating conditions can be defined. If cmin ≫ Km we capture oxygen uptake using zero-order kinetics and proceed analytically. This enables operating equations to be developed that allow the user to choose the medium flow rate, lumen length, and ECS depth to provide a prescribed value of cmin. When , we use numerical techniques to solve full Michaelis–Menten kinetics and present operating data for the bioreactor. The strategy presented utilizes both analytical and numerical approaches and can be applied to any cell type with known oxygen transport properties and uptake kinetics.


A Analytical Model Reduction
The full modelling equations throughout the module are given by ∂c l ∂t +∇·(uc l ) = D l ∇ 2 c l in the lumen, ∂c m ∂t = D m ∇ 2 c m in the membrane, ∂c e ∂t = D e ∇ 2 c e −R (c e ) in the ECS (1) where the reaction term is given by On the lumen/membrane and membrane/ECS boundaries we prescribe continuity of concentration and flux, so that c l = c w and D l ∇c l · n = D w ∇c w · n on the lumen/membrane boundary, c w = c e and D w ∇c w · n = D e ∇c e · n on the membrane/ECS boundary, where n is the unit outward pointing normal to the relevant surface. Finally, we impose a prescribed concentration c in at the lumen inlet, and no diffusive flux of concentration out of the outer ECS boundary, c l = c in on z = 0, and D e ∇c e · n = 0 on the outer ECS boundary.
We mathematically reduce the full system of equations given by equations (1)-(5) to determine analytical expressions for the oxygen concentrations in each of the lumen, membrane and ECS.
First of all, we assume that the fibre is positioned symmetrically in the extra-lumen space and move to a radially symmetric setup described in cylindrical polar coordinates (defined by x = r cos θ, y = r sin θ and z = z), neglecting θ-dependence so that c = c (r, z).
We non-dimensionalise the system above to reduce the number of parameters and to estimate the relative importance of the various terms. We set where d is the lumen radius, L is the length of the lumen, U is the mean velocity in the lumen and c in is the lumen inlet oxygen concentration. From here onwards we drop the 'hats' on the dimensionless variables. The system now becomes Pe ∂c l ∂t + 2 1 − r 2 ∂c l ∂z = 1 r ∂ ∂r r ∂c l ∂r + 2 ∂c l ∂z in 0 < r < 1, PeD l D w with boundary conditions c l = c w and ∂c l ∂r where all parameters are defined in the main paper.

A.1 Solution in the Membrane and ECS
Next we solve the system described by equations (7)-(13) in steady state to predict the oxygen distribution throughout the module. We integrate the steady-state versions of equations (8) and (9) for the concentration in the wall and ECS to give where A (z), B (z), C (z) and D (z) are arbitrary functions of z. Inputting the no-flux boundary condition out of the ECS outlet (from equation (12)), together with the continuity of concentration and flux boundary conditions on r = R w /d (equation (11)) yields where B (z) remains unknown. To determine B, and the solution in the lumen, we must now solve equation (7) (in steady-state), subject to the boundary conditions where γ is defined by and we also require that ∂c l ∂r is bounded at r = 0.

A.2 Solution in the Lumen
We perform a change of variables to give a homogeneous derivative boundary condition on r = 1. We let upon which the system in the lumen becomes ∂c ∂r finite as r → 0,

A.2.1 Solution of the Homogeneous Counterpart to (20)-(22)
First of all, we consider the homogeneous counterpart to (20) given by and solve via separation of variables by seeking a solution of the form c (r, z) = f (r) g (z). This yields The LHS of equation (25) is dependent on r only, whilst the RHS is dependent on z only. Therefore both sides must constant, and indeed this constant must be negative (the concentration must decay as z increases from the inlet at z = 0). Defining this constant to be −λ 2 , the system for f (r) becomes with boundary conditions df dr (1) = 0 and df dr finite as r → 0.
This is a Sturm-Liouville problem -for further details please refer to Abramowitz & Stegun (1965). To solve for f (r), we perform two transformations which transform (26) to which is known as the Kummer equation. The Kummer equation is normally written in the form here we have ν = 1 and µ = 1 2 − λ 4 . The Kummer equation has two linearly independent solutions, denoted KummerM (µ, ν, x) and KummerU (µ, ν, x) (both hypergeometric functions), and so our solution here can be written as a linear sum of the two, where E and F are integration constants to be determined. A brief summary of the properties of the KummerM and KummerU functions required here is given below; for further details please refer to Abramowitz & Stegun (1965).
The requirement that df /dr is finite as r → 0 yields F = 0 (to avoid a 1/ √ X singularity as X → 0) -this uses property 1. The boundary condition at r = 1 from (27) provides the following eigenvalue equation for λ, The determination of this equation uses property 2. above to evaluate the derivative, and 3. to reduce oscillations in the eigenvalue equation. Equation (32) has infinitely many solutions for λ which may be determined by solving Equation (32) numerically; we denote them by λ n , where n = 0 . . . ∞. Each eigenvalue λ n corresponds to a different solution for f (an eigenfunction), which we now denote by f n . These functions f n are given by where Sturm-Liouville theory dictates that the coefficients E n may be fixed used the normalization condition 1 r=0 s (r) f 2 n (r) dr = 1, for each n = 0 . . . ∞, where s (r) = r 1 − r 2 are the weight functions of the Sturm-Liouville equation (26). Therefore,

A.2.2 Solution of the Non-Homogeneous Problem (20)-(22)
Now we return to the non-homogeneous problem given by (20)- (22), and pick up the non-homogeneous term γ/r 1 − r 2 using the radial solution, f . Sturm-Liouville theory states that this non-homogeneous term can be expressed as a linear sum of the eigenfunctions f n , which act as basis functions (analogous to a Fourier series expansion). Therefore, where the F n are constants given by Next we express the solution for c as a linear sum over all the eigenvalues and eigenfunctions (admissible because the equations (20)-(23) are linear), where the g n (z) are undetermined. Substituting the two sum expressions (36) and (38) into the Equation (20) for c and using Equation (26) to simplify the derivative terms in f n yields the following ordinary differential equation for each g n (n = 0 . . . ∞) 2Pe r dg n dz + λ 2 n g n (z) = F n .
The solution to (39) can be written in the form where the integration constants G n are undetermined. Indeed, substituting for c from Equation (38) Expressing g n (0) as a Sturm-Liouville sum yields where the expression for g n given by Equation (40) means that Therefore the constants G n are given by The final solution for c l is summarized in the main paper.