Fluid Dynamic Modeling to Support the Development of Flow-Based Hepatocyte Culture Systems for Metabolism Studies

Introduction

There has been a coordinated, international call to increase the use of in vitro and in silico methods for hazard identification and risk assessment (Dix et al., 2007; Williams et al., 2009; Adeleye et al., 2015; De Boer et al., 2015). In response, scientists have developed advanced three-dimensional (3D) cell cultures with more physiologically relevant microstructures and environments that better match those in the intact organism (Marga et al., 2012; Kostadinova et al., 2013; Shulman and Nahmias, 2013; Ballard et al., 2015). Compared with traditional two-dimensional (2D) cultures, 3D cultures hold the potential to generate biological responses from chemical exposures that are more predictive of human in vivo responses (Zare-Mehrjardi et al., 2011; Pineda et al., 2013).

The liver plays an important role in the systemic clearance of xenobiotics. Currently, freshly isolated primary hepatocytes are regarded as the most relevant experimental system to study hepatic metabolism, as well as metabolism-mediated effects of chemicals and pharmaceuticals (Lecluyse and Alexandre, 2010). Freshly isolated primary hepatocytes express most of the proteins found in the human liver, including those involved in metabolism, membrane transport, and receptor-mediated processes (Li et al., 1995; Vildhede et al., 2015; Yang et al., 2015). However, one of the major drawbacks with primary hepatocytes is the rapid change in phenotype observed in culture. Several techniques have recently been developed to stabilize hepatic phenotype over longer periods, which may enhance the predictive power of primary cell culture for in vivo responses (Vinci et al., 2010; Ballard et al., 2015). These advanced cell culture models have shown their ability to better mimic physiological microenvironments in the human liver in healthy and diseased states (Chan et al., 2013; Ballard et al., 2015). Recently, a dynamic culture device with a 3D scaffold for cell attachment [RealBio^® (RB)] was tested for long-term primary hepatocyte culture (Choi et al., 2014). In this method, cultured primary hepatocyte suspensions are seeded into a scaffold where they form 3D structures. The cultures remained viable for over 7 weeks and maintained some of their liver-specific functionality, e.g., sustained albumin production (Choi et al., 2014). Other similar advances have also been observed with several other culture methods (Dash et al., 2013; Ukairo et al., 2013; Ballard et al., 2015).

There is a growing appreciation that the use of these advanced culture systems can be beneficial for a variety of endpoints within many fields of research (Breslin and O’Driscoll, 2016; Clevers, 2016; Leek et al., 2016). In toxicology, bioreactors have the potential to improve chemical and drug safety assessment as they are expected to better recapitulate biological responses to chemical exposures. The combination of exposure control and ease of sampling with more organotypic cell cultures also make bioreactors attractive for kinetics and metabolism investigations. The potential ability to reproduce exposure profiles of parent chemical and metabolites combinations within culture systems can also improve concentration-effect responses and in vitro to in vivo extrapolations.

Many of the bioreactor systems are flow-based culture systems, and computational fluid dynamics (CFD) modeling is well suited to optimize bioreactor design as it provides a tool to describe fluid dynamics in the engineered culture devices (Kang et al., 2013; Hsu et al., 2014). CFD modeling allows parameter values to be varied and the effects on the system to be characterized rapidly. This includes effects that are easier to describe mathematically than to measure experimentally. For instance, CFD modeling can generate detailed descriptions of fluid flow and nutrient transport within bioreactor systems. It can also be used to describe spatiotemporal distributions of a test compound or oxygen within the system. CFD and other computational methods have previously been applied to enhance our understanding of system parameters fundamental for cell responses within complex 3D culture environments (Kang et al., 2013; Hsu et al., 2014).

In this study, CFD modeling was used to assist the development of a flow-based hepatocyte culture system for long-term hepatocyte culture suitable for extended investigations of metabolic clearance or metabolite profile identification of low clearance compounds. Previous data indicate that the scaffold-based RB system supports long-term hepatocyte viability but not maintained metabolic capabilities. We hypothesize that increased cell–cell interactions fostered in other types of 3D culture conditions in combination with dynamic flow culture could improve the phenotypic stability. Two other bioreactors, QuasiVivo^® (QV) and fluidized bed (FB), were tested for this reason. They do not require a specialized scaffold but are compatible with some existing 3D methods for long-term hepatocyte culture, including alginate-encapsulated hepatocytes. CFD models were developed for the RB, QV, and FB bioreactors with the aim to identify a bioreactor suitable for alginate bead culture.

Materials and Methods

Materials

William’s E medium was purchased from Life Technologies (Grand Island, NY, USA). Trypan blue (TB) and 7-ethoxycoumarin were purchased from Sigma-Aldrich (St. Louis, MO, USA). RB bioreactors were purchased from RealBio^® Technology, Inc. (Kalamazoo, MI, USA). Quasi-Vivo^® bioreactors were provided by Kirkstall Ltd. (Rotherham, UK). The FB bioreactor was manufactured by Professor Jeffrey M. Macdonald and Vaibhav Hans, University of North Carolina (Chapel Hill, NC, USA).

Bioreactor Design

Three bioreactor systems (RB, FB, and QV) were selected for inclusion in the current study. Diagrams of these systems are shown in Figure 1.

FIGURE 1

Figure 1. 2D (A,C,E) and 3D (B,D,F) geometry models of the bioreactors used in the simulations. 2D models are from COMSOL; 3D models are from Solidworks. The three bioreactor types are RealBio (A,B), fluidized bed (C,D), and QuasiVivo (E,F).

The RB bioreactor (Figures 1A,B) is a scaffold-based system where hepatocyte suspensions are seeded and cultured on a 3D woven polycarbonate scaffold (Pfund et al., 2013). The 1-mm thick scaffold is maintained between two media compartments with dynamic flow. Both media compartments are in equilibrium with gas compartments through gas-permeable membranes (Figures 1A,B). All flows were modeled as concurrent, and all media inlets and outlets were of 1-mm internal diameter (ID). The total volume of the RB bioreactor is 7.0 mL.

The FB bioreactor (Figures 1C,D) was custom developed specifically for use with alginate-encapsulated primary hepatocytes. This bioreactor was machined from polysulfone and consists of three main parts: top and bottom pieces creating two chambers for media and a spacer piece wherein the cell compartment is located (Figures 1C,D). The FB bioreactor has a total volume of 4.8 mL. The cell compartment is separated from the media compartments with two polypropylene filters with 100-μm pore size. The media inlet was 1.5-mm ID, and the outlet was 3-mm ID.

The QV bioreactor (Figures 1E,F) has a single compartment with a single media inlet (1-mm ID) and outlet (2-mm ID) (Gordon et al., 2014). The volume of the chamber is 4.0 mL, and it can accommodate 2D or 3D cell cultures. 2D cultures are supported on a slide that rests on six small projections from the floor of the chamber, while alginate-based 3D cultures can be added to the chamber without the slide. The base of the bioreactor is made of Altuglas SG7 acrylic, while the lids are a styrene-based thermoplastic elastomer (TPE). The QV chamber has a polycylindrical shape with a slanted “roof” and offset outlet meant to promote mixing and reduce shear stress along the chamber bottom where the cells are located.

Simulation Software

Bioreactor geometries were described using finite element modeling (FEM) algorithms as implemented in COMSOL Multiphysics 4.3 (The COMSOL Group, Stockholm, Sweden) and Solidworks 2013 (Dassault Systems, Waltham, MA, USA) (Figure 1). Most simulations were performed using COMSOL’s stationary solver that finds time-invariant solutions. Simulations of test compound distribution were performed using COMSOL’s time-dependent solver up to t = 60 min. The remaining simulations were performed using Solidworks.

Fluid Dynamics Modeling

Flow velocity, turbulence, and shear stress profiles in each bioreactor were investigated using 2D and 3D CFD modeling. Fluid motion within the bioreactors was described by a Navier–Stokes equation assuming incompressible, isothermal Newtonian fluids. Equation 1 describes the momentum balance, while Eq. 2 describes continuity within incompressible fluids.

$\mathrm{\rho }\frac{\mathit{\text{∂u}}}{\mathit{\text{∂t}}}-\mathrm{\mu }{\nabla }^{2}u+\mathrm{\rho }\left(u·\nabla \right)u+\nabla P=F$ (1)

$\nabla ·u=0$ (2)

where ρ is density (kilograms per cubic meter), μ is viscosity (kilograms per meter per second), P is pressure (kilograms per meter per square second), F is volume force (kilograms per square meter per square second), u is the 3D velocity field (meter per second), and ▽ is del (aka nabla) operator, $\nabla ={\overrightarrow{e}}_x\frac{\partial }{\mathit{\text{∂x}}}+{\overrightarrow{e}}_y\frac{\partial }{\mathit{\text{∂y}}}+{\overrightarrow{e}}_z\frac{\partial }{\mathit{\text{∂z}}}$ for Cartesian coordinates (Buchwald, 2011). Model parameters used were viscosity = 10⁻³ kg/m/s, fluid density = 1000 kg/m³, temperature = 37°C, and no-slip boundary conditions. The volumetric flow rate was 0.5 mL/min (QV) or 1 mL/min (RB and FB). Outlets were modeled with a pressure of 0 kg/m/s² with no viscous stress or environmental pressure.

Shear stress experienced by the cells was described with Eq. 3. The properties of cell culture media were assumed to be identical to those of water, the reference fluid.

${\mathrm{\tau }}_{\text{wall}}=\mathrm{\mu }{\left(\frac{\mathit{\text{∂u}}}{\mathit{\text{∂y}}}\right)}_{y=0},{\mathrm{\tau }}_{\text{wall}}={\mathrm{\tau }}_{\text{cells}}$ (3)

where τ is shear stress (kilograms per meter per square second), μ is viscosity (kilograms per meter per second), u is flow velocity parallel to the wall (meter per second), and y is the distance to the wall (meters).

Bioreactor fluid dynamic models were developed for 2D and 3D geometries. The pore size of the RB bioreactor scaffold was modeled as 100 μm. For the FB bioreactor, a porosity factor of 0.3 was used to simulate alginate beads in the cell compartment. Sensitivity analysis showed no significant differences in cell compartment flow for porosity values from 0.1 to 0.7. The filters on either side of the FB cell compartment were modeled as a 100-μm thick layer with evenly distributed pores of 100-μm diameter. To investigate the impact of the choice of spatial dimensions for the FB bioreactor, three variations of the width of the FB parabola were modeled: 27.3, 22.2, and 16.7 mm corresponding to the original (Figure 2A), narrow (Figure 2B), and oval (Figure 2C) design, respectively.

FIGURE 2

Figure 2. Three variations on the fluidized bed (FB) bioreactor with different spatial dimensions using Solidworks. (A–C) Flow particle study results and (D–F) flow trajectory study results shown from the top. Option 1: 27.3 mm maximum width (A,D); Option 2: 22.2 mm maximum width (B,E); and Option 3: 16.7 mm maximum width (C,F).

Simulation of Small Molecule and Oxygen Distributions

The CFD models described in Section “Fluid Dynamics Modeling” were used to simulate the spatiotemporal distribution of oxygen and a representative small molecule in each bioreactor. Parameters describing these simulations are listed in Table 1.

TABLE 1

Table 1. Physical constants used in the CFD models to describe the mass transport of oxygen and a small molecule test compound in each bioreactor.

Test Compound Distribution

In our simulations, the test compound enters the bioreactor through the inlet at a fixed concentration, as described in Section “Experimental Studies for Small Molecule Distribution.” Mass transport was by advection and diffusion, as described below (Eq. 4):

$\frac{\mathit{\text{∂c}}}{\mathit{\text{∂t}}}+\nabla ·\left(-D\nabla c\right)=R-u·\nabla c$ (4)

where c is the 3D chemical concentration (moles per cubic meter), D is the chemical-specific diffusion coefficient (square meter per second), R is the reaction rate (moles per cubic meter per second), and u is 3D velocity field (meter per second) (Buchwald, 2011). For the purpose of this study, metabolism was not included and R = 0 mol/m³/s for the test compound was used. The diffusion coefficient for a small molecule was set at 1 × 10⁻⁸ m²/s (Brody and Yager, 1997; Macdonald et al., 1999). A sensitivity analysis using a range of diffusion coefficients from 1 × 10⁻⁶ to 1 × 10⁻¹⁰ m²/s showed minimal influence of this parameter on the time it took to establish a uniform distribution of the test compound in the bioreactors.

Oxygen Distribution

Mass transport of oxygen was handled similarly to that of the test compound. Oxygen diffusion was described using Eq. 4 but using a diffusion coefficient of 3 × 10⁻⁹ m²/s (Buchwald, 2011). The concentration of oxygen in the media at t = 0 was assumed to be in equilibrium with the concentration of oxygen in air, as described by Henry’s law (Eq. 5).

$P_{{\text{O}}_2}=K_{{\text{O}}_2}·C_0$ (5)

where $K_{{\text{O}}_2}$ is the Henry’s law constant and $P_{{\text{O}}_2}$ is the partial pressure of oxygen in air. See Table 1 for a list of parameter values and references. The concentration was also used as the concentration of oxygen in the media entering the bioreactor, and the concentration of oxygen in the media at the air–liquid interface membranes in the RB and FB bioreactors.

In the presence of cells, local oxygen concentrations will depend on mass transport of oxygen in the bioreactor and the rate of oxygen consumption by the cells. The rate of oxygen consumption by the cells was assumed to follow Michaelis–Menten (MM) kinetics (Eq. 6). A step-down function was used to simulate the reduction in the rate of oxygen consumption as the local concentration of oxygen approaches the critical concentration for cell survival (2.8 μM) as previously described (De Groot et al., 1988). This approach has been shown to improve cellular oxygen consumption predictions at low oxygen concentrations.

$R=V_{\text{max}}·\frac{c}{c+k_{\text{m}}}·\mathrm{\delta }\left(c>c_{\text{cr}}\right)$ (6)

$\mathrm{\delta }\left(c>c_{\text{cr}}\right)=\left\{\begin{array}{cc}0\hfill & \text{for\hspace{0.17em}}c<-c_{\mathit{\text{cr}}}\hfill \\ 0.5+0.75\left(\frac{c}{c_{\text{cr}}}\right)-0.25{\left(\frac{c}{c_{\text{cr}}}\right)}^3\hfill & \text{for}-c_{\mathit{\text{cr}}}\le c\le c_{\mathit{\text{cr}}}\hfill \\ 1\hfill & \text{for\hspace{0.17em}}c>c_{\mathit{\text{cr}}}\hfill \end{array}\right.$ (7)

where V_max is the maximal oxygen consumption rate (moles per cubic meter per second), k_m is the MM constant (moles per cubic meter), c_cr is the critical oxygen concentration (moles per cubic meter), and δ is the step-down function used to simulate the decline in oxygen consumption by the hepatocytes when local oxygen concentrations approach c_cr (Eq. 7). See Table 1 for a list of parameter values and references. A cell-normalized oxygen consumption rate of 0.4 nmol/s/10⁶ cells was used, along with a k_m value of 5.6 mmHg (Rotem et al., 1992; Balis et al., 1999). The total number of cells in the RB, FB, and QV bioreactors was estimated at 15, 70, and 1.5 million cells, respectively. V_max values were calculated by multiplying the cell-normalized oxygen consumption rate by the number of cells in the bioreactor and then dividing by the volume of the cell compartment. Cell numbers were estimated based on seeding density (RB) or the volume of alginate beads added to the bioreactor (FB and QV). Three bead diameters were considered: 250, 500, and 1000 μm. The alginate beads were assumed to be 50% (v/v) hepatocytes, and a hepatocyte diameter of 20 μm (Kuntz and Kuntz, 2008) was used to convert to total cell number. The FEM algorithm as applied in COMSOL was used to describe local oxygen concentrations within the beads. All oxygen transport simulations were run using the stationary solver.

Experimental Studies for Small Molecule Distribution

Two small molecules were used as test compounds, 7-ethoxycoumarin (7EC) and TB. Each was tested separately. The compound was added to the cell-free system using matched flow and other incubation conditions used for fluid dynamic modeling described above. The media was continuously recirculated through the system with a peristaltic pump. William’s E medium with 10 μM 7EC or 0.04% TB was flowed through each bioreactor in the absence of cells. The media leaving the reactor was sampled up to 60 min (RB) or 30 min (FB and QV) to determine compound distribution kinetics. Compound sample concentrations were determined by fluorescence (7EC, excitation at λ = 323 nm, emission at λ = 386 nm) or absorbance (TB, λ = 580 nm). Results were normalized to the maximum concentration that was measured in the time course.

Results

Optimizing Chamber Design in the FB Bioreactor

As the FB bioreactor was developed in-house, we were able to alter chamber geometry as needed. Three alternatives based on a shape with a parabolic cross-section were considered. The “original,” “narrow,” and “oval” designs had maximum widths of 27.3, 22.2, and 16.7 mm, respectively. The models were otherwise identical and were investigated for the impact of these design changes to flow pattern and turbulence. An uneven flow pattern with higher cross-flow through the cell compartment near the bioreactor outlet was observed in the original FB chamber (Figure 2A). With the narrow parabola design, the cross-flow pattern shifted to be stronger nearer the inlet (Figure 2B). The uneven flow pattern near the cell compartment was observed again near the outlet for the oval design (Figure 2C). No turbulence was observed in any of the FB chamber designs (Figures 2D–F).

Bioreactor Flow Velocity, Turbulence, and Shear Stress

Computational fluid dynamics modeling simulated flow velocities, turbulence, and shear stress patterns within the RB, FB, and QV bioreactors. Both 2D and 3D models of velocity were examined. Turbulence and shear stress were analyzed using 3D models.

RB Bioreactor

Two-dimensional simulations of the RB bioreactor showed peak flow velocities of 2.9 × 10⁻³ m/s and 6.2 × 10⁻⁴ m/s at the chamber inlet and scaffold level, respectively (Figure 3A). 3D simulations showed peak flow velocity at the chamber inlet to be 0.0137 m/s, approximately five times higher than the 2D predictions. 3D simulations also showed velocities at the scaffold ranging from 1 to 8 × 10⁻⁵ m/s, depending on the position (Figure 4A). Flow was predicted to be laminar throughout the bioreactor, with turbulence possible only near the inlet and outlet where there are large constrictions or expansions (Figure 4D). Media cross-flow through the scaffold was observed between the two compartments, with predicted scaffold shear stress levels between 34.6 and 484 μPa (Figure 4G).

FIGURE 3

Figure 3. 2D simulation of flow velocity in the bioreactors using COMSOL: RealBio (A), fluidized bed (B), and QuasiVivo (C). The volumetric flow rate was 1 mL/min (RealBio and fluidized bed) or 0.5 mL/min (QuasiVivo).

FIGURE 4

Figure 4. 3D simulation of flow velocity and shear stress in the bioreactors. (A–C) Velocity profiles of bioreactors shown in the vicinity of the cell compartment/scaffold (COMSOL). (D–F) Flow trajectories, to identify turbulence within the chambers (Solidworks). (G–I) Shear stress profiles along the cell compartment/scaffold surface (Solidworks). Three bioreactors are shown: RealBio (A,D,G), fluidized bed (B,E,H), and QuasiVivo (C,F,I).

FB Bioreactor

Two-dimensional simulations of the FB bioreactor showed a peak inlet flow velocity of 8.02 × 10⁻³ m/s and cell compartment velocities ranging from 1 to 2 × 10⁻³ m/s (Figure 3B). 3D simulations showed a peak inlet velocity of 5.76 × 10⁻³ m/s and cell compartment flow velocities of 1–9 × 10⁻⁴ m/s (Figure 4B). The velocity in the cell compartment near the inlet was 9.8 × 10⁻⁴ m/s. No turbulence was observed in the FB bioreactor (Figure 4E). Cell compartment shear stress ranged from 1070 to 4260 μPa (Figure 4H).

QV Bioreactor

Two-dimensional simulations of the QV bioreactor showed a peak inlet flow velocity of 4.45 × 10⁻³ m/s and velocities near the alginate beads at the bottom of the chamber of 2 × 10⁻⁷ m/s (Figure 3C). 3D simulations showed a peak inlet velocity of 0.022 m/s, almost five times higher than the 2D predictions. The velocity near the alginate beads varied by position in the range of 3–15 × 10⁻⁶ m/s (Figure 4C). Some turbulence was observed in the QV bioreactor in the downward flow paths after the inlet, as well as at the periphery of the bottom of the chamber (Figure 4F). QV shear stress near the alginate beads varied by position in the range of 4–57 μPa (Figure 4I). These results are in good agreement with those previously reported for a variation of this bioreactor from the same manufacturer (Mazzei et al., 2010).

Test Compound Mass Transport in the Bioreactors

CFD Modeling Results

Two and three-dimensional CFD modeling was used to simulate mass transport of a representative small molecule in the three bioreactors. The 2D model predicted uniform compound distributions at 18, 6, and 6 min in the RB, FB, and QV bioreactors, respectively. Figures 5A–C show the concentration of a small molecule after 6 min. The corresponding 3D models predict it would take 40, 8, and 17.5 min for a small molecule to distribute uniformly in the RB, FB, and QV bioreactors, respectively. Figures 5D–F depict the concentration of a small molecule after 10 min. When compared, the FB bioreactor models resulted in similar predications using either 2D or 3D models. However, for the RB and QV bioreactors, the predicted time to equilibration was two to three times higher in the 3D predictions compared to the 2D.

FIGURE 5

Figure 5. Simulation of compound distribution in the bioreactor using Solidworks. Chemical species transport simulated with 2D models at 6 min (A–C) or 3D models at 10 min (D–F).

Experimental Results

7-ethoxycoumarin and trypan blue were added to the media at the bioreactor inlet, and the media was sampled at the bioreactor outlet. The concentration at the outlet was normalized to the maximum concentration measured. The results of this experiment show that it took 25 min (RB), 8 min (FB), and 10 min (QV) for the outlet concentration to reach equilibrium (Figures 6A–C). In the QV bioreactor, the outlet concentration rose rapidly to a maximum and then began falling again (Figure 6C). This result is consistent with equilibration in a recirculating system where the flow initially carries the chemical through the bioreactor in a current along the top of the chamber, followed by slower back-mixing inside the chamber causing a drop in the outlet concentration over time. In this scenario, using the initial peak would underestimate the time to system equilibration.

FIGURE 6

Figure 6. Experimentally measured concentrations of 7-ethoxycoumarin (blue diamonds) or trypan blue (red squares) exiting the bioreactor, expressed as a percentage of the maximum detected concentration. Results for RealBio (A), fluidized bed (B), and QuasiVivo (C).

For the FB bioreactor, both 2D and 3D CFD predictions were found to be in agreement with the experimental data. For the RB bioreactor, 2D predictions resulted in a 28% underestimation of the time to equilibration, while the corresponding 3D prediction resulted in a 60% overestimation. A similar pattern was observed for the QV bioreactor, and time to equilibration was underestimated by 40% by the 2D predictions and overestimated by 75% by the 3D predictions. In general, the 2D models averaged at a 1.6-fold underestimation of the time to equilibration, while the 3D models resulted in an average of 1.4-fold overestimation of the time to equilibration.

Oxygen Distribution in the Bioreactors

Distribution of oxygen within the RB, FB, and QV bioreactors was simulated using 2D and 3D models. For the RB bioreactor, the 2D model predicted oxygen levels of 160 μM, well above the critical oxygen concentration, for the majority of the scaffold over the first hour in culture (Figure 7A). Using the 3D CFD model, oxygen concentrations down to 3.7 μM, less than twice as high as the critical oxygen concentration of 2.8 μM, were predicted closer to the chamber outlet and for parts of the center of the scaffold (Figure 8A).

FIGURE 7

Figure 7. 2D simulation of oxygen concentration in the bioreactors using COMSOL. Results for RealBio (A), fluidized bed (B), and QuasiVivo (C).

FIGURE 8

Figure 8. 3D simulation of oxygen concentration in the bioreactors using COMSOL. Results for RealBio (A), fluidized bed (B), and QuasiVivo (C).

In the FB bioreactor, 2D simulations indicate an average oxygen concentration of 130 μM in the area of the cell compartment after 60 min, with a moderate gradient ranging from 220 μM near the inlet to 100 μM at the outlet (Figure 7B). A low-oxygen area immediately above the inlet had a local velocity and oxygen level of 10⁻⁹ m/s and 1.9 μM, respectively (Figures 3B and 7B). However, the 3D simulation predicted a steady-state oxygen concentration of 4 μM throughout a majority of the cell compartment (Figure 8B). A higher oxygen of 12 μM was predicted around the inlet.

For QV, the 2D simulation results in oxygen levels above 200 μM for a majority of the cell compartment. An area with oxygen concentrations down to 7.5 μM in the area below the inlet was observed after 60 min (Figure 7C). The corresponding 3D simulation resulted in local oxygen concentrations as low as 1.4 μM at the cell level of the bioreactor (Figure 8C).

Optimization of Flow Rate and Oxygen Partial Pressure in the QV Bioreactor

Simulations of the QV bioreactor showed favorable results for compound distribution and shear stress levels close to the alginate beads. However, predicted oxygen levels were too low to support the alginate-encapsulated 3D hepatocyte culture. Different media flow rates and external oxygen concentrations were therefore modeled to optimize the QV system for applications involving alginate-encapsulated 3D hepatocyte culture. As presented in Section “Oxygen Distribution in the Bioreactors,” the local oxygen concentration at the bottom of the QV bioreactor was predicted to be 1.4 μM at a flow rate of 0.5 mL/min. This local oxygen concentration increased to 1.7 and 33 μM at flow rates of 1.5 and 3 mL/min, respectively (Figure 9). Variations in local oxygen concentrations at the bottom of the QV bioreactor were investigated at the atmospheric level of 21% oxygen, as well as at 35 and 50% oxygen. When modeled using a 1.5 mL/min initial flow rate, the local oxygen concentration increased from 1.7 μM to 43 and 180 μM at 35 and 50%, respectively (Figure 10). The results suggest that increased flow rate and higher external oxygen concentration would be beneficial to increase QV applicability for alginate-encapsulated 3D culture of primary hepatocytes.

FIGURE 9

Figure 9. 3D simulation of oxygen concentration in the QuasiVivo for different flow rates using COMSOL. Results for 0.5 mL/min (A), 1.5 mL/min (B), and 3 mL/min (C).

FIGURE 10

Figure 10. 3D simulation of oxygen concentration in the QuasiVivo for different oxygen percentages in the gas phase at a flow rate of 1.5 mL/min (COMSOL). Results for 21% O₂ (A), 35% O₂ (B), and 50% O₂ (C).

Oxygen Distribution within Alginate Beads

Local oxygen concentrations within alginate and hepatocyte beads were investigated for 250, 500, and 1000 μm diameter beads using 2D and 3D simulations. The initial media concentration was set to 214 μM, as calculated using Eq. 5.

The 2D simulations show steady-state oxygen concentrations at the center of the beads of 175, 25, and 4 μM for bead diameters of 250, 500, and 1000 μm, respectively (Figure 11, top row). All oxygen concentrations were predicted to be above the critical oxygen concentration of 2.8 μM, although by less than a factor of 2 for the 1000 μm beads. Corresponding 3D simulations resulted in steady-state oxygen concentration at the center of the beads of 200, 150, and 6 μM for bead diameters of 250, 500, and 1000 μm, respectively (Figure 11, bottom row).

FIGURE 11

Figure 11. Simulation of oxygen concentration inside beads of alginate-encapsulated hepatocytes using COMSOL. 2D model results (top row) and 3D model results (bottom row) for the RealBio (left column), fluidized bed (center column), and QuasiVivo (right column) bioreactors.

Discussion

Over the past decade, there has been an increase in the use of flow-based, long-term culture systems (i.e., “bioreactors”) for many purposes, including the investigation of organ function and cell metabolism. Several hepatic bioreactor systems have been developed with the aim to recapitulate parts of the complex in vivo microenvironment to better maintain the in vivo hepatocyte phenotype over time in culture. Typical parts of this microenvironment include extracellular factors, such as cell–cell interactions, cell–matrix interactions, chemical microenvironment, or physical parameters, such as shear stress (Kuntz and Kuntz, 2008; Lecluyse et al., 2012). Examples of advanced models include scaffold-based 3D cultures, spheroid cultures, or alginate-supported 3D culture. Several of these culture methods maintain hepatocyte morphology and viability over many weeks and previous results indicate that the inclusion of dynamic media flow extends the primary hepatocyte viability even further (Dash et al., 2013; Ukairo et al., 2013; Choi et al., 2014). With the expansion of the use of flow-based systems, CFD modeling has been an important tool to describe internal bioreactor parameters. CFD modeling enables rapid assessment of bioreactor parameters such as internal flow velocity, shear stress, and oxygen concentration that are difficult to determine experimentally. In the present study, CFD modeling was used in three different liver bioreactors as a long-term culture system for primary hepatocytes on 3D scaffolds or in alginate-encapsulated 3D beads. Bioreactor and operational parameters were optimized for use in chemical and drug safety assessments, more specifically for studying xenobiotic metabolism. It should be noted that all of the bioreactors investigated here were designed and optimized by their inventors for purposes other than long-term hepatocyte culture. Our goal was to identify design criteria for optimizing the bioreactors for assessing metabolism and to use computational modeling to explore the advantages and disadvantages of each system when used for our purposes.

An initial exploration was conducted to examine the effect of altering the spatial dimensions of the FB bioreactor to see if improvements could be made to the existing system. Examination of the results from the original system (Figures 2A,D) versus the two alternative systems (Figures 2B,E,F) showed no major improvements. Because no advantage could be observed, the original design was carried forward for future modeling studies.

Two- and three-dimensional models describing the geometry of the three bioreactors were built using two software platforms, COMSOL and Solidworks (Figure 1). Velocity fields were generated using these models to better understand how fluid moves through each system. The RB and FB bioreactors both had simple velocity fields, with fluid moving from inlet to outlet past and through the hepatocyte layer (Figures 4D–E). The RB bioreactor contains a scaffold where hepatocytes are seeded directly, flanked on both sides by two media compartments that are in communication with reservoirs of air via gas-permeable membranes to facilitate greater media oxygenation. Each media compartment possesses its own inlet and outlet, but they are angled so as to encourage a small amount of cross-flow through the scaffold. The FB bioreactor has one inlet and one outlet and is seeded with alginate supported 3D cell cultures in the cell compartment situated between the inlet and outlet. Media entering the bioreactors must flow through the alginate beads before exiting from the opposite side of the bioreactor. In both cases, the flow is almost entirely unidirectional.

In this study, the QV bioreactor is used for the culture of alginate-encapsulated 3D hepatocytes, just as the FB bioreactor. The outlet is offset from the inlet resulting in the media flow tangential to the cell compartment at the bottom of the chamber. A substantial portion of the flow establishes a circulatory pattern with subsequent mixing in the lower half of the chamber (Figure 4F). Unlike the RB or FB bioreactors, the QV system had areas with negligible flow (Figures 4D–F).

Flow regime is another factor to be considered during the optimization process. Turbulent flow can affect local shear stress as well as introduce a chaotic element into the modeling process. No turbulence was predicted for the FB bioreactor and was minimal for the RB and QV bioreactors at the volumetric flow rates under consideration. As such, it was not a major consideration when choosing between the three systems.

Shear stress experienced in the vicinity of the cells was also modeled (Figures 4G–I). A low level of shear stress provides a better representation of the in vivo microenvironment and may improve hepatocyte viability and phenotype stability (Dash et al., 2013; Choi et al., 2014). Conversely, too much shear stress will introduce cell stress that will damage or even kill the cells. The QV bioreactor had the lowest shear stress, followed by the RB bioreactor, with the FB bioreactor having the highest levels.

Accurate kinetic modeling requires a high degree of certainty regarding both the magnitude and the timing of exposure. In the context of this study, it means that a test compound added to the system should, in the ideal case, rapidly establish a uniform concentration throughout the system. Through CFD modeling (Figures 5D–F) and experiment (Figure 6), it was determined that the FB and QV bioreactors both achieved this goal within the first 10 min after dosing. The RB bioreactor took more than twice as long to establish a uniform concentration.

In order to support cellular respiration, sufficient oxygen concentrations are required. CFD modeling showed that both the RB and FB bioreactors had a minimal oxygen concentration above the critical level throughout (Figures 8A,B). The lowest predicted oxygen concentrations for the RB were less than twice the critical oxygen concentration and were observed near the outlet of the chamber and along the center of the scaffold were the media flow was less. In our experiments with the RB system, a lower density of cells was observed in the central region and near the outlet (Figure S1 in Supplementary Material). These areas were predicted to have local oxygen concentrations close to the critical oxygen concentration and almost a complete lack of shear stress, which could have impeded hepatocyte viability. In the QV bioreactor, the oxygen concentration dropped below this level, to 1.4 μM (Figure 8C). However, by increasing the flow rate to 3 mL/min (Figure 9C), or by increasing the flow rate to 1.5 mL/min and supplementing additional oxygen gas to the system above 35% (Figures 10B,C), the minimum oxygen concentration was brought above the critical level. Even in the presence of sufficient oxygen in the surrounding media, diffusion and consumption of oxygen within the alginate beads themselves could lead to insufficient oxygen levels near their centers. Modeling of the beads (Figure 11) showed that a diameter ranging from 250 to 500 μm is preferable. Cell viability of alginate supported 3D cultures of rat hepatocytes was maintained through day 31 of culture in the lab (Figure S2 in Supplementary Material).

The main design parameters under consideration were (1) minimizing shear stress on the cells/alginate beads, (2) minimizing the time to establish a uniform distribution of test compound in the bioreactor, and (3) maintaining oxygen concentrations above the critical level. Taken together, the QV system was selected because it showed the lowest shear stress of any system under consideration, a short time scale for test compound distribution, and sufficient oxygen availability at higher flow rates. The low shear stress was attributable to the placement of the cells at the bottom of the bioreactor, away from the high velocity current at the top of the chamber. The rapid distribution of test compound was due to the flow pattern in the chamber, which promoted mixing down to the bottom without high velocities that could cause shear stress on the alginate beads.

In conclusion, CFD modeling was successfully used to analyze the applicability of three different bioreactors for extended hepatocyte culture, with the potential for improved phenotypic stability. The results show that CFD modeling correctly describes and improves rapid assessments and optimizations of bioreactor design as well as the identification of suitable systems for specific applications. In vitro long-term chemical metabolism and kinetic investigations in the QV bioreactor are ongoing.

Author Contributions

JP designed experiments, acquired the in-lab experimental data, interpreted the simulation results, and wrote the manuscript. Y-SS ran the computational models and assisted with the in-lab experiments. VH designed the fluidized bed bioreactor. MP, JM, GW, MA, HC, and MY contributed to study design, interpretation of the simulation data, and revision of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank Dr. Kristina Wolf for providing us with freshly isolated rat primary hepatocytes for our study, Dr. Andrey Tikunov for preparing the alginate beads, and Mr. David Billings for reviewing the manuscript. We would like to thank RealBio^® Technology, Inc. and Kirkstall Ltd. for providing their products for evaluation in the current study.

Funding

This research was funded by the Long-Range Research Initiative of the American Chemistry Council.

Supplementary Material

The Supplementary Material for this article can be found online at https://www.frontiersin.org/article/10.3389/fbioe.2016.00072

References

Adeleye, Y., Andersen, M., Clewell, R., Davies, M., Dent, M., Edwards, S., et al. (2015). Implementing toxicity testing in the 21st century (TT21C): making safety decisions using toxicity pathways, and progress in a prototype risk assessment. Toxicology 332, 102–111. doi:10.1016/j.tox.2014.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, J. W., and Bhatia, S. N. (2003). Formation of steady-state oxygen gradients in vitro: application to liver zonation. Biotechnol. Bioeng. 82, 253–262. doi:10.1002/bit.10569

PubMed Abstract | CrossRef Full Text | Google Scholar

Balis, U. J., Behnia, K., Dwarakanath, B., Bhatia, S. N., Sullivan, S. J., Yarmush, M. L., et al. (1999). Oxygen consumption characteristics of porcine hepatocytes. Metab. Eng. 1, 49–62. doi:10.1006/mben.1998.0105

PubMed Abstract | CrossRef Full Text | Google Scholar

Ballard, T. E., Wang, S., Cox, L. M., Moen, M. A., Krzyzewski, S. A., Ukairo, O., et al. (2015). Application of a micropatterned co-cultured (MPCC) hepatocyte system to predict preclinical and human specific drug metabolism. Drug Metab. Dispos. 44, 172–179. doi:10.1124/dmd.115.066688

CrossRef Full Text | Google Scholar

Breslin, S., and O’Driscoll, L. (2016). The relevance of using 3D cell cultures, in addition to 2D monolayer cultures, when evaluating breast cancer drug sensitivity and resistance. Oncotarget. doi:10.18632/oncotarget.9935

CrossRef Full Text | Google Scholar

Brody, J. P., and Yager, P. (1997). Diffusion-based extraction in a microfabricated device. Sens. Actuators A Phys. 58, 13–18. doi:10.1016/S0924-4247(97)80219-1

CrossRef Full Text | Google Scholar

Buchwald, P. (2011). A local glucose-and oxygen concentration-based insulin secretion model for pancreatic islets. Theor. Biol. Med. Model. 8, 20. doi:10.1186/1742-4682-8-20

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, T. S., Yu, H., Moore, A., Khetani, S. R., and Tweedie, D. (2013). Meeting the challenge of predicting hepatic clearance of compounds slowly metabolized by cytochrome P450 using a novel hepatocyte model, HepatoPac. Drug Metab. Dispos. 41, 2024–2032. doi:10.1124/dmd.113.053397

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, K., Pfund, W. P., Andersen, M. E., Thomas, R. S., Clewell, H. J., and Lecluyse, E. L. (2014). Development of 3D dynamic flow model of human liver and its application to prediction of metabolic clearance of 7-ethoxycoumarin. Tissue Eng. C Methods 20, 641–651. doi:10.1089/ten.TEC.2013.0562

PubMed Abstract | CrossRef Full Text | Google Scholar

Clevers, H. (2016). Modeling development and disease with organoids. Cell 165, 1586–1597. doi:10.1016/j.cell.2016.05.082

PubMed Abstract | CrossRef Full Text | Google Scholar

Dash, A., Simmers, M. B., Deering, T. G., Berry, D. J., Feaver, R. E., Hastings, N. E., et al. (2013). Hemodynamic flow improves rat hepatocyte morphology, function, and metabolic activity in vitro. Am. J. Physiol. Cell Physiol. 304, C1053–C1063. doi:10.1152/ajpcell.00331.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

De Boer, J., Fritsche, E., Schoeters, G., and Kimber, I. (2015). The European Long-range Research Initiative (LRI): a decade of contributions to human health protection, exposure modelling and environmental integrity. Toxicology 337, 83–90. doi:10.1016/j.tox.2015.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

De Groot, H., Littauer, A., and Noll, T. (1988). “Metabolic and pathological aspects of hypoxia in liver cells,” in Oxygen Sensing in Tissues, ed. H. Acker (Berlin, Heidelberg: Springer), 49–64.

Google Scholar

Dix, D. J., Houck, K. A., Martin, M. T., Richard, A. M., Setzer, R. W., and Kavlock, R. J. (2007). The ToxCast program for prioritizing toxicity testing of environmental chemicals. Toxicol. Sci. 95, 5–12. doi:10.1093/toxsci/kfl103

PubMed Abstract | CrossRef Full Text | Google Scholar

Foy, B. D., Rotem, A., Toner, M., Tompkins, R. G., and Yarmush, M. L. (1994). A device to measure the oxygen uptake rate of attached cells: importance in bioartificial organ design. Cell Transplant. 3, 515–527.

PubMed Abstract | Google Scholar

Gordon, P., Wilkinson, M., Ward, W., and Davies, D. (2014). Bioreactor Chamber. Available at: http://google.com/patents/WO2012156683A1?cl=ar

Google Scholar

Hsu, M. N., Tan, G. D., Tania, M., Birgersson, E., and Leo, H. L. (2014). Computational fluid model incorporating liver metabolic activities in perfusion bioreactor. Biotechnol. Bioeng. 111, 885–895. doi:10.1002/bit.25157

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, C. W., Wang, Y., Tania, M., Zhou, H., Gao, Y., Ba, T., et al. (2013). Computational fluid modeling and performance analysis of a bidirectional rotating perfusion culture system. Biotechnol. Prog. 29, 1002–1012. doi:10.1002/btpr.1736

PubMed Abstract | CrossRef Full Text | Google Scholar

Kostadinova, R., Boess, F., Applegate, D., Suter, L., Weiser, T., Singer, T., et al. (2013). A long-term three dimensional liver co-culture system for improved prediction of clinically relevant drug-induced hepatotoxicity. Toxicol. Appl. Pharmacol. 268, 1–16. doi:10.1016/j.taap.2013.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuntz, E., and Kuntz, H.-D. (2008). Hepatology. Berlin, Germany: Spring- Verlag.

Google Scholar

Lecluyse, E. L., and Alexandre, E. (2010). Isolation and culture of primary hepatocytes from resected human liver tissue. Methods Mol. Biol. 640, 57–82. doi:10.1007/978-1-60761-688-7_3

PubMed Abstract | CrossRef Full Text | Google Scholar

Lecluyse, E. L., Witek, R. P., Andersen, M. E., and Powers, M. J. (2012). Organotypic liver culture models: meeting current challenges in toxicity testing. Crit. Rev. Toxicol. 42, 501–548. doi:10.3109/10408444.2012.682115

PubMed Abstract | CrossRef Full Text | Google Scholar

Leek, R., Grimes, D. R., Harris, A. L., and Mcintyre, A. (2016). Methods: using three-dimensional culture (spheroids) as an in vitro model of tumour hypoxia. Adv. Exp. Med. Biol. 899, 167–196. doi:10.1007/978-3-319-26666-4_10

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, A. P., Rasmussen, A., Xu, L., and Kaminski, D. L. (1995). Rifampicin induction of lidocaine metabolism in cultured human hepatocytes. J. Pharmacol. Exp. Ther. 274, 673–677.

PubMed Abstract | Google Scholar

Macdonald, J. M., Griffin, J. P., Kubota, H., Griffith, L., Fair, J., and Reid, L. M. (1999). “Bioartificial livers,” in Cell Encapsulation Technology and Therapeutics, eds W. M. Kühtreiber, R. P. Lanza, and W. L. Chick (New York, NY: Springer Science+Business Media), 252–286.

Google Scholar

Marga, F., Jakab, K., Khatiwala, C., Shepherd, B., Dorfman, S., Hubbard, B., et al. (2012). Toward engineering functional organ modules by additive manufacturing. Biofabrication 4, 022001. doi:10.1088/1758-5082/4/2/022001

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazzei, D., Guzzardi, M. A., Giusti, S., and Ahluwalia, A. (2010). A low shear stress modular bioreactor for connected cell culture under high flow rates. Biotechnol. Bioeng. 106, 127–137. doi:10.1002/bit.22671

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfund, W. P., Noll, L. A., and Martin, G. A. (2013). In Vitro Tumor Metastasis Model. Available at: https://www.google.com/patents/WO2013017282A1?cl=en

Google Scholar

Pineda, E. T., Nerem, R. M., and Ahsan, T. (2013). Differentiation patterns of embryonic stem cells in two- versus three-dimensional culture. Cells Tissues Organs 197, 399–410. doi:10.1159/000346166

CrossRef Full Text | Google Scholar

Rotem, A., Toner, M., Tompkins, R. G., and Yarmush, M. L. (1992). Oxygen uptake rates in cultured rat hepatocytes. Biotechnol. Bioeng. 40, 1286–1291. doi:10.1002/bit.260401020

PubMed Abstract | CrossRef Full Text | Google Scholar

Shulman, M., and Nahmias, Y. (2013). Long-term culture and coculture of primary rat and human hepatocytes. Methods Mol. Biol. 945, 287–302. doi:10.1007/978-1-62703-125-7_17

PubMed Abstract | CrossRef Full Text | Google Scholar

Ukairo, O., Kanchagar, C., Moore, A., Shi, J., Gaffney, J., Aoyama, S., et al. (2013). Long-term stability of primary rat hepatocytes in micropatterned cocultures. J. Biochem. Mol. Toxicol. 27, 204–212. doi:10.1002/jbt.21469

PubMed Abstract | CrossRef Full Text | Google Scholar

USGS. (2016). Dissolved Oxygen Solubility Tables [Online]. United States Geological Survey, Office of Water Quality. Available at: http://water.usgs.gov/software/DOTABLES/

Google Scholar

Vildhede, A., Wisniewski, J. R., Noren, A., Karlgren, M., and Artursson, P. (2015). Comparative proteomic analysis of human liver tissue and isolated hepatocytes with a focus on proteins determining drug exposure. J. Proteome Res. 14, 3305–3314. doi:10.1021/acs.jproteome.5b00334

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinci, B., Cavallone, D., Vozzi, G., Mazzei, D., Domenici, C., Brunetto, M., et al. (2010). In vitro liver model using microfabricated scaffolds in a modular bioreactor. Biotechnol. J. 5, 232–241. doi:10.1002/biot.200900074

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, E. S., Panko, J., and Paustenbach, D. J. (2009). The European Union’s REACH regulation: a review of its history and requirements. Crit. Rev. Toxicol. 39, 553–575. doi:10.1080/10408440903036056

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, K., Pfeifer, N. D., Kock, K., and Brouwer, K. L. (2015). Species differences in hepatobiliary disposition of taurocholic acid in human and rat sandwich-cultured hepatocytes: implications for drug-induced liver injury. J. Pharmacol. Exp. Ther. 353, 415–423. doi:10.1124/jpet.114.221564

PubMed Abstract | CrossRef Full Text | Google Scholar

Zare-Mehrjardi, N., Khorasani, M. T., Hemmesi, K., Mirzadeh, H., Azizi, H., Sadatnia, B., et al. (2011). Differentiation of embryonic stem cells into neural cells on 3D poly (d, l-lactic acid) scaffolds versus 2D cultures. Int. J. Artif. Organs 34, 1012–1023. doi:10.5301/ijao.5000002

PubMed Abstract | CrossRef Full Text | Google Scholar