Analysis of oxygen transport in microfluidic bioreactors for cell culture and organ‐on‐chip applications

In the recent decade, development of microfluidic bioreactors and organ‐on‐chip platforms for drug screening and disease modeling has been rising significantly. Prediction of oxygen level within the microfluidic bioreactors to create physiologically relevant oxygen tension for realistic cellular behavior is of critical importance. This article presented an analytical method to calculate oxygen tension in microchannel parallel plate bioreactors. Two‐dimensional convection‐diffusion equation was solved analytically with considering diffusion in two directions. Cellular oxygen consumption was assumed to follow Michaelis‐Menten kinetics. Effects of oxygenator design, gas permeability of the microfluidic channels, as well as flow rate and cell density on the oxygen distribution within the bioreactor were examined. In addition, a mathematical model was developed to predict oxygen tension in a fluidic circuit containing two interconnected bioreactors with and without recirculation of the media for the culture of hepatocytes and cardiomyocytes. The oxygenators were used to maintain normoxia in the bioreactors independently. The model allowed for the prediction and independent modulation of oxygen tension in the two bioreactors through changing the length of the oxygenators. In overall, the obtained results provide critical insights required for the design and operation of microfluidic bioreactors with desired levels of oxygen tension.


INTRODUCTION
Discovery of new drugs is generally very costly and requires long-term research studies. Drug development processes include biological cell-based tests, animal tests, and clinical trials. Eventually, an average of less than 11% of animal tests on human samples provide successful outcome to get approval from Food and Drug Administration. 1,2 Nonclinical/clinical safety is the cause of more than half of the failures while the efficacy is responsible for more than 10% of the failures. 3,4 The main reason for high attrition rates of drug candidates is that the animal model responds differently compared to its human counterparts. 5 Besides, in vitro cell-based drug tests are usually performed using either the immortalized cell lines or only limited cell types. Thus, such tests cannot fully reflect the body's response. 6 In this way, there is a gap between the cell-based/animal studies and the human-based trials. To address such issues, microfluidic devices and organ-on-chip (OoC) systems have been recently developed to mimic the cellular microenvironment for more realistic predictions of cell-drug responses at in vitro conditions. 1 Microfluidic OoC systems are mainly comprised of interconnected bioreactors containing cells to recapitulate the minimal organs functions. 7 Various parameters such as cellular matrix stiffness, pH, oxygen tension, and mechanical forces such as flow-induced shear stress should be deliberately tuned within the bioreactor to mimic the physiological cellular microenvironment conditions.
Oxygen tension within the microfluidic bioreactors is a critical parameter that has significant effects on the regulation of cellular behaviors. 8 In particular, lack of sufficient oxygen (hypoxia) can cause changes in physiological pathways, 9 cellular metabolism, 10 and remodeling of tissue. 11 In addition, high oxygen tension (hyperoxia) may result in the generation of elevated level of reactive oxygen species in pathological pathways. 12 Thus, the balance of oxygen availability and oxygen consumption is necessary to maintain normoxia and to prevent hyperoxia and hypoxia or even complete absence of oxygen (anoxia). 13,14 Importantly, in multiple-OoC systems, two or more bioreactors, which might require different levels of oxygen tension, are interconnected. In this way, the prediction and control of oxygen tension in each bioreactor is a critical design consideration. 15 Delivery of sufficient oxygen to microfluidic bioreactors mainly occur through (i) oxygen diffusion from the permeable walls of the fluidic system or (ii) an oxygenator that is mainly made from a thin oxygen-permeable membrane implemented in the wall of the fluidic system or the bioreactor. The former oxygen transport mainly occurs for bioreactors made in PDMS that has high oxygen permeability. The latter oxygen delivery is designed for the microfluidic systems require localized oxygen transport (gradient of oxygen tension) or within the microfluidic systems made from low-oxygen permeability materials such as certain types of thermoplastics including PMMA and PC. For such devices, the oxygenation is provided through an oxygen-permeable membrane that can be made from high-oxygen-permeable materials such as PDMS or TPU. 16 The easiest way to control oxygen tension in a bioreactor is through the modulation of the flow rate. However, changing the flow rate could affect the flow-induced shear stresses experienced by the cultured cells, pressure drop, and the transport characteristics of biomolecular species that consequently limits the practical domain of the flow rate. 17,18 Moreover, by changing the cell density and microchannel dimensions, the oxygen tension within the microchannel could be modulated. 17,19 Thus, to design a bioreactor with accurate levels of desired oxygen tension, effective parameters should be taken into consideration.
In this paper, a mathematical model to determine oxygen transport within microchannel parallel plate bioreactor with and without an internal oxygen-permeable membrane was developed. Oxygen diffusion in both directions, including parallel to the flow direction and perpendicular to the flow direction due to the cellular oxygen uptake were considered. Effects of various parameters such as flow velocity, cell density, and the microchannel length on the oxygen tension of the bioreactors with and without the use of an oxygenator were investigated. Two modes of operation for bioreactors were studied including the single-pass perfusion and the closed-loop recirculation perfusion. The results provide critical insights required for the design and operation of microfluidic bioreactors with desired levels of oxygen tension. In overall, the proposed method could be employed to evaluate oxygen tension in bioreactors for accurate oxygen modulation in microfluidic cell culture applications.

MATERIALS AND METHODS
The aim of this study is to determine the spatial oxygen concentration within microchannel bioreactors to create normoxia conditions for cell culture applications. Figure 1 shows the bioreactor design with and without internal oxygen-permeable membrane for oxygenation. These bioreactors are modeled as two-dimensional rectangular microchannels with pressure-driven flow. Cell culture medium is considered as the perfusion liquid. Cellular monolayer (rat hepatocytes) is assumed to adhere and cover the bottom of the channel. The bioreactor consists of a microchannel with depth of 100 μm and length of 20 mm. The geometric ratio of the channel is defined as γ = H/L, where H is the channel height and L is the channel length.
The flow profile is solved explicitly using the steady state Navier-Stokes equation with incompressible flow assumption. The equations of continuity and steady state conservation momentum are u refers to the average velocity P is the pressure, ρ is the density, and μ is the viscosity. Application of no-slip condition over the wall and uniform specified axial velocity at the entrance following velocity profile is obtained: The distribution of oxygen concentration in the culture medium is described by the mass conservation equation at steady state condition R net refers to volumetric oxygen consumption, D e is the oxygen diffusion coefficient in culture medium that is considered as a fixed parameter. c refers to oxygen concentration within bioreactor. Since, the thickness of the cellular monolayer is very low, R net is considered as zero, and the effect of cellular uptake is considered as the boundary condition over the bottom wall. Eventually, oxygen concentration is obtained from the solution of convection-diffusion equation Neglecting the variation of oxygen concentration along the depth, the equation is solved in two dimensions. c(x, y) shows the oxygen concentration within the bioreactor. The length of the bioreactor, L, is much larger than its height, L/H >> 1. The diffusive term in the axial direction is included because of low flow rates. This term has been neglected in the previous studies. 18 However, the diffusive term can be neglected for relatively large flow rates (u >> D e ∕L).
The boundary conditions of the bioreactor for solving Equation (5) are as follows. 1. Oxygen is entered at a constant concentration (C in ) at the inlet of the microchannel as 2. Cellular oxygen uptake ( ) follows Michaelis-Menten kinetics that is independent of cell density [19][20][21] V max and K M are uptake parameters that V max is the maximum oxygen uptake rate of cell and K M is the oxygen concentration at which = 0.5 V max . The total uptake flux of oxygen over the bottom wall is assumed to be a sum of the individual uptake rates. Cell density (α), defined as the number of cells per unit area of the channel, is employed to obtain the total oxygen consumption ) .
3. The oxygen flux through the oxygen-permeable membrane at the top of the channel is modeled by constant gas phase partial pressure of oxygen For the bioreactor with oxygen-impermeable membrane, the right side of Equation (9) is set to zero.
The following assumptions are considered in order to simplify the analytical study. 1. Average flow velocity was used for the flow velocity in Equation (5). In the other words, velocity is uniform, constant, and equal to axial average velocity.
2. Since the value of K M = 0.005 mM in comparison with the concentration is negligible, it is logical to suppose c(x, ) K M +c(x, ) = 1. Thus, Equation (9) is reduced to 17 Table 1 shows the dimensionless numbers employed in this study. To determine oxygen tension in the bioreactor with oxygen-permeable membrane, as shown in Figure 1A, the following equations along with the boundary conditions were solved:

Bioreactor with an internal oxygen-permeable membrane
Superposition method was employed to solve the equation. Although convection-diffusion equation is a homogeneous linear differential equation, three boundary conditions are nonhomogeneous. 22 Thus, a change of variable was employed to eliminate heterogeneity in the oxygen-permeable membrane boundary condition It is known that the general solution of C 1 (X, Y) could be written by the superposition of two particular solutions. Thus, the equation could be separated into two parts where C 2 (Y) satisfies the following equations with one heterogeneity in Y-direction: Solving Equations (18)-(20) gives Equation (21) as The other nonhomogeneity of boundary conditions belongs to the channel inlet that yields to the second problem C 3 (X, Y) in the superposition Equation (22) was solved by the method of separation of variables. For the solution of C 3 (X, Y), the eigenfunction obtained as λ n tan(λ n ) = Sh and the related eigenvalues (λ n ) are the roots of this transcendental equation, whereas this equation does not have a precise analytical solution, Maple is used for each Sherwood number. The first 21 eigenvalues for a specific Sherwood number was determined. Finally, by inserting Equation (21) into Equation (23), the dimensionless concentration distribution of oxygen is obtained as

Bioreactor with an oxygen-impermeable membrane
Oxygen flux through the top flat plate (Y = 1) of a bioreactor with no permeable membrane is equal to zero (Sh = 0) ( Figure 1B). Since the result of Equation (26) yields to infinite once the Sherwood number is zero, a new solution for this boundary condition is considered. To this end, superposition is implemented in three different concentrations as C 1 , C 2 , and C 3 : Heterogeneity along the Y-direction for cell surface boundary condition was applied to solve C 2 (Y). By the use of Equation (29), C 3 (X) was obtained as an exponential function C 2 (Y) and C 3 (X) functions are inserted in the nonhomogeneous boundary condition of C 1 (X, Y) Eigenvalues are obtained through λ n = nπ for the solution of C 1 (X, Y). Thus, dimensionless concentration distribution of oxygen within the bioreactor with oxygen-impermeable membrane is given by

Oxygenator
The function of oxygenator is to supply oxygen to its adjacent microchannel. As for the bioreactor with an internal oxygen-permeable membrane, oxygen is directly supplied into the bioreactor through its top wall. Thus, there is an oxygen flux in Y = 1 through the oxygen-permeable membrane with no cellular oxygen uptake. For a bioreactor with no oxygenation, the oxygen flux in Y = 1 is set to be zero.

Operating conditions
The main objective of this study is to provide a mathematical model to determine oxygen tension in microchannel bioreactors with different oxygenation arrangements. To perform the analysis, three dimensionless numbers including Pe, Da, and Sh were employed. Sherwood number represents permeation through the oxygenator membrane to the diffusion rate in the fluid, 8 while Damkohler number indicates the rate of cellular oxygen consumption that is affected by cell density and cell type. Peclet number as a delegate of medium velocity varies between 5 and 250 that is corresponding to velocities between 0.1 and 5.25 mm/s. Cell density is fixed at 10 5 cells/cm 2 in accordance with what is commonly used in monolayer cell culture experiments. 18,23 Since the oxygen diffusion coefficient in the culture medium is D e = 2.1 × 10 −9 m 2 /s 17 and the value of V max for rat hepatocytes is 4.3 × 10 −16 mol/s/cell, 24,25 the dimensionless number Da varies from 0.01 to 1. For instance, the value of Da = 0.13 for oxygen transport corresponds to the rat hepatocytes cultured at a cell density of 1.38 × 10 5 cells/cm 2 . Sherwood number in case of bioreactor with oxygen-impermeable membrane and oxygen-permeable membrane is set to be 0 and 0.27, respectively. 18

Bioreactor design
Hepatic cells and cardiomyocytes were assumed to be cultured in the bioreactors. Among all cases of drug failure, it is estimated that cost of safety liabilities for cardiovascular systems is 45%, while 32% are due to hepatotoxicity. 26 Moreover, almost all drugs interact with the liver, as the largest organ of the body, due to its metabolism and detoxification functions. Therefore, liver toxicity plays a dramatic role in the development of the new drug. 27 Furthermore, most developing drugs are often tested for cardiotoxicity with cardiomyocytes, because the heart is the core organ of the circulatory system. 26 Normal values of oxygen partial pressure (PaO 2 ) in the cellular microenvironment are determined in various human tissues based on organ function and cellular behavior. 14 Oxygen pressure in liver is 5.4 ± 0.7%, 14,28,29 meaning that, if the nondimensionalized oxygen concentration is less than 0.22, the hepatocytes experience hypoxia. Oxygen tension is nondimentionalized by 21% as inlet oxygen. Most liver functions are achieved by hepatocytes, which account for approximately 80% of the total volume of the liver. 1 In addition, approximately 70% of the heart cells are made up of cardiomyocytes, which are considered as a representative of the heart cells. 30,31 The value of V max for human-induced pluripotent stem cells derived cardiomyocytes (hiPSC-CMs) is 1.97 × 10 −16 mol/s/cell. 25 Value of oxygen pressure for arterial blood and venous blood, respectively, are 13.2% and 5.3%. 14 Thus, left atrium and left ventricle are exposed to a blood flow with the nondimensionalized oxygen tension of 0.25, whereas right atrium and right ventricle experience oxygen tension of around 0.62. Figure 2 shows the effect of Pe number (Pe of 5, 25, and 125) on the distribution of oxygen concentration over the cellular monolayer for bioreactors without oxygen-permeable membrane (Sh = 0) and with oxygen-permeable membrane (Sh = 0.27). Once Pe/γ > 1 and Da > 1, diffusive gradient through Y-direction finds important role on oxygen concentration distribution. 19 The average oxygen concentration in each cross section, C avg = ∫ 1 0 C d , is not a suitable representative of oxygen availability for the cells. Thus, the oxygen concentration over the cellular layer, C s = C(x, 0), was considered as proper criterion to evaluate the oxygen availability within the bioreactor. Oxygen concentration is nondimensionalized by C sat = 0.225 mM that is equivalent for 21% oxygen partial pressure in aqueous solutions at 37 • C. 24 In the case of the bioreactor without an oxygen-permeable membrane, only less than 20% of the microchannel length receives sufficient oxygen for the survival of the cells at Pe = 5 ( Figure 2). In addition, at Pe = 25, approximately 15% of cell surface area is oxygen depleted and anoxia occurs. Use of oxygen-permeable membrane increases the oxygen supply to the bioreactor. As for Pe = 5, there is a balance between oxygen supply and oxygen uptake after the second half of the channel length. At high Pe number of 125, the effect of oxygen-permeable membrane on channel oxygenation was minor due to high flow velocity and short resident time. Figure 3A shows the effect of Pe number on outlet oxygen concentration (C out = C(1, 0)). The effect of Da (representation of cell density) on oxygen distribution at Peclet numbers of 5 and 250 for two designs of the bioreactors are shown in Figure 3B. Linear drops for both Peclet numbers are observed as the cell density increases. Figure 3C displays the effects of channel length on the output concentration for three different Peclet numbers. At Pe = 5, the oxygen concentration drop reaches to zero for the bioreactor with no oxygen-permeable membrane, if the channel length is more than 20 mm while the output concentration in bioreactor with membrane remains unchanged because of the balance between the cellular oxygen consumption and oxygen supply.

Oxygen concentration in a two-bioreactor single-pass perfusion circuit
Oxygen tension in a single-pass perfusion circuit with two interconnected bioreactors was determined. The circuit was modeled as two microchannels resembling the bioreactors interconnected through a microchannel (Figure 4). Four cases of fluidic channel were considered.
In Case 1, the whole fluidic channel including the bioreactors and the interconnecting channel is oxygen impermeable. In Case 2, the bioreactors have oxygen-permeable membrane. For Case 3, the two bioreactors have oxygen-impermeable membrane while the interconnecting channel has an oxygen-permeable membrane. For Case 4, all three microchannels have oxygen-permeable membrane. The interconnecting channel has no cultured cells. This channel is implemented within the design of the fluidic circuit to connect the two bioreactors and also create a uniform oxygen tension within the flowing medium before entering the second bioreactor. Owing to the fact that the inlet boundary condition for the analytical solution should be constant concentration, placement of L c is inevitable. For the four cases, it is assumed that the first bioreactor contained hepatocytes and the second bioreactor contained hiPSC-CMs. Two different cell types were selected to investigate the effect of cell type (indicated by V max ) on the oxygen tension. The effect of oxygenator on the distribution of oxygen tension in the single-pass fluidic circuit was investigated ( Figure 5). As for all cases, the increase of Pe number generally enhances the oxygen tension within the fluidic circuit. As for Case 1, no oxygen would reach the second bioreactor at Pe numbers of 5 and 25. In addition, it is quite clear that the use of oxygen-permeable membrane within the fluidic circuit could enhance oxygen supply to the cellular layer.
Moreover, the rate of oxygen tension drop within the first bioreactor is more than the second bioreactor that is mainly due to higher oxygen uptake of the hepatic cells (higher V max ). The findings reveal that the oxygen-permeable membrane has a very distinct role for the oxygenation of the bioreactors at lower Pe number flows.
To provide a better understanding on the role of oxygenator, the effect of Pe number (flow rate) and the cell density on the oxygen tension of three points within the Case 3 was studied ( Figure 6A). At low Peclet number, the diffusion transport is dominant. The oxygen concentration at Point 1 is zero for Peclet numbers of lower than 29; that is indicative of inadequate oxygen supply by the flow, ie, convective transport ( Figure 6B). However, at higher Peclet numbers, the oxygen tension increase due to dominant convective delivery of oxygen. In contrary to Point 1, the zero-oxygen tension in the second bioreactor is limited to the Peclet numbers within the range of 5 to 8. This is mainly associated with the role of oxygenator located before the bioreactor and lower oxygen consumption (V max ). The concentration of Point 2 is controlled by diffusion from the oxygenator across the channel at low Pe numbers. However, for Pe numbers higher than 29, oxygen supply to Point 3 mainly occurs by convective transport. Figure 6C shows the effect of cell density on oxygen tension in the three points. As the cell density within the bioreactors increases, the oxygen tension at all points decreases. The findings could provide critical insights for the design and operation of microfluidic bioreactors made from oxygen-impermeable materials. Such bioreactors require oxygenators for sufficient oxygen supply.
Oxygenator has a critical role on oxygen supply to the second bioreactor. Modulation of the oxygenator length is of critical important to regulate the oxygen tension in the bioreactor. To this end, an analysis was performed to determine the

Oxygen concentration in a two-bioreactor closed-loop perfusion circuit
An analysis was performed to determine the oxygen tension distribution in two bioreactors interconnected in a closed-loop recirculation fluidic circuit. Such fluidic configuration is mainly used for OoC applications once the cell-secreted biomarkers and drug metabolites are desired to be maintained and accumulated in the flowing culture medium. To this end, a fluidic circuit with two bioreactors and two oxygenators were considered ( Figure 8A). The oxygenators should maintain a constant oxygen tension for both bioreactors over time. To perform the analysis, the oxygen tension at the outlet of the second bioreactor was considered as 0.3 that is equal to the oxygen tension at the entrance of the first oxygenator. The criterion to determine the length of the first oxygenator is to maintain an oxygen tension of 0.22 at the outlet of the first bioreactor. The length of the second oxygenator was determined to produce oxygen tension of 0.3 at the outlet of the second bioreactors. Figure 8B shows the optimum length for the oxygenators.

CONCLUSIONS
This article presented an analytical solution to determine oxygen tension in microchannel bioreactors with and without oxygenators. Axial diffusion through the channel and the diffusion along the channel height were considered. The results show that the oxygen tension within a bioreactor is influenced by the flow rate of the culture medium, cell type and cell density, and the bioreactor dimensions. In a fixed flow rate, the amount of oxygen supply could be modulated through changing the oxygenator length. From the design point of view, in multicell fluidic systems, the difference of oxygen uptake rate for each cell type should be kept into considerations for proper determination of the oxygenator length.