Numerical Investigations of Hepatic Spheroids Metabolic Reactions in a Perfusion Bioreactor

Introduction

Drug development researches require strict toxicology that postulates the implication of numerous in vivo and in vitro studies (Esch et al., 2011). However, manufacturing effective drugs is often hampered by weak predictive preclinical experiments to accurately and reliably model the effects of drug compounds inside the human body, leading to low efficacy of drug development processes (Esch et al., 2011; Huh et al., 2012). Due to the role of the liver in the metabolism of drugs and exogenous chemical compounds, miniaturized human cell-based bioreactors mimicking the environment of the human's liver has emerged as a promising alternative tool for complementing the existing drug discovery approaches (Gebhardt, 1992; Groneberg et al., 2002; Baudoin et al., 2011). Apart from vital detoxification effects, the liver has other indispensable functions such as glucogenesis, lipid metabolism, and protein metabolism. Blood plasma protein, particularly albumin, is synthesized mostly by the liver (Gebhardt, 1992; Baudoin et al., 2011). This protein is essential in oncotic pressure maintenance and also participate as a transporter for fatty acids and steroid hormones (Baudoin et al., 2011). Furthermore, many substantial catabolites removal happens in the liver. For instance, highly toxic ammonia breakdowns to urea primarily occur in the liver, which is then excreted as urine by kidneys (Gebhardt, 1992). Majority of the aforementioned functions are carried out by the main parenchymal cell type of the liver, i.e., hepatocytes (Barrett et al., 2006). Metabolism for the vast majority of drugs and chemicals happens in the liver in which toxic metabolites are converted into relatively less harmful compounds for eventual urinary excretion (Barrett et al., 2006). It is noteworthy that acute drug-induced liver failure is one of the major reasons for post-market drug removal. Consequently, persistent exposure of hepatocytes to exogenous compounds makes them optimal candidates for toxicology and drug screening studies.

Thus, engineering liver-on-a-chip platforms for the purpose of drug discovery and toxicity analysis has been the subject of various investigations (Khetani and Bhatia, 2008; Baudoin et al., 2011). On the contrary to current in vitro models utilizing two dimensional (2D) cellular cultures, liver on-chip platforms typically accommodate 3D cell cultures (organoids), are perusable and can utilize human cells (Maguire et al., 2009; Ramaiahgari et al., 2014; Ahmed et al., 2017). Liver slices or biopsies may be the first 3D versatile model which intimately resemble the complex in vivo cytostructure but they rapidly lose their viabilities for about 1 day (van Midwoud et al., 2010; Soldatow et al., 2013). Furthermore, the rate of drug metabolism in liver slices are lower in comparison to the isolated hepatocytes (Ekins et al., 1995; Yu et al., 2001). Therefore, liver slices may not be a suitable alternative to be used for a high-throughput drug screening. Isolated hepatocytes entrapped in collagen sandwiches (Dunn et al., 1991; LeCluyse et al., 1994; Kern et al., 1997), hollow fiber scaffolds (Bao et al., 2011; Skardal et al., 2012; Godoy et al., 2013), hydrogels (Hou et al., 2010, 2011b; Kim et al., 2011) and hepatospheres are well-known 3D models of improving hepatocyte viability and increasing the cell polarization toward creating a more complex network leading to more liver-like models. The use of hepatic spheroids has been shown to maintain liver viability and liver-specific functions such as albumin production, urea synthesis for longer time intervals respect to the other aforementioned cultures (Ramaiahgari et al., 2014). In spheroids, cells are closely in contact with each other and produce their own extracellular matrix (ECM), which result in retaining most of their cell-cell and cell-ECM contacts in comparison to the majority of the three-dimensional culture methods (Godoy et al., 2013).

Despite the efforts in making liver-on-a-chip platforms, there are still no guidelines for determining the proper design parameters, including perfusion rate, number of cells, size of the spheroids, etc. One of the critical challenges of 3D hepatic cultures is supplying adequate oxygen and nutrition to the cell surface (Godoy et al., 2013; Hsu et al., 2014). The rate of oxygen consumption of hepatocyte in comparison to the other cell types is high due to the high metabolic activity of these cells (Ramaiahgari et al., 2014). Also, oxygen solubility in culture media at 37°C is low (Yu et al., 2001). Supplying sufficient oxygen to the hepatic spheroids becomes more severe because the abovementioned oxygen transport obstacles are accompanied by diffusion dominant transport mechanism inside the spheroid, which becomes particularly serious in spheroids with a larger diameter. Consequently, determining the oxygen distribution within liver-on-a-chip models comprising cells in the form of hepatic spheroids has remained an important factor, which is extremely hard to optimize experimentally.

Dynamic systems such as the perfusion bioreactors can be a great help since perfusion continuously supplies oxygen and nutrition to the cells and simultaneously removes metabolic waste, leading to improve viability and life span of the hepatocytes. However, perfusion of media inside the bioreactor will impose undesirable flow-induced shear stress on the cells, which can be detrimental or even cause cell detachment (Hou et al., 2011a). Therefore, there must be a balance between the culture medium flow rate and the oxygen supply for the cells (Hsu et al., 2011; Dumont et al., 2018).

Till now, several mathematical methods have been developed to investigate the effect of shear stress and mass transfer at the cellular surface (Curcio et al., 2007; Hsu et al., 2014; Bhise et al., 2016; Khakpour et al., 2017). Transport of oxygen inside the spheroids with different sizes cultured in a rotating wall system under both gas permeable and gas impermeable conditions was investigated by Curcio et al. (2007). Their results showed that spheroids with sizes larger than 200 μm were suffered from oxygen depletion. Khakpour et al. (2017) investigated the oxygen transfer to hepatospheres inside the hollow fiber membrane bioreactor, and the effect of various parameters such as perfusion rate, porosities, spheroid sizes, etc. were investigated. Optimal flow condition and oxygen concentration gradient were calculated numerically in a bioreactor containing hepatic spheroids encapsulated in hydrogel under continuous flow by Bhise et al. (2016). Spheroids in their bioreactor showed considerable viability and functionality over 30 days of culture periods (Bhise et al., 2016). Their results obtained from analyzing and monitoring the concentration of the secreted biomarkers such as albumin, alpha-1 antitrypsin (A1AT) indicated that the hepatic spheroids in their bioreactor remained functional over 30 days of culture periods (Bhise et al., 2016). In another study, Hsu et al. developed a mathematical model of a bioreactor in which a monolayer of HepG2 cells was cultured at the bottom (Hsu et al., 2014). Hydrodynamic equations, along with mass transfer equation were solved to obtain flow characteristics and metabolites concentration, respectively. Despite abovementioned efforts to model the shear stress and oxygen levels at the cellular surface, there is still lack of enough literature to mathematically incorporating intracellular metabolic activities of hepatocytes to their related metabolites in the extracellular environment via governing mass transfer in a bioreactor to quantitatively investigate the effect of different parameters. Such models provide an analytical platform for design, optimization of the bioreactor, and help researchers in perceiving better the interplay between bioreactor design parameters and liver-specific functions. This mathematical system also provides rapid analyses which can be helpful in the disciplines like forensic bioengineering where immediate and accurate results of trace sampling are crucial. Results from the presented numerical platform can also be used in better interpreting and certifying the outcomes of investigations. Various drug pathways like depressants, stimulants, analgesics, and psychotomimetics can be added to this numerical liver-on-a-chip; and hence, providing the results which can be used as a guide for testifying the related findings of drug abuse and forensic toxicology investigations. Also, the numerical system can be useful in prediction and estimation on the concentration of the species wherein some cases, only small amount of samples might be gathered in the victim scene. High-speed result of different concentrations of the substance can be provided using this numerical liver-on-a-chip platform.

In this study, we developed a numerical model based on the 3D culture of hepatic encapsulated spheroids in miniaturized chambers under continuous perfusion (Figure 1A). Optimized working conditions of the bioreactor were obtained in terms of sufficient oxygen delivery to the cellular surface and protecting the cells from high levels of shear stress. In the first step, fluid flow distribution inside the bioreactor is calculated. Transport phenomena are considered in different sections of the bioreactor, i.e., free flow, porous flow, spheroids, and single cell-level, providing a powerful tool mimicking a liver-on-a-chip system prior to any fabrications. Through the present computational multiscale model, researchers can test multiple designs and observe the effect of various parameters with real-time monitoring of the pivotal phenomena of dynamic nature of model under investigation in cost and time-effective manner. Additionally, mass transfer equations are incorporated into hydrodynamics equations to measure the concentration of different metabolites. Oxygen concentration and other main metabolite functions of the hepatocytes such as albumin production, glucose, and glutamine consumptions, ammonia and urea generation are obtained in subcellular level inside the implicitly modeled hepatocytes aggregated as a spheroid configuration using a mathematical model. Since sufficient oxygen availability is a necessary condition for the cellular viability, metabolites concentrations were calculated by incorporating the effect of oxygen tension on metabolite concentration. The above metabolites are diffused into or out of the spheroids through the hydrogel. The concentration distribution of each metabolite obtained in each spheroid and the whole domain. The model shows a great promise as the simulation results calculated for albumin over 4 weeks of culture showed good compatibility with experimental data of liver albumin production level. Moreover, in the final section, a parametric study has been carried out and the effect of different prevailing parameters on oxygen transfer to hepatocyte spheroids are investigated. Oxygen concentration is chosen as the main optimization criterion among other investigated metabolites for the parametric study because hepatocytes are highly susceptible to the oxygen level and they drastically lose their liver-specific functionalities under the effect of inadequate oxygen concentration.

FIGURE 1

Figure 1. General configuration of liver spheroid-on-a-chip model. (A) Schematic of the bioreactor. Bioreactor composed of two spheroid chambers and a media channel. Spheroids were placed at the certain position inside the hydrogel. Flow direction is depicted with blue arrows. Flow enters from the middle channel and permeates through both upper and lower hydrogel chambers. The bioreactor was connected to the reservoir, and the whole system was run using a peristaltic pump. (B) Dimensions of the different section of bioreactor. The dimension of the upper and lower sections is similar.

Methodology

The methodology used for modeling distribution of mass, flow, and concentration is illustrated in the following sections. Since the problem under investigation is composed of multi-physics, governing equations of the fluid flow and concentrations are coupled. The solution of continuity, hydrodynamic, and mass transfer equations are obtained using 3D finite element method. The Galerkin method was used in order to approximate the non-linear partial differential governing equations with a system of ordinary differential equations. The equivalent weak formulations of the Navier–Stokes, and convection-diffusion equations were achieved by applying the boundary conditions, followed by discretization of the weak forms in a finite-dimensional subspace. An Eulerian approach was chosen for the one-dimensional, time-dependent domain.

The geometry of the device is composed of three different chambers; two side chambers containing hydrogel in which hepatic spheroids are encapsulated and one middle chamber through which medium and nutrition flow. Schematic view of the bioreactor, including dimensions of different parts is depicted in Figure 1B. The overall bioreactor was connected to the reservoir, and the whole system was run using a peristaltic pump (Figure 1A).

Modeling methodology of the presented simulation is shown in a schematic configuration in Figure 2 denoting the coupling between governing equations describing the multiphysical nature of the problem. Fluid flow dynamics is coupled with mass transport and reactions inside the bioreactor. First, equations of continuity and momentum are solved, and fluid velocity distribution inside each section of the bioreactor is obtained. Then, for each concentration, mass transfer equation is computed along with related metabolic reactions.

FIGURE 2

Figure 2. Block diagram of the modeling process. First, fluid flow equations were solved, and velocity distribution and wall shear stress inside the bioreactor were obtained. Next, mass transfer equations were solved for oxygen, albumin, glucose, glutamine, ammonia, urea, and their related metabolites.

Fluid-Dynamics Model

Flow containing oxygen and nutrition enters via the inlet port of the bioreactor at a prescribed volumetric flow rate. The top and the bottom surface of the fluidic channel, depicted in Figure 3A, are in contact with hydrogel. The hydrogel in this model is considered as a uniform, isotropic porous medium. Oxygen and nutrition will permeate via constant perfusion throughout the hydrogel due to its porous characteristics. Thereupon, the governing equation of fluid flow, mass transfer, and cell kinetics are coupled. The hydrodynamic modeling of the bioreactor is divided into two parts, i.e., free flow and porous flow regimes, based on their governing equations (Figure 3A and Figure S1). Inlet volumetric flow rate is 1 mlit/min, declaring that Reynolds number (Re) <2000, and the fluid flow regime is laminar inside the channel (For further explanation, please check Supplementary Materials). Based on the assumptions mentioned above, the governing equation demonstrating the flow of fluid and two hydrogel chambers will be Navier-Stokes and Stokes-Brinkman equations, respectively (Truskey et al., 2004; Figure S1A). Equations of continuity, momentum for free and porous flow regimes can be written as:

$\begin{array}{lll}\nabla ·\vec{\text{u}}& =0& (1)\end{array}$

$\begin{array}{lll}\rho (\vec{\text{u}}·\nabla \vec{\text{u}})& =\nabla ·\left[-\text{pI}+\mu (\nabla \vec{\text{u}}+{(\nabla \vec{\text{u}})}^{\text{T}})\right]& (2)\end{array}$

$\begin{array}{l}\frac{\mu }{K}\overrightarrow{u}=\frac{1}{\rho }\nabla ·[-pI+\frac{\mu }{\epsilon }\left(\nabla \overrightarrow{u}+{\left(\nabla \overrightarrow{u}\right)}^T\right)\\ -\left(\nabla ·\overrightarrow{u}\right)I]& (3)\end{array}$

where ρ, μ, and p are fluid density, viscosity and pressure, respectively. K is the permeability, and ε is porosity in the porous flow regime, i.e., representing hydrogel and spheroids sections. $\vec{u}$ denotes three-dimensional fluid velocity vector $\vec{u}=$(u, v, w). It is assumed that the solution is an incompressible, isothermal, Newtonian fluid. Since all of the nutrition inside the cell culture medium is dissolved in aqueous-based liquid, all of the fluid properties are assumed to be similar to the properties of water at the same temperature. The schematic of the boundary conditions used in the present simulation is depicted in Figure 3B. The above equations are coupled via two stress boundary conditions at the interface of the fluidic and two hydrogel chambers (For further explanation, please check Supplementary Materials). Here, in this model, hepatic spheroids with different diameters of 100 to 400 μm are located at the specific places inside the hydrogels. Spheroids are modeled again as a uniform, homogeneous porous medium. Boundary condition on each wall of the bioreactor is assumed to be zero velocity, representing no-slip boundary condition (Figure 3B). Since each spheroid is modeled as a porous sphere, hence, the velocity distribution inside the spheroid is obtained again by solving Stokes-Brinkman equation similar to the hydrogel compartments. Porosity and permeability of the spheroids are calculated based on their cellular density, which will be described further in details. Constant inlet flow velocity $(\overrightarrow{u}=u_{in})$at the inlet and outflow boundary condition $(\partial \vec{u}/\partial n=0)$, where n is a normal vector, at the outlet of the fluid flow channel are assumed for the hydrodynamic boundary conditions. The boundary condition at the outlet implicates that flow is fully developed (Figure 3B).

FIGURE 3

Figure 3. Concept schematic of the modeled geometry. (A) Schematic front view (Z-Y) of the bioreactor. The bioreactor is composed of two hydrogel chambers containing spheroids and a central main chamber. Medium flows into the main-stream and then spreads to the hydrogel compartments. (B) The imposed boundary conditions used for current simulation. Constant inlet velocity and constant initial concentration for species are set as the inlet boundary conditions for hydrodynamics and mass transfer equations, respectively. No slip (zero velocity) and no flux boundary conditions are applied on the boundaries marked by black dashed lines, since there is no penetration at the solid side walls of the bioreactor. The dashed red lines show the interface boundaries between free and porous flow domains and hydrogel-porous spheroid domain where species and stress continuity boundary conditions are applied. Outflow boundary condition is set at the outlet.

The shear stress on the cell surface can be measured by

$\begin{array}{ll}{\tau }_s=\mu \frac{\partial \vec{\text{u}}}{\partial \text{n}}& (4)\end{array}$

n is a normal vector of the cell surface. Parameters used for the simulation of fluid flow are given in Table 1.

TABLE 1

Table 1. Material properties and constants used in hydrodynamic equations.

Mass Transfer Modeling

In this study, the distribution of each metabolite's concentration inside the spheroids, hydrogel and in the main section are obtained using mass transfer equation, i.e.,

$\begin{array}{ll}\frac{\partial {\text{C}}_i}{\partial \text{t}}+\vec{\text{u}}·\nabla {\text{C}}_{\text{i}}=\nabla ·({{\text{D}}_{\text{i}}\text{C}}_{\text{i}})& (5)\end{array}$

for the hydrogel compartments and the main section and

$\begin{array}{ll}\frac{\partial {\text{C}}_{\text{i}}}{\partial \text{t}}+\vec{\text{u}}·\nabla {\text{C}}_{\text{i}}=\nabla ·({{\text{D}}_{\text{i}}\text{C}}_{\text{i}})+{\text{R}}_{\text{i}}& (6)\end{array}$

for the spheroid parts, where C_i is the concentration of species, D_i represents species diffusion coefficient inside substrate solution, $\vec{\text{u}}$ is the velocity vector obtained from the hydrodynamic equations, i corresponds to the specific species and, R_i indicates the reaction rates of each species which specifies respected metabolic reactions of the cells (Figure S1B). The diffusion coefficient of metabolites inside the hydrogel is correlated to its counterparts in the medium substrate by an equation, which will be defined later. It is shown that the metabolite diffusion coefficient inside the hydrogel is increasing directly with porosity and has an inverse relation with the tortuosity of the hydrogel (Hou et al., 2010). The species modeled in this study are oxygen, albumin, glucose, glutamine, ammonia, and urea. At the inlet, constant inlet concentration C_i=C_i.in and outflow boundary condition at the outlet;

$\begin{array}{ll}-n·{{\text{D}}_{\text{i}}\text{C}}_{\text{i}}=0& (7)\end{array}$

respect to each metabolite is considered as boundary conditions for mass transfer equation. Respected values for each metabolite is given in each section. No flux boundary condition is assumed for all walls, as well (Figure S1A). Parameters used for the simulation of mass transfer are given in Table 2.

TABLE 2

Table 2. Parameters constants used in the mass transfer modeling.

The hepatic cell density of the spheroid is a function of its porosity where a high packed spheroid containing a large number of hepatic cells will have the smaller void fraction or lower porosity; hence

$\begin{array}{ll}{\rho }_{\text{cell}-\text{spheroid}}=\frac{1-{\epsilon }_{\text{sph}}}{{\text{V}}_{\text{cell}}}& (8)\end{array}$

ε_sph is the spheroid porosity, and V_cell is the volume of the cells (Khakpour et al., 2017). Furthermore, the effective diffusion coeffcient for each metabolite inside the porous media, D_eff, is corresponded to its respected diffusion coefficient in the liquid phase,

$\begin{array}{ll}{\text{D}}_{\text{eff}.\text{i}}=\frac{\epsilon }{\tau }{\text{D}}_{\text{i}}& (9)\end{array}$

which can be written for hydrogel sections and spheroids. τ and ε are represented tortuosity and porosity, respectively. Here Bruggeman (Vafai, 2015; Equation 10) with ε_sph> 0.2 and Wakao-Smith (Merdaw et al., 2010; Equation 11) relations are used for calculating the tortuosity of the hepatic spheroid and hydrogel parts, respectively, i.e.:

$\begin{array}{ll}{\tau }_{\text{sph}}=\frac{1}{{\epsilon }_{\text{sph}}^{0.5}}& (10)\end{array}$

$\begin{array}{ll}{\tau }_{\text{hyd}}=\frac{1}{{\epsilon }_{\text{hyd}}}& (11)\end{array}$

The permeability of media in the hepatic spheroids is modeled using Carman Kozeny model (Shipley and Waters, 2011), in which hepatocytes are considered as a sphere with the diameter of d_hep = 16 μm,

$\begin{array}{ll}{\text{K}}_{\text{s}\text{p}\text{h}}=\frac{{\epsilon }_{\text{s}\text{p}\text{h}}^3{\text{d}}_{\text{h}e\text{p}}^2}{180\times {(1-{\epsilon }_{\text{s}\text{p}\text{h}})}^2}& (12)\end{array}$

Also, the hydrogel permeability is calculated via Merdaw approximation (Merdaw et al., 2010) in which permeability of water inside the hydrogel was modeled. For the presented study, the fluid properties of culture media considered to be similar to the fluid properties of the water at 37°C. Hence, the Merdaw approximation for the presented model (Merdaw et al., 2010) becomes;

$\begin{array}{ll}{\text{K}}_{\text{hyd}}=\frac{\text{DM}\mu }{\rho {\text{R}}_{\text{g}}\text{T}}{\epsilon }_{\text{hyd}}^2+\frac{{\text{d}}_{\text{hyd}}^2}{11.25}\frac{{\epsilon }_{\text{hyd}}^3}{180\times {(1-{\epsilon }_{\text{hyd}})}^2}& (13)\end{array}$

where D, μ, M, and ρ are the water self-diffusivity, viscosity, molecular weight, and density at 37°C, respectively. R_g is the universal gas constant, T is temperature, ε_hyd is hydrogel porosity, and d_hyd is mean pore diameter of the hydrogel. These parameters are given in Table 1.

Oxygen Concentration

Inlet oxygen concentration inside the medium is obtained via Henry's law, which is relating dissolved oxygen concentration inside the medium to the oxygen partial pressure at 37°C and atmospheric pressure condition (Barrett et al., 2006; Khetani and Bhatia, 2008):

$\begin{array}{ll}{p}_{o_2}=k_H{C}_{o_2}& (14)\end{array}$

where p_O2 denotes partial pressure of the oxygen, k_H is Henry's constant and C_O2 represents the concentration of the oxygen inside the solution (Maguire et al., 2009). This concentration is assumed to be constant at the inlet since the PDMS (Polydimethylsiloxane) is permeable to the oxygen and oxygen pressure will always equilibrate with oxygen existed in the surrounding environment, i.e., inside the incubator. It is possible to modify the inlet oxygen concentration by changing the oxygen supply level of the incubator. The effect of this parameter will be investigated in the parametric study section. Three layers of the bioreactor are fixed via two PMMA (Polymethylmethacrylate) sheets; therefore, there is no oxygen diffusion from the top and bottom surface of the bioreactor.

Based on the previous studies (Gebhardt, 1992), it is shown that the oxygen consumption rate of hepatocytes depends on the available oxygen concentration in the vicinity of the cell surface (Gebhardt, 1992). Oxygen consumption rate is governed by the first order Michaelis-Menten equation, which is the best-known model relating the rate of hepatic oxygen consumption to the oxygen concentration inside the medium (Gebhardt, 1992):

$\begin{array}{ll}{\text{R}}_{{\text{o}}_2}={\text{N}}_{\text{cell}}\times \frac{{\text{V}}_{\text{m}}\times {\text{C}}_{{\text{o}}_2}}{{\text{K}}_{\text{m}}+{\text{C}}_{{\text{o}}_2}}& (15)\end{array}$

where N_cell is cell number for each spheroid, V_m is maximum oxygen uptake rate per cell and K_m denotes Michaelis-Menten constant indicating oxygen partial pressure at the point where oxygen uptake rate reaches half of its maximum amount. Parameters used in this section are given in Table 2.

Albumin Concentration

It is mentioned that the albumin production of HepG2 cells is related to the local oxygen concentration at the cell surface and shear stress sensed by the cell due to fluid flow. Albumin production rate can be modeled via

$\begin{array}{l}{\text{R}}_{Alb}={\text{N}}_{cell}\times ({\text{e}}_{1}ln({\tau }_{\text{s}})+{\text{e}}_2+{\text{e}}_3{\text{C}}_{{\text{o}}_2})\times ({\text{C}}_{{\text{o}}_2}>0.03)\\ \times ({\text{C}}_{Glu}>0)\times ({\text{C}}_{Glut}>0)& (16)\end{array}$

where R_Alb is albumin production rate of hepatic cells and τ_s is the shear stress felt by cells (Hsu et al., 2014). The last three terms in the above equation indicating that albumin production requires sufficient oxygen and nutrition concentrations. Therefore, in a situation where the oxygen concentration is below 0.03 mol/m³, or concentration of glucose or glutamine is zero, albumin production will be zero. Constants used in Equation (16) are given in Table 3.

TABLE 3

Table 3. Constants used in mass transfer modeling of Albumin, Glucose, and Glutamine (Hsu et al., 2014).

Glucose and Glutamine Consumption

HepG2 cells consume Glucose and Glutamine of the nutrition media. It is shown that (Barrett et al., 2006) Glucose consumption rate of the HepG2 cells is inversely related to the oxygen availability of the cells where at high oxygen concentration, Glucose consumption of the HepG2 is low while at lower oxygen tension, it is found that these cells will take up more Glucose (Barrett et al., 2006). Therefore, the Glucose consumption rate of the HepG2 can be written as

$\begin{array}{ll}{\text{R}}_{\text{Glu}}={\text{N}}_{\text{cell}}\times ({\text{e}}_4{\text{C}}_{{\text{o}}_2}+{\text{e}}_5)\times ({\text{C}}_{\text{Glu}}>0)& (17)\end{array}$

where e₄ and e₅ are given in Table 3.

Glutamine consumption rate (Figure 4A) of HepG2 is depended on the ammonia concentration inside the media. It is shown that (Baudoin et al., 2011) at the ammonia concentration of <10 mM, the glutamine consumption rate of the cells is 1.43 × 10⁻¹⁶ mols⁻¹cell⁻¹. At ammonia concentration of higher than 10 mM, the glutamine consumption rate will be 2.85 × 10⁻¹⁶ mols⁻¹cell⁻¹ (Baudoin et al., 2011). Hence, Glutamine consumption rate can be obtained via

$\begin{array}{l}{\text{R}}_{\text{Glut}}={\text{N}}_{\text{cell}}\times \left[\left({\text{e}}_6\times \left({\text{C}}_{{\text{NH}}_3}<10\right)\right)+\left({\text{e}}_7\times \left({\text{C}}_{{\text{NH}}_3}>10\right)\right)\right]\\ \times \left({\text{C}}_{\text{Glut}}>0\right)& (18)\end{array}$

The above constants, e₆, and e₇, are given in Table 3.

FIGURE 4

Figure 4. Schematic diagram of the urea cycle, ammonia detoxification, and glutamine synthesis inside a hepatocyte. (A) Periportal hepatocytes receive blood with high oxygen supply and nutrition from the portal vein and hepatic artery. In these hepatocytes, the ammonia detoxification is represented by the synthesis of the urea. Perivenous hepatocytes, aligned near the central vein, are supplied with less nutrition and low oxygenated blood, the ammonia detoxification here is accompanied by glutamine synthesis. (B) Ammonia enters the hepatocytes urea cycle as carbamoyl phosphate, passes through four enzymatic reactions and at final stage urea excretes to the blood.

Ammonia Production and Consumption

Ammonia and glutamate are generated by degradation of glutamine. Also, in the following steps, Glutamate dehydrogenase can lead to the formation of another molecule of ammonia from glutamate (Figure 4A; Baudoin et al., 2011). Therefore, during the degradation process of the glutamine, two molecules of ammonia will be generated. Based on the mentioned relation between glutamine degradation and ammonia formation, the rate of ammonia production is twice the rate of glutamine consumption (Hsu et al., 2014). Also, liver cells convert the available ammonia to urea through the urea cycle (Figure 4B). The detailed urea cycle reactions are given in Supplementary Materials. Hence,

$\begin{array}{ll}{\text{R}}_{{\text{NH}}_3}={2\times \text{R}}_{\text{Glut}}\times ({\text{C}}_{{\text{o}}_2}>0.03)-{\text{R}}_{{\text{NH}}_3.\text{c}}& (19)\end{array}$

where the last term on the right-hand side of the equation denotes ammonia conversion to urea, which is modeled based on four main enzymatic reactions in the urea cycle (Kuchel et al., 1977). The detailed enzymatic reactions used for modeling urea generation are given in Supplementary Materials in detail.

Results and Discussion

Validation of the Numerical Method

The oxygen concentration profile along the spheroid diameter is plotted for different mesh numbers (Figure 5A), providing the mesh size independency of the obtained results of the simulation. From this figure, it can be inferred that results obtained from the grid with 589,006 elements are suitable for further studies and finer grid numbers do not show any significant discrepancies. Also, to validate the finite element code used for the present study, the generated code was used for replicating of the previous study remodeling carried out by Tourlomousis and Chang (2016b). Computed results of the fluid velocity inside the channel are compared with those reported by Tourlomousis and Chang (2016b) showing excellent compatibility (Figure 5B). Also, the results of the predicted metabolites participating in the urea cycle were modeled and compared with the results given by Kuchel et al. (1977). Three random metabolite concentrations are chosen and depicted in Figure 5C. Obtained results have <10% differences in comparison to the data given by Kuchel et al. (1977). For validating the oxygen concentration data predicted by the present simulation, calculated oxygen concentrations are compared with the reported experimental data for three different spheroid sizes, Figure 5D, representing a good agreement between computed oxygen concentration and experimental data given by Murphy et al. (2017).

FIGURE 5

Figure 5. Validation and Grid study for evaluation of the presented code used in this study. (A) Oxygen concentration profile along the spheroid diameter for different mesh numbers. (B) Comparison of computed results of the fluid velocity inside the channel to those reported by Tourlomousis and Chang (2016b). (C) Comparison of the predicted metabolites participating in the urea cycle with the results given by Kuchel et al. (1977) (D) Oxygen concentrations are compared with the experimental data obtained by Murphy et al. (2017) for three different spheroid sizes.

Hydrodynamic and Shear Stress

The main purpose of the perfusion system inside each microbioreactor is providing nutrition and oxygen to the cultured cells (Tanaka et al., 2006). Availability of a certain amount of oxygen is of great importance for cells to have normal metabolic activity. According to the literature (Evenou et al., 2010; Hsu et al., 2014), hepatocytes consume the highest amount of oxygen in comparison to the other cell types due to their high metabolic activities. Based on the low oxygen solubility and diffusivity in the culture medium, the concentration of oxygen is a critical factor in designing the bioreactors. Perfusion of the culture medium can enhance and control the oxygen supply to the cells by regulating the medium flow rate, which in turn imposes shear stresses on them. It is mentioned that shear stress affects metabolic functions and morphological configuration of the hepatocytes (Tanaka et al., 2006). Accordingly, there should be a balance between oxygen delivery to the cells and perfusion rate of the medium flow since high perfusion rate of culture media implies the excess amount of shear stress to the cells, which is also detrimental. Several studies indicated that hepatocytes viability under high shear stress conditions fell below that of hepatocytes cultured in static conditions (Park et al., 2008). It should be mentioned that the shear stress would affect hepatocyte functionality at values > 0.03 Pa (Tilles et al., 2001). Upon the hepatocyte's sensitivity to oxygen concentration values, optimization and parametric study carry out here for perfusion bioreactor has considered different aspects affecting oxygen supply to the hepatocytes. Perfusion rate is one of the efficacious parameters, where high flow rates have been proved to provide ample oxygen supply to the cells and on the other hand, imposes deleterious shear stress on the cell surface. Thereupon, the optimization in terms of the abovementioned parameters, which will be discussed further in the parametric study, signifies a delight balance between oxygen supply and shear stress.

Culture medium velocity distribution is displayed in Figure 6A. The arrows are aligned from the inlet to the outlet of the bioreactor, tangential to the velocity of each point which length is corresponding to the velocity magnitude. It can be observed that the velocity vector at the centerline of the bioreactor is higher than the rest of the vectors, indicating a parabolic profile for the velocity inside the main channel of the bioreactor. Fluid flow velocity inside the hydrogel part of the bioreactor has small values, leading to low-value shear stresses applied on the spheroids hosted inside the hydrogel. Thus, the perfusion rate of culture medium for encapsulated cells can be higher than those cells directly exposed to the fluid flow domain. Shear stress distribution inside the bioreactor and shear stresses imposed on spheroids are shown in Figure 6B. The allowable shear stress tolerated by cells or required as a desirable factor inducing mechanotransduction effects is an intricate issue depending on many factors including cell type, cultural condition, culture period and elevated hepatic functions (Wang et al., 2010; Ebrahimkhani et al., 2014). The critical information from shear stress distribution (Figure 6B) is that shear stresses inside the hydrogels compartments and on spheroids are at values below the reported threshold of 0.03 Pa (Tilles et al., 2001). Such analysis is indispensable in perfusion bioreactors as cells are prone to a high level of shear stress.

FIGURE 6

Figure 6. Results of the numerical simulation of fluid flow inside the bioreactor. (A) Fluid velocity distribution inside the bioreactor. Red vectors show the direction of the flow, which their length is correlated to the velocity magnitude (B) Shear stress distribution inside the whole bioreactor system. Shear stress imposed by the perfusion of the culture medium is below the threshold.

Oxygen Concentration

The bioreactor is perfused with the medium under the in vitro standard conditions, i.e., 20% oxygen partial pressure at 37°C. Since spheroids are considered to be non-vascularized, it is possible to increase inlet oxygen concentration to prevent cells from oxygen deprivation. Such higher oxygen concentration is not detrimental to the hepatocytes since the oxygen concentration in the liver is between 4 and 9% (42–91 μmol/L) (Park et al., 2008). The effect of higher inlet oxygen concentration will be investigated in the parametric study section. Figure 7A shows oxygen concentration distribution inside the bioreactor after 5 min of perfusion, indicating that oxygen concentration in the vicinity of the hepatic spheroids is sufficient for the cells. Oxygen insufficiency at the center part of the spheroids is common since the transport mechanism for oxygen delivery to hepatocytes inside the spheroid is mainly based on diffusion. Due to the effect of the perfusion, oxygen supply at spheroid surfaces is improved, but it still is possible that in some part of the hepatic spheroids, oxygen level falls below the minimum oxygen supply necessary for the cells. In Figure 7B oxygen concentration distribution inside the first, fourth, and seventh spheroid at t = 5 min is presented, denoting that in the whole spheroids adequate oxygen supply is available for the hepatocytes. Oxygen is continuously consumed by the cells; therefore, oxygen tension might fall below the critical value [0.03 mol/m³ (Jungermann and Kietzmann, 2000)]. Oxygen distribution throughout the bioreactor and inside the hepatic spheroids at t = 10 min are displayed in Figures 7C,D, respectively which reflects that oxygen tension is also in an adequate level for the cells showing that perfusion can maintain sufficient oxygen supply to the hepatocytes resulting in preserving hepatocytes functionality. These results clarify that the hepatocytes aggregates inside the perfusion bioreactor are under the normal oxygenation. Greater oxygen tension can be supplied by applying higher oxygen level at the inlet or increasing the perfusion rate in this bioreactor since hydrogel in the cellular section of the bioreactor preserves the cells from a high level of shear stress exposure. A higher level of oxygen at the inlet can be obtained by increasing the oxygen supply inside the incubator. Furthermore, oxygen penetration from the PDMS sections of the bioreactor incessantly keeps oxygen tension at the maximum level, which also has a profound effect on sustaining hepatocyte-specific functions. The effect of this parameter will be analyzed in the parametric study.

FIGURE 7

Figure 7. Results of the Oxygen concentration distribution inside the bioreactor. (A) The oxygen concentration in the bioreactor at t = 5 min. (B) Oxygen distribution inside the 1st, 4th, and 7th spheroids, t = 5 min (C) The oxygen concentration in the bioreactor at t = 10 min. (D) Oxygen distribution inside the 1st, 4th, and 7th spheroids t = 10 min.

Albumin Concentrations

Albumin secretion analysis is generally one of the critical factors used by many studies in evaluating the functionality of the hepatocyte bioreactor (Glicklis et al., 2004; Anada et al., 2012; Bhise et al., 2016; Ahmed et al., 2017). It is proved that (Dullien, 2012) providing higher oxygen level to the hepatocytes will increase their albumin production rates whereas increasing shear stress levels imposed on cells by medium flow will decrease metabolic functionality and thereby the albumin production rate (Tan et al., 2013; Hsu et al., 2014). Hence, albumin production rate under the influence of two effective parameters, i.e., oxygen tension at the vicinity of the cells and shear stress sensed by them, is analyzed here via applying the method developed by Hsu et al. (2014). Results of the hepatocyte albumin production mapped on the central surface of the bioreactor at t = 10 min is displayed in Figure 8A. From this figure, it can be observed that albumin concentration is raised from zero to 2 × 10⁻¹⁴ mol/m³ in about 10 min of the culture duration. Also, the cumulative concentration of the albumin secretion at the outlet of the bioreactor over 7 days of the culture period is calculated (Figure 8B) and results are compared with experimental data of Ahmed et al. (2017), representing good agreement with the experimental data, though experimental results are slightly higher than the predicted value obtained by numerical simulation. This discrepancy might be based on the simplification made for generating the albumin production rate equation (Equaton 16), in which some other parameters affecting albumin synthesis in hepatocytes like nutritional and hormonal factors are ignored (Rothschild et al., 1977; Hutson et al., 1987).

FIGURE 8

Figure 8. Results of the Albumin, Ammonia, Glucose, Glutamine, and Urea concentration inside the bioreactor. (A) Albumin concentration distribution inside the bioreactor at t = 10 min. (B) Cumulative albumin concentration at 7th day, showing a good agreement with the experimental data calculated by Ahmed et al. (2017). (C) Calculated ammonia concentration distribution. (D) Predicted glucose concentration distribution, (E) Glutamine concentration distribution. (F) Predicted urea production distribution at day 2 inside the bioreactor.

Glutamine, Glucose, Ammonia, and Urea Concentrations

Ammonia synthesis is one of the hepatocytes critical functions in which hepatocytes eliminate ammonia from the circulating system through the urea cycle leading to a non-toxic compound called urea. Figure 8C shows the predicted ammonia concentration distribution inside the incubator. Ammonia is removed via the urea cycle and generated in glutamine synthesis. Glucose is a necessary compound in culturing the hepatocytes. The distribution of glucose concentration in Figure 8D is also indicating the lower level of glucose respect to the inlet concentration. The negligible decrease of glucose concentration implies the low level of glucose consumption rate of the cells. This small low diminution of glucose is compensated by the relatively high concentration of the glucose at the inlet, providing sufficient glucose level for the cells. Distribution of the glutamine concentration is shown in Figure 8E, demonstrating the glutamine consumption at the spheroid surface. The glutamine level falls from its initial value due to the cellular intake. This small gradual decrease continues during the culture period and controls mainly by ammonia concentration inside the bioreactor (Equation 18). Urea production through urea synthesis is shown in Figure 8F denoting the small value of urea generated during 2 days of the culture period. Ammonia to urea conversion happens at the subcellular level, which is incorporated here for each hepatocyte and modeled through four enzymatic reactions (Supplementary Materials). Despite various assumptions for obtaining the aforesaid metabolites, predicted results provide precious insights into the different parameters for the researchers before performing any experiments.

Both calculated glucose and glutamine distributions are obtained based on the simplified methods developed by Hsu et al. (2014). More detailed metabolic relations should be included to have a more accurate and sensitive prediction, similar to the method used here for urea cycle prediction, so that their respected subcellular functions are incorporated. Therefore, glutamine and glucose distribution though providing a general reasonable overview and initial estimation, may not be considered as a sensitive tool in the evaluation and analysis of bioreactor design. On the contrary, results of oxygen, albumin, ammonia, and urea concentration are showed the promising indexes which can be utilized in assessing the efficiency and performance of the perfusion bioreactor. In this regard, based on the results of the oxygen distribution, a parametric study is conducted in the next section, and the effect of various design and working parameters are analyzed. Providing such numerical studies is of great importance in analyzing and optimizing of the bioreactor and can be beneficial for drug development industries by reducing time and design cost of the predesign process.

Parametric Study

In this section, the result of changing different parameters affecting the efficiency of the perfusion bioreactor is investigated. Several parameters like perfusion rate, intra spheroid porosity, hydrogel porosity, spheroid diameter, saturated dissolved oxygen concentration, V_m and K_m in Michaelis-Menten kinetics are varied individually over the specified range while keeping other parameters at the fixed baseline values mentioned in Table 4. Porosity, considered as a microstructural property of the porous media, modifies the transport phenomena (Equation 3), effective diffusion coefficient (Equation 9) and permeability (Equations 12 and 13). Evaluation of in vitro hepatic spheroid porosity is quite challenging, and there is little information about that (Moussy, 2003; Khakpour et al., 2017). Additionally, the hepatic spheroid porosity is not a permanent value and varies as cells proliferate and rearrange during the culture period. Several methods have been developed for estimating the porosity of packed beds of monodispersed rigid spheres in which mean porosity values are associated with the processes of formation (Dullien, 2012). The range of porosities covers the range from, thinnest regular packing 0.48, very loose random packing 0.44, loose random packing 0.40–0.41, poured random packing 0.375–0.391 to minimum porosities in close random packing 0.36–0.37 and 0.26 in densest regular packing (hexagonal close pack) (Dullien, 2012). Porosity in hepatocytes spheroid can be lower than 0.26 since the cells are not rigid and can adhere together firmly (Khakpour et al., 2017). Hence, the presented parametric study is performed for three different spheroid porosity, ε_sph, namely 0.2, 0.3, 0.5, which ε_sph = 0.2 is considered as the baseline value. Furthermore, since in larger spheroid diameter, cells especially in the core sections of the aggregates may suffer from severe oxygen limitation, oxygen concentration inside the different spheroid diameters is evaluated, and feasibility of reaching physiological oxygen concentration is analyzed. The parametric study was carried out for the last spheroid, which is considered as the severe case of study since its location is in the furthest distance from the inlet.

TABLE 4

Table 4. Parameters, their respected default value and their range over which they were investigated in the parametric study.

Dimensional Analysis

Dimensional analysis is carried out to investigate the dominance of parameters participating in the transport mechanism. The parameters are defined in Table 5, and their values are given in Table 6. The Peclet number compares the rate of convection to the diffusion for oxygen. In the present study, three Peclet numbers are defined: Pe_sph, Peclet number for the intra spheroid section, Pe_hyd, Peclet number for hydrogel part and Pe_main, Peclet number for the main-stream section. The oxygen diffusion coefficient is acquired by Equation (9) based on the baseline porosity value (3 × 10⁻¹⁰ m²/s for ε = 0.2), which increases with an increase in porosity. Based on the low-velocity rate inside the spheroids and spheroid diameter of 200 μm (baseline values), the Peclet number is generally in a range from O(10⁻²) to O(10⁻⁶) obtained based on the reported range of perfusion rate and aggregate diameter (Table 4). Such low values of the Peclet number indicates that the transport mechanism in intra spheroid space is mainly diffusion dominant.

TABLE 5

Table 5. Dimensional parameters used for evaluation of oxygen concentration distribution inside the bioreactor.

TABLE 6

Table 6. Values obtained for dimensionless parameters considering baseline values in Table 4.

The Peclet number for hydrogel section is computed mostly around O(1) by considering the aforementioned range for perfusion rate, element size 2 mm and the oxygen diffusion coefficient in hydrogel over 1.6 × 10⁻⁹ m²/s. The O(1) is denoting the comparable effects of diffusion and convection rate. Additionally, the Peclet number [O(100)] is also estimated for the main-stream section by considering 4 × 10⁻³ m as element size and 3.4 × 10⁻⁹ m²/s for the oxygen diffusion coefficient, showing the dominance of convection transport mechanism. Higher perfusion rate will decrease the oxygen concentration gradient inside the hydrogel section, which in turn results in the higher level of oxygen supply at the vicinity of the hepatic spheroids (Figure S2). The order of magnitude of Peclet numbers is O(10⁻³), O(1), and O(10²) for intra spheroid, hydrogel, and main-stream, respectively, representing the dominance of diffusion rate inside the spheroid and prevalence of the convection mechanism in the main-stream.

The relation between the oxygen reaction rate and diffusion rate is investigated through the Damkohler number, which compares the rate of oxygen diffusion rate from hydrogel to spheroid to the maximum reaction rate V_m. By considering the aforesaid parameters, Da number varies from 3.2 to 5.5.

Besides, two more indexes are introduced here to show the conditions in which bioreactor provide a suitable microenvironment for the hepatocytes. Biomimicry index extracted from a similar method of previous studies carried out by Dulong and Legallais (2007) and Khakpour et al. (2017), indicating the ratio of regions in which hepatocytes are exposed to physiological level of oxygen supply (4–9%) (Khakpour et al., 2017), though providing a higher or lower level of oxygen supplies does not necessarily indicate total loss of hepatocyte functionalities (Tourlomousis and Chang, 2016a). Hypoxia index also illustrates the region in the spheroids where oxygen levels fall below 4% (42 μmol/L), which can be helpful in choosing appropriate design and working parameters of the bioreactor.

Next, the effect of changing various operational and design parameters affecting the efficiency of the bioreactor will be analyzed. Operational parameters like inlet oxygen concentration and perfusion rate can easily be modified to provide an optimized oxygen supply for the hepatocytes. Changing design parameters such as spheroid diameter, hydrogel porosity, and kinetic parameters is not straightforward since hepatocytes diameter is changing over culture time. Also, it is possible to choose different biocompatible hydrogel for encapsulating the cells, but hydrogel porosity ranges are not large (between 0.3 and 0.7) (Tan et al., 2013; Mai et al., 2014; Pacelli et al., 2018). Kinetic parameters like (V_m, K_m) are intrinsic to the cell type, and therefore, these parameters cannot be altered (Curcio et al., 2014).

Oxygen concentration at the inlet can be modified by modifying the oxygen saturation level in the perfusion media via changing the oxygen partial pressure of the incubator (Equation 14). Figure 9A displays the effect of changing oxygen partial pressure inside the incubator where intra spheroid oxygen level follows the similar increment of oxygen partial pressure, especially at higher porosities. Higher oxygen pressure will increase intra spheroid oxygen level. Hence in terms of dimensionless parameters, an increase of oxygen partial pressure from 100 to 300 μmol/L will result in elevated intra spheroid oxygen concentration with the almost similar linear trend. Additionally, supplying hepatocytes with higher oxygen level will push down the Da. A gradual increase in oxygen concentration from 100 to 300 μmol/L at ε_sph = 0.2 will cause a reduction in both Da, which falls from 8.3 to 3.3.

FIGURE 9

Figure 9. Analyzing the effect of different parameters on spheroid oxygen concentration. Measuring oxygen concentration level inside the spheroid for different (A) Dissolve oxygen concentration (B) Perfusion rate (C) Spheroid diameter (D) Da (E) Hydrogel porosity (F) Intra spheroid porosity.

Perfusion rate, as the second operational parameters studied here, can be easily tuned by modifying the rotational speed of the peristaltic pump. Increasing the perfusion rate is accompanied by higher oxygen delivery to the cells (Figure 9B), but there should be a limit for the maximum value of the perfusion rate based on the shear stress threshold tolerated by hepatocytes. It can be noticed spheroid oxygen concentrations do not change significantly under the studied perfusion rate.

Spheroid diameter is changing, resulting from hepatocytes proliferation during the culture period. Hepatocytes in larger hepatic spheroids are more prone to lower oxygen levels considering the diffusion as the dominant oxygen transport mechanism. Figure 9C depicts the oxygen concentration level for spheroids of different sizes. Oxygen supply level in the smaller spheroids even at low porosities, based on the short diffusive penetration depth of oxygen, is almost in normoxic conditions. Increasing spheroid diameter from 100 to 400 μm will enhance β_Hpx from 0 to 0.64 showing that spheroids with >200 μm, especially at lower porosities are at the high risk of hypoxic conditions.

Measuring K_m value is quite challenging, which requires complex evaluation; hence, data reported for different values of K_m does not cover a large domain. Three values of 0.7, 1, and 7.8 μmol/L are chosen based on the values given by Hay (Dunn et al., 1991) and Khakpour et al. (2017). Figure 9D shows the dependency of hepatocytes oxygen concentration on Da variation for different aggregates porosities. Oxygen concentration level inside the spheroid is highly dependent on V_m values, where an increase from 5 to 20 nmol/(s-cm3) will contribute to around 30 % raise in β_Hpx values at a porosity of 0.2. The variation trend slows down at higher porosities. For instance, ε_sph = 0.5 is accompanied by 10% raise for β_Hpx.

It is also possible to modify the porosity of the hydrogel by using different appropriate hydrogels or applying various conditions of hydrogel cross-linking (Tourlomousis and Chang, 2016b). Modifying the membrane porosity does not affect the perfusion rate of the medium, provided that the peristaltic pump circulates the media at the constant flow rate. However, a hydrogel with high porosity will reduce pressure drop requires for flow perfusion (Equation 13). Hydrogel with higher porosity provides less oxygen transport resistance, and thereby higher level of oxygen supply to the cells (Figure 9E) showing that changing hydrogel porosities will modify both maximum and minimum of oxygen concentration inside the hepatocytes. Increasing hydrogel porosity will cause a slight increment in biomimicry index (12%) at higher porosities and remain almost similar at ε = 0.2.

Effects of different spheroid porosities on the intra spheroid oxygen concentration are also investigated for two K_m values, i.e., 0.7, and 7.8 μmol/L representing the minimum, the maximum values of the oxygen concentration, respectively. Figure 9F represents the effects of changing intra spheroid porosity on the oxygen level inside the spheroid, where the oxygen level is much higher in aggregates with larger porosities, especially in minimum values. Higher spheroid porosity is correlated with a less number of cells in the aggregates with the specific diameter (Equation 8), which in turn leading to a less oxygen demand from the cells inside the spheroid. Furthermore, the effective diffusion coefficient increases with an increase in the porosity of the spheroids (Equations 9–11), resulting in less mass transfer (oxygen delivery) resistance to the cells. Increasing ε_sph from 0.2 to 0.5 resulted in increasing of Pe_sph level (from 2 × 10⁻³ to 7 × 10⁻²) and decreasing of β_Bm (from 0.13 to 0).

Conclusion

Here in the presented study, a 3D computational model of perfusion bioreactor containing hepatocyte spheroids encapsulated in the hydrogel was investigated, and the influence of different parameters affecting the intra spheroid oxygen concentration was inspected. A full hydrodynamic model of medium flow inside the bioreactor confirmed that shear stress imposed on hepatocytes was far below the detrimental threshold (0.03 Pa). Hosting spheroids inside the hydrogel would protect cells from shear stress, and hence, in conditions in which higher oxygen level would be required, such as aggregates with the larger diameters or the higher number of cells, the higher level of flow rates could be applied.

Additionally, the presented model incorporated mass transfer of various metabolites (i.e., oxygen, albumin, glucose, glutamine, ammonia, and urea) known as crucial parameters in liver specific-functions. The distribution of the oxygen concentration inside the bioreactor was obtained. Albumin, glucose, glutamine, ammonia, and urea distribution inside the bioreactor were also evaluated. Predicted results of the albumin production represented a good agreement with the experimental data. Results of glucose and glutamine predictions would be of great importance in generating insight for researchers by providing a baseline index prior to any fabrication. For example, obtained results indicated the low glucose consumption rate of the hepatocytes, would be beneficial in estimating the minimum necessary glucose requirement of the hepatocytes. Ammonia is a toxic compound which is eliminated by hepatocytes from the circulating system through the urea cycle. Ammonia in the presented model had two terms of elimination and generation. It was removed by converting to urea and generated in the glutamine synthesis. Ammonia to urea conversion was modeled through four enzymatic reactions. The detailed concentration of each participating metabolite was obtained by solving 12 coupled differential equations. Estimated results obtained from urea synthesis could be further applied to predict the diseases relating to the urea cycle, such as inborn urea cycle disorders.

Since oxygen tension is a critical parameter in viability and maintaining the liver-specific functionalities of the cells, in the next step of this study, a parametric study was carried out based on the oxygen concentration and the effect of different design conditions on the oxygen level inside the hepatic spheroids was investigated. Oxygen tension inside the hepatic aggregates were obtained for the different perfusion rates. Since the dominant transport mechanism in the main-stream is convection [Pe_main O(100)], higher perfusion rates would provide higher oxygen supply to the vicinity of the spheroids. Among different studied design parameters, saturation oxygen concentration at the inlet was the most influential parameter which could be tuned based on the required level of the oxygen in the system.

Numerical results showed that spheroids with smaller diameters were less prone to hypoxia conditions based on their short diffusive penetration depth and for the spheroids with the diameter larger than 200 μm, the precise control of oxygen delivery must be applied. Transport mechanism inside the spheroids was diffusion dominant [Pe_sph O(10⁻³)]. Maximum oxygen consumption rate, V_m, also played an important role in determining the oxygen distribution inside the spheroids. High V_m values could lead to anoxia condition, which was accompanied by a high level of Da. Spheroid porosity was shown to be an essential factor affecting both the number of the hepatocytes and the oxygen concentration distribution inside the spheroid. It is worth mentioning that hydrogel porosity noticeably affected the intra spheroid oxygen concentration. Both convection and diffusion had a comparable dominance in the hydrogel section [Pe_hyd O(1)]; hence, mass transfer resistance caused by hydrogel porosity was found to be a pivotal parameter affecting the intra spheroid oxygen concentration. This multi-scale computational system allows the study of some main liver-specific functions and optimization of bioreactor characteristics based on efficient oxygen transport. We also envision to develop a more complex model of metabolisms and include the other liver-specific metabolisms to the presented model.

Data Availability

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

FS and BF conceived of the presented idea. FS developed the theory and performed the computations. She developed the theoretical formalism, performed the analytic calculations, and performed the numerical simulations. BF and KF verified the analytical methods and supervised the findings of this work. FS wrote the manuscript under the supervision of BF and KF. All authors discussed the results and contributed to the final manuscript.

Funding

We gratefully acknowledge the financial support (Grant# G940502) from Deputy of Research office of Sharif University of Technology.

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2019.00221/full#supplementary-material

References

Ahmed, H. M. M., Salerno, S., Piscioneri, A., Khakpour, S., Giorno, L., and De Bartolo, L. (2017). Human liver microtissue spheroids in hollow fiber membrane bioreactor. Colloids Surf. B. Biointerf. 160, 272–280. doi: 10.1016/j.colsurfb.2017.09.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Anada, T., Fukuda, J., Sai, Y., and Suzuki, O. (2012). An oxygen-permeable spheroid culture system for the prevention of central hypoxia and necrosis of spheroids. Biomaterials 33, 8430–8441. doi: 10.1016/j.biomaterials.2012.08.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Bao, J., Shi, Y., Sun, H., Yin, X., Yang, R., Li, L., et al. (2011). Construction of a portal implantable functional tissue-engineered liver using perfusion-decellularized matrix and hepatocytes in rats. Cell Transplant. 20, 753–766. doi: 10.3727/096368910X536572

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrett, K. E., Barrett, K., and Barrett, K. E. (2006). Gastrointestinal Physiology. New York, NY: Mcgraw-Hill.

Google Scholar

Baudoin, R., Griscom, L., MatthieuProt, J., Legallais, C., and Leclerc, E. (2011). Behavior of HepG2/C3A cell cultures in a microfluidic bioreactor. Biochem. Eng. J. 53, 172–181. doi: 10.1016/j.bej.2010.10.007

CrossRef Full Text | Google Scholar

Bhise, N. S., Manoharan, V., Massa, S., Tamayol, A., Ghaderi, M., Miscuglio, M., et al. (2016). A liver-on-a-chip platform with bioprinted hepatic spheroids. Biofabrication 8:014101. doi: 10.1088/1758-5090/8/1/014101

PubMed Abstract | CrossRef Full Text | Google Scholar

Bodenberger, N., Kubiczek, D., Abrosimova, I., Scharm, A., Kipper, F., Walther, P., et al. (2016). Evaluation of methods for pore generation and their influence on physio-chemical properties of a protein based hydrogel. Biotechnol. Rep. 12, 6–12. doi: 10.1016/j.btre.2016.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Curcio, E., Piscioneri, A., Morelli, S., Salerno, S., Macchiarini, P., and De Bartolo, L. (2014). Kinetics of oxygen uptake by cells potentially used in a tissue engineered trachea. Biomaterials 35, 6829–6837. doi: 10.1016/j.biomaterials.2014.04.100

PubMed Abstract | CrossRef Full Text | Google Scholar

Curcio, E., Salerno, S., Barbieri, G., De Bartolo, L., Drioli, E., and Bader, A. (2007). Mass transfer and metabolic reactions in hepatocyte spheroids cultured in rotating wall gas-permeable membrane system. Biomaterials. 28, 5487–5497. doi: 10.1016/j.biomaterials.2007.08.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Dullien, F. A. (2012). Porous Media: Fluid Transport and Pore Structure. San Diego, CA: Academic press.

Google Scholar

Dulong, J. L., and Legallais, C. (2007). A theoretical study of oxygen transfer including cell necrosis for the design of a bioartificial pancreas. Biotechnol. Bioeng. 96, 990–998. doi: 10.1002/bit.21140

PubMed Abstract | CrossRef Full Text | Google Scholar

Dumont, C., Piselli, J., Temple, S., Dai, G., and Thompson, D. M. (2018). Endothelial cells exposed to fluid shear stress support diffusion based maturation of adult neural progenitor cells. Cell. Mol. Bioeng. 11, 117–130. doi: 10.1007/s12195-017-0516-5

CrossRef Full Text | Google Scholar

Dunn, J. C., Tompkins, R. G., and Yarmush, M. L. (1991). Long-term in vitro function of adult hepatocytes in a collagen sandwich configuration. Biotechnol. Prog. 7, 237–245. doi: 10.1021/bp00009a007

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebrahimkhani, M. R., Neiman, J. A., Raredon, M. S., Hughes, D. J., and Griffith, L. G. (2014). Bioreactor technologies to support liver function in vitro. Adv. Drug Deliv. Rev. 69, 132–157. doi: 10.1016/j.addr.2014.02.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ekins, S., Murray, G. I., Burke, M. D., Williams, J. A., Marchant, N. C., and Hawksworth, G. M. (1995). Quantitative differences in phase I and II metabolism between rat precision-cut liver slices and isolated hepatocytes. Drug Metabol. Dispos. 23, 1274–1279.

PubMed Abstract | Google Scholar

Esch, M. B., King, T. L., and Shuler, M. L. (2011). The role of body-on-a-chip devices in drug and toxicity studies. Ann. Rev. Biomed. Eng. 13, 55–72. doi: 10.1146/annurev-bioeng-071910-124629

PubMed Abstract | CrossRef Full Text | Google Scholar

Evenou, F., Fujii, T., and Sakai, Y. (2010). Spontaneous formation of highly functional three-dimensional multilayer from human hepatoma Hep G2 cells cultured on an oxygen-permeable polydimethylsiloxane membrane. Tissue Eng. Part C: Methods. 16, 311–318. doi: 10.1089/ten.tec.2009.0042

PubMed Abstract | CrossRef Full Text | Google Scholar

Fournier, R. L. (2017). Basic Transport Phenomena in Biomedical Engineering. Roca Raton, FL: CRC press.

Google Scholar

Gebhardt, R. (1992). Metabolic zonation of the liver: regulation and implications for liver function. Pharmacol. Therapeut. 53, 275–354. doi: 10.1016/0163-7258(92)90055-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Glicklis, R., Merchuk, J. C., and Cohen, S. (2004). Modeling mass transfer in hepatocyte spheroids via cell viability, spheroid size, and hepatocellular functions. Biotechnol. Bioeng. 86, 672–680. doi: 10.1002/bit.20086

PubMed Abstract | CrossRef Full Text | Google Scholar

Godoy, P., Hewitt, N. J., Albrecht, U., Andersen, M. E., Ansari, N., Bhattacharya, S., et al. (2013). Recent advances in 2D and 3D in vitro systems using primary hepatocytes, alternative hepatocyte sources and non-parenchymal liver cells and their use in investigating mechanisms of hepatotoxicity, cell signaling and ADME. Arch. Toxicol. 87, 1315–1530. doi: 10.1007/s00204-013-1078-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Groneberg, D. A., Grosse-Siestrup, C., and Fischer, A. (2002). In vitro models to study hepatotoxicity. Toxicol. Pathol. 30, 394–399. doi: 10.1080/01926230252929972

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, H. W., Lee, W. C., Leong, M. C., Sonam, S., Vedula, S. R. K., Lim, C. T., et al. (2011a). Microfluidics for applications in cell mechanics and mechanobiology. Cell. Mol. Bioeng. 4, 591–602. doi: 10.1007/s12195-011-0209-4

CrossRef Full Text | Google Scholar

Hou, Y. T., Ijima, H., Matsumoto, S., Kubo, T., Takei, T., Sakai, S., et al. (2010). Effect of a hepatocyte growth factor/heparin-immobilized collagen system on albumin synthesis and spheroid formation by hepatocytes. J. Biosci. Bioeng. 110, 208–216. doi: 10.1016/j.jbiosc.2010.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, Y. T., Ijima, H., Takei, T., and Kawakami, K. (2011b). Growth factor/heparin-immobilized collagen gel system enhances viability of transplanted hepatocytes and induces angiogenesis. J. Biosci. Bioeng. 112, 265–272. doi: 10.1016/j.jbiosc.2011.05.003

PubMed Abstract | CrossRef Full Text | 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

Hsu, S., Chu, J. S., Chen, F. F., Wang, A., and Li, S. (2011). Effects of fluid shear stress on a distinct population of vascular smooth muscle cells. Cell. Mol. Bioeng. 4, 627–636. doi: 10.1007/s12195-011-0205-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Huh, D., Leslie, D. C., Matthews, B. D., Fraser, J. P., Jurek, S., Hamilton, G. A., et al. (2012). A human disease model of drug toxicity–induced pulmonary edema in a lung-on-a-chip microdevice. Sci. Transl. Med. 4:159ra147. doi: 10.1126/scitranslmed.3004249

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutson, S. M., Stinson-Fisher, C., Shiman, R., and Jefferson, L. S. (1987). Regulation of albumin synthesis by hormones and amino acids in primary cultures of rat hepatocytes. Am. J. Physiol. Endocrinol. Metabol. 252, E291–E298. doi: 10.1152/ajpendo.1987.252.3.E291

PubMed Abstract | CrossRef Full Text | Google Scholar

Jungermann, K., and Kietzmann, T. (2000). Oxygen: modulator of metabolic zonation and disease of the liver. Hepatology. 31, 255–260. doi: 10.1002/hep.510310201

PubMed Abstract | CrossRef Full Text | Google Scholar

Kern, A., Bader, A., Pichlmayr, R., and Sewing, K. F. (1997). Drug metabolism in hepatocyte sandwich cultures of rats and humans. Biochem. Pharmacol. 54, 761–772. doi: 10.1016/S0006-2952(97)00204-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Khakpour, S., Di Renzo, A., Curcio, E., Di Maio, F. P., Giorno, L., De Bartolo, L., et al. (2017). Oxygen transport in hollow fibre membrane bioreactors for hepatic 3D cell culture: a parametric study. J. Membrane Sci. 544, 312–322. doi: 10.1016/j.memsci.2017.09.024

CrossRef Full Text | Google Scholar

Khetani, S. R., and Bhatia, S. N. (2008). Microscale culture of human liver cells for drug development. Nat. Biotechnol. 26, 120–126. doi: 10.1038/nbt1361

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, B. S., Park, K. I., Hoshiba, T., Jiang, H. L., Choi, Y. J., Akaike, T., et al. (2011). Design of artificial extracellular matrices for tissue engineering. Prog. Poly. Sci. 36, 238–268. doi: 10.1016/j.progpolymsci.2010.10.001

CrossRef Full Text | Google Scholar

Kuchel, P. W., Roberts, D. V., and Nichol, L. W. (1977). The simulation of the urea cycle: correlation of effects due to inborn errors in the catalytic properties of the enzymes with clinical-biochemical observations. Immunol. Cell Biol. 55:309. doi: 10.1038/icb.1977.26

PubMed Abstract | CrossRef Full Text | Google Scholar

Laemmle, A., Gallagher, R. C., Keogh, A., Stricker, T., Gautschi, M., Nuoffer, J. M., et al. (2016). Frequency and pathophysiology of acute liver failure in ornithine transcarbamylase deficiency (OTCD). PLoS ONE 11:e0153358. doi: 10.1371/journal.pone.0153358

PubMed Abstract | CrossRef Full Text | Google Scholar

LeCluyse, E. L., Audus, K. L., and Hochman, J. H. (1994). Formation of extensive canalicular networks by rat hepatocytes cultured in collagen-sandwich configuration. Am. J. Physiol. Cell Physiol. 266, C1764–C1774. doi: 10.1152/ajpcell.1994.266.6.C1764

PubMed Abstract | CrossRef Full Text | Google Scholar

Longsworth, L. (1953). Diffusion measurements, at 25, of aqueous solutions of amino acids, peptides and sugars. J. Am. Chem. Soc. 75, 5705–5709. doi: 10.1021/ja01118a065

CrossRef Full Text | Google Scholar

Maguire, T. J., Novik, E., Chao, P., Barminko, J., Nahmias, Y., Yarmush, M. L., et al. (2009). Design and application of microfluidic systems for in vitro pharmacokinetic evaluation of drug candidates. Curr. Drug Metabol. 10, 1192–1199. doi: 10.2174/138920009790820093

PubMed Abstract | CrossRef Full Text | Google Scholar

Mai, M., Hempel, U., Hacker, M. C., and Dieter, P. (2014). Effects of HyStem™-HP hydrogel elasticity on osteogenic differentiation of human mesenchymal stromal cells. Cell. Mol. Bioeng. 7, 155–164. doi: 10.1007/s12195-013-0314-7

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

Merdaw, A., Sharif, A., and Derwish, G. (2010). Water permeability in polymeric membranes, Part I. Desalination 260, 180–192. doi: 10.1016/j.desal.2010.04.042

CrossRef Full Text | Google Scholar

Moussy, Y. (2003). Convective flow through a hollow fiber bioartificial liver. Arti. Org. 27, 1041–1049. doi: 10.1046/j.1525-1594.2003.07074.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Munson, B. R., Rothmayer, A. P., Okiishi, T. H., and Huebsch, W. W. (2014). Fundamentals of Fluid Mechanics. Jefferson City, MO: John Wiley & Sons.

Google Scholar

Murphy, K. C., Hung, B. P., Browne-Bourne, S., Zhou, D., Yeung, J., Genetos, D. C., et al. (2017). Measurement of oxygen tension within mesenchymal stem cell spheroids. J. R. Soc. Interf. 14:20160851. doi: 10.1098/rsif.2016.0851

PubMed Abstract | CrossRef Full Text | Google Scholar

Pacelli, S., Basu, S., Berkland, C., Wang, J., and Paul, A. (2018). Design of a cytocompatible hydrogel coating to modulate properties of ceramic-based scaffolds for bone repair. Cell. Mol. Bioeng. 11, 211–217. doi: 10.1007/s12195-018-0521-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, J., Li, Y., Berthiaume, F., Toner, M., Yarmush, M. L., and Tilles, A. W. (2008). Radial flow hepatocyte bioreactor using stacked microfabricated grooved substrates. Biotechnol. Bioeng. 99, 455–467. doi: 10.1002/bit.21572

PubMed Abstract | CrossRef Full Text | Google Scholar

Picioreanu, C. M., van Loosdretch, M. C. M., and Heijnen, J. J. (1997). Modelling the effect of oxygen concentration on nitrite accumulation in a biofilm airlift suspension reactor. Water Sci. Technol. 36:147. doi: 10.2166/wst.1997.0034

CrossRef Full Text | Google Scholar

Ramaiahgari, S. C., den Braver, M. W., Herpers, B., Terpstra, V., Commandeur, J. N., van de Water, B., et al. (2014). A 3D in vitro model of differentiated HepG2 cell spheroids with improved liver-like properties for repeated dose high-throughput toxicity studies. Arch. Toxicol. 88, 1083–1095. doi: 10.1007/s00204-014-1215-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Rothschild, M. A., Oratz, M., and Schreiber, S. S. (1977). “Albumin synthesis,” in Albumin: Structure, Function and Uses, eds V. M. Rosenoer, M. Oratz, and M. A. Rothschild (Oxford, UK: Elsevier, 227–253.

Google Scholar

Shipley, R. J., and Waters, S. L. (2011). Fluid and mass transport modelling to drive the design of cell-packed hollow fibre bioreactors for tissue engineering applications. Mathemat. Med. Biol. 29, 329–359. doi: 10.1093/imammb/dqr025

PubMed Abstract | CrossRef Full Text | Google Scholar

Skardal, A., Smith, L., Bharadwaj, S., Atala, A., Soker, S., and Zhang, Y. (2012). Tissue specific synthetic ECM hydrogels for 3-D in vitro maintenance of hepatocyte function. Biomaterials 33, 4565–4575. doi: 10.1016/j.biomaterials.2012.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Soldatow, V. Y., Lecluyse, E. L., Griffith, L. G., and Rusyn, I. (2013). In vitro models for liver toxicity testing. Toxicol. Res. 2, 23–39. doi: 10.1039/C2TX20051A

PubMed Abstract | CrossRef Full Text | Google Scholar

Stein, W. D., and Litman, T. (2014). Channels, Carriers, and Pumps: an Introduction to Membrane Transport. London Elsevier.

Google Scholar

Tan, G. D., Toh, G. W., Birgersson, E., Robens, J., van Noort, D., and Leo, H. L. (2013). A thin-walled polydimethylsiloxane bioreactor for high-density hepatocyte sandwich culture. Biotechnol. Bioeng. 110, 1663–1673. doi: 10.1002/bit.24822

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanaka, Y., Yamato, M., Okano, T., Kitomori, T., and Sato, K. (2006). Evaluation of effects of shear stress on hepatocytes by a microchip-based system. Measure. Sci. Technol. 17:3167. doi: 10.1088/0957-0233/17/12/S08

CrossRef Full Text | Google Scholar

Tilles, A. W., Baskaran, H., Roy, P., Yarmush, M. L., and Toner, M. (2001). Effects of oxygenation and flow on the viability and function of rat hepatocytes cocultured in a microchannel flat-plate bioreactor. Biotechnol. Bioeng. 73, 379–389. doi: 10.1002/bit.1071

PubMed Abstract | CrossRef Full Text | Google Scholar

Tourlomousis, F., and Chang, R. C. (2016a). Numerical investigation of dynamic microorgan devices as drug screening platforms. Part I: macroscale modeling approach & validation. Biotechnol. Bioeng. 113, 612–622. doi: 10.1002/bit.25822

PubMed Abstract | CrossRef Full Text | Google Scholar

Tourlomousis, F., and Chang, R. C. (2016b). Numerical investigation of dynamic microorgan devices as drug screening platforms. Part II: microscale modeling approach and validation. Biotechnol. Bioeng. 113, 623–634. doi: 10.1002/bit.25824

PubMed Abstract | CrossRef Full Text | Google Scholar

Truskey, G. A., Yuan, F., and Katz, D. F. (2004). Transport Phenomena in Biological Systems. Upper Saddle River, NJ: Pearson Prentice Hall.

Google Scholar

Vafai, K. (2015). Handbook of Porous Media. New York, NY: CRC Press.

Google Scholar

van Midwoud, P. M., Merema, M. T., Verpoorte, E., and Groothuis, G. M. (2010). A microfluidic approach for in vitro assessment of interorgan interactions in drug metabolism using intestinal and liver slices. Lab Chip. 10, 2778–2786. doi: 10.1039/c0lc00043d

PubMed Abstract | CrossRef Full Text | Google Scholar

Vogel, S. (1988). Life's Devices: the Physical World of Animals and Plants. Princeton, NJ: Princeton University Press.

Google Scholar

Wang, Y., Susando, T., Lei, X., Anene-Nzelu, C., Zhou, H., Liang, L. H., et al. (2010). Current development of bioreactors for extracorporeal bioartificial liver. Biointerphases. 5, FA116–FA131. doi: 10.1116/1.3521520

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, L. J., Matias, J., Scudiero, D. A., Hite, K. M., Monks, A., Sausville, E. A., et al. (2001). P450 enzyme expression patterns in the NCI human tumor cell line panel. Drug Metabol. Dispos. 29, 304–312.

PubMed Abstract | Google Scholar

Yue, K., Trujillo-de Santiago, G., Alvarez, M. M., Tamayol, A., Annabi, N., and Khademhosseini, A. (2015). Synthesis, properties, and biomedical applications of gelatin methacryloyl (GelMA) hydrogels. Biomaterials 73, 254–271. doi: 10.1016/j.biomaterials.2015.08.045

PubMed Abstract | CrossRef Full Text | Google Scholar