Oxygen Consumption Characteristics in 3D Constructs Depend on Cell Density

Introduction

Engineered tissues have many applications, ranging from in-vitro models of human pathophysiology to platforms for drugs and treatment testing (Lee et al., 2005), to providing an alternative to the shortage of human donor organs for transplants (Mattei et al., 2017b, 2018). However, the supply of oxygen and other nutrients as well as the removal of toxic waste products, which typically occur via passive diffusion in in-vitro cellular constructs, often limit cell viability in thick (millimeter-sized) engineered tissues (Mattei et al., 2014). The typical diffusion limit of cell-rich tissues, such as skeletal muscle or liver, is considered to be of ~200 μm, which defines the side of the so-called functional unit, i.e., the smallest autonomous organ cubic voxel which can survive in the absence of blood vessels (Helmlinger et al., 1997; Mattei et al., 2014). In-vivo, nutrient supply and waste removal for thicker tissues are enhanced by convection through the vascular system, which is typically absent in-vitro (Folkman and Moscona, 1978; Rivron et al., 2008). Oxygen is considered the limiting factor when culturing three-dimensional (3D) cell constructs in-vitro, especially in the case of non-porous systems where its mass transport relies only on gradient-driven passive diffusion. The reason for this is its poor solubility in culture media (typically ~0.2 mM, when atmospheric oxygen is used), making it difficult to provide a sufficient resource supply to cells in the core of the engineered 3D constructs. In fact, although oxygen is typically consumed at a similar molar rate per cell as glucose and has a ~4-fold higher diffusion coefficient in aqueous media, this is more than offset by its very low solubility with respect to that of glucose, which results in an O₂ availability of about 0.05–0.2 mM under physiologically relevant conditions, vs. 3–15 mM of glucose (Martin and Vermette, 2005). Moreover, using pure oxygen instead of air or increasing gas pressure without an appropriate carrier such as hemoglobin to increase oxygen concentration in the culture medium has been shown to induce the presence of free radicals, which are cytotoxic (Freshney, 2015). For this reason, the culture medium is often re-circulated and re-oxygenated by passing through an in-line gas exchanger (Mattei et al., 2017a; Ahluwalia et al., 2018). As oxygen is not only crucial for cell survival but also a determinant for cell fate and function, a better understanding of its diffusion and consumption in in-vitro constructs is critical to enable improvements in tissue engineering.

In most studies on oxygen utilization, the parameter reported is the average cellular oxygen consumption rate or OCR, expressed as [mol·cell⁻¹ · s⁻¹], typically estimated by measuring the overall oxygen consumption rate (in mol·s⁻¹) in the system and dividing by the total number of cells present. The range of OCR is usually between 1 × 10⁻¹⁶ and 1 × 10⁻¹⁸ mol · cell⁻¹·s⁻¹ (or 100 to 1 amol·cell⁻¹ · s⁻¹) (Wagner et al., 2011), but the fact that it is an average consumption rate per cell and not an absolute value is often overlooked. Noticeably, the value varies according to experimental conditions and cell type; for instance, the OCR measured for a 2D monolayer of hepatocytes (Smith et al., 1996; Balis et al., 1999)—which are known to have oxygen consumption rates of about 10–100 times higher than most other cell types (Cho et al., 2007)—appears to be an order of magnitude greater than that measured for 3D cultures (Shatford et al., 1992; Nyberg et al., 1993, 1994; Sielaff et al., 1997). Moreover, hepatocyte OCR seems to decrease with increasing cell density in 3D cell cultures. Patzer attributed this phenomenon to the fact that denser cultures better approximate the physiological density of native liver, thus cells are less metabolically “stressed” (Patzer, 2004). Interestingly, cells in-vivo generally have much lower OCRs than their in-vitro counterparts (Glazier, 2015; Ahluwalia, 2017; Magliaro et al., 2019).

Computational fluid dynamic (CFD) models have been largely used to investigate whether cell culture parameters (such as culture medium flow) and configuration (e.g., volume and geometry of both fluidic compartment and cell construct) are appropriate for a given cell type. The models enable quantitative estimates of nutrient supply and waste removal capacity as well as flow-induced shear stress. Several mathematical models which describe oxygen consumption and transport in engineered tissues have been published and the resulting spatial gradients of oxygen concentration have been correlated with either cell density (Lewis et al., 2005; Demol et al., 2011), viability (Radisic et al., 2006; Cheema et al., 2012), organization (Malda et al., 2004; Brown et al., 2007), metabolism (Zhou et al., 2008) or growth factor expression (Mac Gabhann et al., 2007; Cheema et al., 2008).

In most of these models, the volumetric oxygen consumption rate within the cell construct (R in [mol·m⁻³·s⁻¹]) is described by the well-known Michaelis-Menten law (Equation 1):

$\begin{array}{ll}R\left(x,y,z,t\right)=\frac{V_{max}·C(x,y,z,t)}{K_m+C(x,y,z,t)}& (1)\end{array}$

Note that R depends on the local oxygen concentration (C) that is generally a function of both space (lateral dimensions x and y, height z) and time (t), i.e., C(x, y, z, t). V_max and K_m are respectively, the maximum volumetric oxygen consumption rate and the Michaelis-Menten constant (corresponding to the oxygen concentration at which consumption drops to half of its maximum and expressed in [mol·m⁻³]). The terms V_max and K_m are considered as constants for a given system. Since C is the only variable, the cellular oxygen consumption rate is “adaptive” only with respect to the local oxygen concentration.

Using the above equation, Ahluwalia (2017) showed that the OCR is identical for each cell plated in a monolayer as they are all the same distance from the oxygen supply and hence they perceive the same oxygen concentration at any instant (i.e., C is the same for all x and y, while z is constant). Conversely, in 3D constructs as in-vivo, cells close to the media or to a capillary perceive high levels of oxygen, and so consume at a higher rate, while those further away perceive lower levels and hence consume at a lower rate. Thus, the average oxygen consumption rate per cell is lower in 3D than in 2D.

Besides cell type and 2 or 3D configuration, the oxygen consumption characteristics may also depend on other factors such as the cell density or even variations in intrinsic consumptions constants such as K_m. In this study we therefore investigated oxygen utilization in 3D constructs as a function of cell density (ρ_cell in [cells·m⁻³]) through a combined experimental-computational approach. In particular, we focused on 3D hydrogels laden with hepatocytes at different densities. Parameters typically considered as constants in mathematical models of oxygen transport (i.e., diffusion coefficient, D) and consumption (i.e., V_max and K_m) were treated as variables in this study and derived by fitting experimental oxygen measurements with a time-dependent computational model.

Materials and Methods

Oxygen Sensing Setup

Experiments were performed in a commercial polydimethylsiloxane (PDMS) bioreactor (LiveBox1, IVTech s.r.l., Italy) featured with a glass top and bottom, which provides an optically transparent cylindrical cell culture chamber (15 mm diameter−16 mm height). A 15 mm diameter RedEye^® Fospor oxygen sensor patch (Ocean Optics, Germany) was attached to the inner side of bioreactor glass bottom to continuously monitor the oxygen concentration at the base of the cell culture chamber. This sensor uses dynamic fluorescence quenching of a ruthenium complex in a sol-gel and allows for non-invasive/non-destructive measurements through a transparent material without consuming oxygen (measurement range: 0–21% O₂; accuracy: ±0.01% O₂; response time: <5 seconds). Real time oxygen measurements were made by reading oxygen sensor fluorescence with a RE-BIFBORO-2 optical fiber (Ocean Optics) coupled to the outside of the bioreactor glass bottom and connected to a NEOFOX-GT phase fluorometer. Partial oxygen pressure was obtained from fluorescence through the Stern-Volmer equation and then converted to a concentration using Henry's Law (Curcio et al., 2010). A two-point calibration curve was obtained first by reading the fluorescence after pipetting 500 μL of fresh Eagles Minimal Essential Medium (EMEM) culture medium into the bioreactor culture chamber. This corresponds to 20% oxygen partial pressure or a concentration of 0.2 mol·m⁻³ (Curcio et al., 2010). The second point was obtained by adding 500 μL of EMEM containing 1% w/v sodium bisulfite (NaHSO3, Sigma-Aldrich), which establishes the point at 0% O₂ partial pressure, i.e., 0 mol·m⁻³ oxygen concentration (Tremper et al., 1986).

Cell Source

Human hepatoma HepG2 cells were purchased from ICLC (Genova, Italy) and cultured according to the supplier's instructions. Cells were grown in EMEM with Earle's Balanced Salts (EBSS) supplemented with 10% v/v fetal bovine serum, 2 mM glutamine and 1% v/v non-essential amino acids in a 37°C/5% CO₂ humidified incubator. A trypsin/EDTA solution was used to detach confluent cultures. The number of viable cells was assessed by manual cell count in a Burker chamber with the trypan blue exclusion test. After the count, cells were centrifuged (300 × g, 5 min, room temperature) and resuspended in media at the required cell density.

Cell-Laden Hydrogel Construct

Hydrogels laden with HepG2 cells at different densities (i.e., 0.5·10¹², 1·10¹², and 5·10¹² cells·m⁻³) were prepared using a sterile 3 mg/mL collagen solution from bovine skin (C4243, Sigma-Aldrich). The collagen solution was mixed in an 8:1 volume ratio with M199 10x medium (Sigma) containing HepG2 cells at 9 times the final desired density. A 1N NaOH solution was added drop-wise to adjust the pre-gel solution pH to 7.4. The preparation was carried out on ice to prevent collagen gelation. Then, 220 μL of the cell-containing pre-gel solution were pipetted onto the bioreactor glass bottom (previously covered with the oxygen sensor patch) and incubated for 1 h at 37°C to allow the collagen polymerization, obtaining a 1.25 mm thick cell-laden collagen construct. All samples were tested immediately after polymerization in order to guarantee a nearly constant initial value of ambient oxygen concentration (20%) throughout their thickness. A sample of the gels was also stained with Live-Dead^®Viability/Cytotoxicity Kit, for mammalian cells (L3224 - ThermoFisher, Waltham, MA USA). The solution was added to sample media according to manufacturer's instructions and incubated for 20 min at room temperature. Gels were then transferred onto glass microscope slides and imaged with an Olympus IX81 inverted microscope with a 10x objective. Figure 1 shows the cells within the hydrogel scaffolds at the different cell densities investigated.

FIGURE 1

Figure 1. HepG2 cells encapsulated in the gels, stained with Live-Dead. Cell densities: (A) 0.5·10¹², (B) 1·10¹², and (C) 5·10¹² cells·m⁻³. Scale bar: 50 μm.

Oxygen Consumption Measurements

Oxygen measurements were performed by placing the setup inside a 37°C/5% CO₂ humidified incubator. At the beginning of the experiment (t = 0), hydrogel constructs polymerised and maintained in the 20% O₂ incubator (thus containing 20% O₂ concentration, in equilibrium with the overlying culture medium) were exposed to a step change in the O₂ concentration of their overlying culture medium. Hypoxic culture conditions were obtained by replacing the 500 μL of EMEM culture medium with the EMEM containing 1% w/v sodium bisulfite. Even though the medium was the same as that used for sensor calibration, we found that—in presence of the hydrogel—the oxygen concentration in the overlying medium was maintained at a constant level of 1·10⁻² mol·m⁻³, instead of the expected 0 mol·m⁻³, during the whole experiment (i.e., up to 3 h, measured at different times and positions with a commercial Ocean Optics needle oxygen sensor; data not shown). Experiments at the three different cell densities were performed in triplicate while keeping all the other experimental parameters (e.g., oxygen sensing setup, hydrogel composition, volume and thickness, culture medium composition, acquisition time, initial and boundary conditions) unaltered. A total of 9 collagen hydrogels were analyzed.

Mathematical Model

Geometry

Figure 2 shows the geometry of the mathematical model describing the oxygen transport throughout the cell-laden construct detailed in section Materials and Methods. The oxygen concentration within the collagen gel of thickness H = 1.25 mm and radius of 7.5 mm is a function of both the vertical coordinate and time, C(z, t). The bulk environment (i.e., culture medium) directly above the collagen gel is assumed to maintain a uniform and time-invariant oxygen concentration, C(H, t) = C_b, representing the hypoxic boundary.

FIGURE 2

Figure 2. Schematic representation of the modeled system consisting of a cell-laden collagen hydrogel polymerised inside the bioreactor chamber on top of the O₂ sensor. The gel is covered with culture medium with a uniform and time-invariant oxygen concentration (C_b). The oxygen sensor enables real-time, non-invasive measurement of O₂ concentration at z = 0 by reading fluorescence via an external optical fiber.

Oxygen Diffusion and Consumption Rate

Given the cylindrical symmetry of the system investigated, mass transfer throughout the cell-laden hydrogel occurs by passive diffusion in the axial direction only, since oxygen concentration is assumed to be independent of the radial and polar dimensions. The mass balance of oxygen in the cell-laden hydrogel is then described by the Fick's second law with a reaction term:

$\begin{array}{ll}\frac{\partial C(z,t)}{\partial t}=D\frac{{\partial }^{2}C(z,t)}{{\partial z}^2}-R(z,t)& (2)\end{array}$

where C [mol^.m⁻³] is the concentration of oxygen as a function of space and time, D [m²·s⁻¹] is the diffusion coefficient of oxygen in the gel (assumed to be spatially homogeneous), and R [mol O₂ consumed·m⁻³·s⁻¹] is the cellular oxygen consumption rate per unit volume of the construct, assumed to be homogeneously populated with cells. Changes in cell number due to proliferation or death during the time scale of the cellular experiments (<1 h) were neglected due to the short duration of the experiment. Cellular oxygen consumption was modeled using the Michaelis-Menten kinetics (Mattei et al., 2014) (Equation 1), by substituting _{Vmax = ρcell}·sOCR and considering only oxygen concentration variations in time and height, i.e., C(z, t), obtaining:

$\begin{array}{ll}\frac{\partial C(z,t)}{\partial t}=D\frac{{\partial }^{2}C(z,t)}{{\partial z}^2}-\frac{{\rho }_{cell}·sOCR·C(z,t)}{K_m+C(z,t)}& (3)\end{array}$

In Equation (3), the sOCR (mol·cell⁻¹·s⁻¹) is the maximum rate of oxygen consumption per single cell. At high oxygen concentrations (i.e., when C(z, t) ≫ K_m), the system is saturated and the rate of reaction is maximal, R(z, t) = _{Vmax = ρcell} · sOCR.

As pointed out in the introduction, most studies on oxygen utilization report the average cellular oxygen consumption rate (OCR, expressed as [mol · cell⁻¹·s⁻¹]), a time-dependent parameter which can be expressed as the volumetric average of the local OCR (dependent on both time and space) within the cell construct of volume V (Equation 4):

$\begin{array}{ll}OCR(t)=\frac{1}{V·{\rho }_{cell}}{\displaystyle \iiint }\frac{{\rho }_{cell}·sOCR·C(x,y,z,t)}{K_m+C(x,y,z,t)}\partial x\partial y\partial z& (4)\end{array}$

Given the cylindrical symmetry of the problem investigated in this work, Equation (4) simplifies to:

$\begin{array}{l}OCR\left(t\right)=\frac{\pi R^2}{V·{\rho }_{cell}}{\displaystyle \int \frac{{\rho }_{cell}·sOCR·C\left(z,t\right)}{K_m+C\left(z,t\right)}}\partial z=\\ \frac{1}{H}{\displaystyle \int \frac{sOCR·C\left(z,t\right)}{K_m+C\left(z,t\right)}\partial z}& (5)\end{array}$

Notably, at high oxygen concentrations, where C(z, t) ≫ K_m ,the OCR tends to the sOCR.

Boundary and Initial Conditions

In each experiment, oxygen concentration measurements started when the EMEM culture medium overlying the hydrogel construct was replaced with that containing 1% w/v sodium bisulfite, setting the experimental time zero (t = 0). Atmospheric air contains 20.9% oxygen, which reduces to 20% O₂ in cell culture incubators due to the presence of 5% CO₂. Because all the hydrogels were prepared under ambient conditions and then polymerised and maintained in a 20% O₂ incubator before experiments, the initial concentration of oxygen in the gel was assumed to be spatially uniform and equal to C(z, t = 0) = C₀ = 0.2 mol·m⁻³. Conversely, the O₂ concentration in the overlying culture medium was experimentally found to be C_b = 1% O₂ (section Oxygen Consumption Measurements), thus giving the boundary condition C(z = H, t) = 1·10⁻² mol·m⁻³. No flux boundary conditions were used for all boundaries except the top surface of the gel (z = H, which was held constant at C_b), implying that:

$\begin{array}{ll}\frac{\partial C(z=0,t)}{\partial z}=0& (6)\end{array}$

Non-dimensionalisation of the Governing Equation

The governing diffusion-reaction equation (Equation 2) was rendered non-dimensional by introducing the following group of variables:

$\begin{array}{lll}\Theta & =\frac{C}{C_0}& (7)\end{array}$

$\begin{array}{lll}\tau & =\frac{t·D}{H^2}& (8)\end{array}$

$\begin{array}{lll}\zeta & =\frac{z}{H}& (9)\end{array}$

where Θ, τ and ζ represent dimensionless concentration, time, and length, respectively.

In case of combined diffusion and reaction phenomena it is useful to introduce the Thiele modulus, a dimensionless parameter representing the ratio between the characteristic reaction rate and diffusion rate of a given species throughout the hydrogel construct, defined as (Cheng et al., 2009):

$\begin{array}{ll}\varphi =H\sqrt{\frac{sOCR·{\rho }_{cell}}{D·K_m}}& (10)\end{array}$

Equation (3) can then be rewritten in the following dimensionless form:

$\begin{array}{ll}\frac{\partial \Theta }{\partial \tau }=\frac{{\partial }^2\Theta }{{\partial \zeta }^2}-{\varphi }^2\frac{\Theta }{1+\frac{\Theta }{K_m}}& (11)\end{array}$

which is subjected to the initial and boundary conditions below:

$\begin{array}{ll}\Theta (0,\zeta )=1& (12)\end{array}$

$\begin{array}{ll}\Theta (\tau ,1)=\frac{C_b}{C_0}& (13)\end{array}$

$\begin{array}{ll}\frac{\partial \Theta (\tau ,0)}{\partial \zeta }=0& (14)\end{array}$

A summary of all the model parameters is presented in Table 1.

TABLE 1

Table 1. Summary of model parameters.

Parameter Estimation

The dynamic process modeling software gPROMS (generalized PRocess Modeling System; PSE; London, UK) was used to fit the oxygen concentrations measured at the base of the hydrogel to the non-dimensional model equations shown in section Non-dimensionalisation of the Governing Equation. The unknown parameters D, V_max and K_m were estimated by fitting experimental data obtained from cell-laden hydrogels at the three different volumetric cell densities (i.e., ρ₁ = 0.5 · 10¹² cells·m⁻³, ρ₂ = 1 · 10¹² cells·m⁻³ and ρ₋₃ = 5 · 10¹² cells·m⁻³). The goodness of fit was assessed by computing the coefficient of variation (CV) for each estimated parameter, defined as the ratio between the standard deviation and the mean of the estimate. In particular, estimates with low values of CV denote a good fit.

Results

Experimental Oxygen Concentration Profiles, Fitting Results, and Estimated Parameters

Figure 3 shows the experimental non-dimensional oxygen concentrations measured at the base of the 3 different cell-laden hydrogel constructs over time, i.e., Θ(z = 0, t). The computed Θ(z = 0, t) temporal profiles—derived by fitting experimental data obtained at different cell densities to the model described in section Mathematical Model—are shown in the same figure along with their respective measured profiles. A good fit was obtained for all the 3 different cell-laden constructs investigated. Estimated parameters and the corresponding CV values are summarized in Table 2. The parameters identified in this study lie within the values reported in the literature (Patzer, 2004; Wagner et al., 2011).

FIGURE 3

Figure 3. Experimental non-dimensional oxygen concentrations vs. time (black lines, average data from n = 3 independent experiments per each cell density investigated) and respective model estimations (dashed red lines) obtained for different cell densities: (A) 0.5·10¹², (B) 1·10¹², and (C) 5·10¹² cells·m⁻³. Error bars—shown every 50 s for the sake of clarity—denote standard deviations of experimental data.

TABLE 2

Table 2. Estimated oxygen diffusion and consumption parameters as a function of cell density.

Notably, the CV values were very low for all the 3 estimated parameters (i.e., D, K_m, and V_max), confirming the goodness of fit and indicating that the Michaelis-Menten law accurately describes the oxygen consumption kinetics in these constructs for each cell density. However, the reaction parameters are not constants but depend on the cell density.

Average vs. Single Cell Oxygen Consumption Rate (OCR vs. sOCR)

After estimating the values of D, K_m, and V_max for each condition, the average OCR of the 3D constructs at different timepoints was obtained using numerical methods through Equation 5. Figure 4 shows the average OCR as a function of cell density, computed at different experimental times (t = 0, 50, 100, 250, 500, 1,000, and 2,000 s).

FIGURE 4

Figure 4. Average cellular oxygen consumption rate (OCR) in hepatocyte laden gels as a function of time for different cell densities. Empty symbols on the y-axis represent the respective (time-independent) sOCRs estimated for different cell densities. The dash-dotted horizontal lines serve only as a guide to the eye for identifying the three sOCR values.

The OCR in the constructs decreases with time and depends on the cell density, as expected from the decrease in oxygen concentration over time and the decrease in sOCR with increasing cell density, respectively. Notably, at time zero, when the oxygen concentration throughout the construct is maximum, the OCR approaches (but does not equal) the sOCR.

Discussion

In this paper we describe a study to accurately evaluate the OCR in a 3D construct without assuming any prior knowledge of the model parameters (D, V_max, K_m). Ehsan and George (2013) report an experimental set up and theoretical approach similar to ours. The authors assume literature values for all the constants and then directly solve for oxygen profiles using finite element software. Here, as far as we know, for the first time the Michaelis- Menten constants are treated as unknowns and an inverse (iterative) problem approach is used for parameter estimation.

Only a small decrease in the oxygen diffusion coefficient throughout the hydrogel construct (D) was observed with increasing cell density, as expected due to the very low volume fraction occupied by cells (<1%). Notably, the expected increase in V_max with increasing cell density was concomitant with a decrease in sOCR and an increase in K_m, both of which suggest the presence of cooperative behavior in the cell oxygen consumption characteristics, which are dependent on cell density. This cooperative behavior is not accounted for in the relatively simple Michaelis-Menten model, which describes cell oxygen consumption kinetics as an adaptive behavior depending only on the local O₂ concentration perceived by cells. All the other parameters (i.e., K_m and sOCR) are considered as constants depending on the specific cell type, but not on their density. The results indicate that cells somehow sense the increase in number and consequently not only reduce their maximal oxygen consumption rate (sOCR), but also their affinity toward oxygen (i.e., decrease the efficiency with which oxygen is captured), as indicated by the increase in K_m. This is also reflected in the fact that the average OCR decreases with increasing cell density as predicted by Ahluwalia (2017). In addition, for a given cell density the (average) OCR decreases with time. If we assume that the cell number is constant (i.e., no cell proliferation or death—reasonable in the short time frame of the 1 h experiments), then clearly the OCR is also a function of oxygen concentration. This may explain why measured OCRs vary so widely in the literature: slight changes in media heights, even in monolayer cultures, can cause large changes in the oxygen concentration perceived by cells.

Beyond these considerations, one should note that in a typical cell culture experiment, oxygen transport and consumption is governed by two Fick's laws (Mattei et al., 2014): (i) one for the “culture medium” domain (Equation 15, subscript m) where oxygen is transported by diffusion (plus convection in the presence of a velocity field u as in bioreactors, or micro-fluidic devices), and (ii) one for the “cell construct” compartment (Equation 16, subscript c), where oxygen is generally transported by passive diffusion and consumed at a rate R, which is typically modeled according to the Michaelis-Menten kinetics, as outlined in the introduction. The internal boundary between culture medium and cell construct domain is characterized by continuity of oxygen concentration and flux.

$\begin{array}{lll}\frac{\partial c_m}{\partial t}& =\nabla ·(D_m\nabla c_m)-u·\nabla c_m& (15)\end{array}$

$\begin{array}{lll}\frac{\partial c_c}{\partial t}& =\nabla ·(D_c\nabla c_c)-R& (16)\end{array}$

Often, the average OCR is derived by measuring the first time-derivative of the oxygen concentration decrease in the “culture medium” compartment (i.e., ∂c_m/∂t) and dividing it by the cell density (or number of cells divided by the volume of media) (Patzer, 2004; Russell et al., 2017), neglecting oxygen gradients in the culture medium and assuming that the rate at which oxygen disappears from the medium is equal to the rate of oxygen utilization by the cell construct.

The assumption that $\frac{\partial c_m}{\partial t}=R$ can be reasonable for 2D cultures, since de facto the cell monolayer represents a boundary of the culture medium compartment, thus the quantity of oxygen disappearing from the media over time is equal to that consumed by cells. However, the same cannot be said in the presence of 3D cultures, where oxygen diffusion throughout the construct cannot be neglected. For a construct in contact with a well-defined medium domain, the rate of decrease in oxygen concentration as measured from the culture medium side is a consequence of both its diffusion [i.e., ∇ · (D_c∇c_c)] and consumption (R) in the 3D cellular domain. Hence, the average OCR cannot be directly derived as ∂c_m/∂t divided by ρ_cell as doing so would likely lead to an overestimation of the cellular uptake rate. Approaches in which the boundary conditions are controllable and known and in which the oxygen consumption is measured within the construct, as described herein are thus required to get meaningful values of average OCR and sOCR.

Since the oxygen consumption kinetics depends (among others) on the available O₂ concentration, it is of interest to consider the following definition of the Thiele modulus (Equation 17), which conveniently includes the transition from a zero- to a first-order consumption kinetics as the O₂ concentration decreases and returns a time-dependent parameter, differently than the classical Thiele modulus which is a constant (Equation 10) (Fink et al., 1973):

$\begin{array}{ll}{\varphi }^{*}(t)=H\sqrt{\frac{sOCR·{\rho }_{cell}}{D·(K_m+C_{av}(t))}}& (17)\end{array}$

In the latter equation, C_av(t) denotes the time-dependent volumetric average of the local oxygen concentration C(z, t), defined as follows for the cylindrically symmetric problem investigated in this work:

$\begin{array}{ll}{C}_{av}(t)=\frac{1}{H}{\displaystyle \int }C(z,t)\partial z& (18)\end{array}$

Figure 5 shows the ϕ^* temporal profiles computed at different experimental times (t = 0, 50, 100, 250, 500, 1,000, and 2,000 s) for each of the three cell densities under study. At each time point, the ϕ^* increases with cell density, as one would expect because of the increasing oxygen demand. At time zero, the oxygen supply meets the cell demand for each of the cell densities investigated, as indicated by ϕ^* ≤ 1. Afterwards, this parameter increases monotonically over time. Notably, the instant at which ϕ^* = 1 (corresponding to the point where oxygen consumption rate equates the oxygen diffusion rate) is reached more rapidly in denser cell-laden hydrogel constructs. Beyond this point, oxygen diffusion becomes the limiting factor in the diffusion-reaction processes occurring within the cellularized hydrogel constructs (Mattei et al., 2014).

FIGURE 5

Figure 5. The time dependent Thiele modulus ϕ* obtained at different experimental times for each of the three cell densities investigated. The dotted horizontal line serves only as a guide to the eye for identifying the instant at which ϕ* = 1.

Conclusions

The oxygen concentration perceived by cells conditions their function and behavior. Therefore, the capacity to monitor oxygen in a culture and to quantify the cell oxygen consumption kinetics is critical for the development of tissue engineered products. However, the measurement of the average cellular consumption rate in 3D constructs, known as the OCR, is difficult because of the lack of tools for the precise quantification of oxygen concentration in space and time in tissues. Our results show that the OCR is neither a constant for a given cell type nor is it constant in a given tissue or tissue construct. Moreover, for given construct dimensions, the Michaelis-Menten constant and the value of V_max as well as the oxygen diffusion coefficient depend on the number of cells in the tissue or scaffold volume. Thus, consumption parameters should not be assumed as literature constants. Given the importance of oxygen consumption in tissues, this study paves the way for a more meaningful interpretation of the significance of OCR, the refinement of analytical techniques and better control of resource supply to 3D in-vitro cell culture systems. The method described in this paper can be used to determine the OCR as a function of cell density in any 3D construct seeded with any cell type. Moreover, it can be also used to calculate the OCR at discrete time points provided the cell density is accurately known at the time points of interest.

Data Availability Statement

The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

Author Contributions

GM and AA conceived and designed the work. FI, AC, and CM performed the experiments and acquired data and images, GM, FI, CM, VP, and AA interpreted and analyzed the data. GM and AA drafted the work. All the authors have approved the submitted version and have agreed both to be personally accountable for their own contributions and to ensure that questions related to the accuracy or integrity of any part of the work, even ones in which the author was not personally involved, are appropriately investigated, resolved, and the resolution documented in the literature.

Funding

CM has received funding from the Fondazione Umberto Veronesi under the Post-Doctoral Fellowship 2019. This work has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 760813 (PATROLS). Funds from the Swiss National Science Foundation, through Project SINERGIA 2019 572 CRSII5_186422/1 are also acknowledged.

Conflict of Interest

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.

References

Ahluwalia, A. (2017). Allometric scaling in-vitro. Sci Rep. 7:42113. doi: 10.1038/srep42113

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahluwalia, A., Misto, A., Vozzi, F., Magliaro, C., Mattei, G., Marescotti, M. C., et al. (2018). Systemic and vascular inflammation in an in-vitro model of central obesity. PLoS ONE 13:e0192824. doi: 10.1371/journal.pone.0192824

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

Brown, D. A., MacLellan, W. R., Laks, H., Dunn, J. C., Wu, B. M., and Beygui, R. E. (2007). Analysis of oxygen transport in a diffusion-limited model of engineered heart tissue. Biotechnol. Bioeng. 97, 962–975. doi: 10.1002/bit.21295

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheema, U., Brown, R. A., Alp, B., and MacRobert, A. J. (2008). Spatially defined oxygen gradients and vascular endothelial growth factor expression in an engineered 3D cell model. Cell Mol. Life Sci. 65, 177–186. doi: 10.1007/s00018-007-7356-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheema, U., Rong, Z., Kirresh, O., MacRobert, A. J., Vadgama, P., and Brown, R. A. (2012). Oxygen diffusion through collagen scaffolds at defined densities: implications for cell survival in tissue models. J. Tissue Eng. Regen. Med. 6, 77–84. doi: 10.1002/term.402

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, G., Markenscoff, P., and Zygourakis, K. (2009). A 3D hybrid model for tissue growth: the interplay between cell population and mass transport dynamics. Biophys. J. 97, 401–414. doi: 10.1016/j.bpj.2009.03.067

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, C. H., Park, J., Nagrath, D., Tilles, A. W., Berthiaume, F., Toner, M., et al. (2007). Oxygen uptake rates and liver-specific functions of hepatocyte and 3T3 fibroblast co-cultures. Biotechnol. Bioeng. 97, 188–199. doi: 10.1002/bit.21225

PubMed Abstract | CrossRef Full Text | Google Scholar

Curcio, E., Macchiarini, P., and De Bartolo, L. (2010). Oxygen mass transfer in a human tissue-engineered trachea. Biomaterials 31, 5131–5136. doi: 10.1016/j.biomaterials.2010.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Demol, J., Lambrechts, D., Geris, L., Schrooten, J., and Van Oosterwyck, H. (2011). Towards a quantitative understanding of oxygen tension and cell density evolution in fibrin hydrogels. Biomaterials 32, 107–118. doi: 10.1016/j.biomaterials.2010.08.093

PubMed Abstract | CrossRef Full Text | Google Scholar

Ehsan, S. M., and George, S. C. (2013). Nonsteady state oxygen transport in engineered tissue: implications for design. Tissue Eng. Part A 19, 1433–1442. doi: 10.1089/ten.tea.2012.0587

PubMed Abstract | CrossRef Full Text | Google Scholar

Fink, D. J., Na, T., and Schultz, J. S. (1973). Effectiveness factor calculations for immobilized enzyme catalysts. Biotechnol. Bioeng. 15, 879–888. doi: 10.1002/bit.260150505

CrossRef Full Text | Google Scholar

Folkman, J., and Moscona, A. (1978). Role of cell shape in growth control. Nature 273, 345–349. doi: 10.1038/273345a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Freshney, R. I. (2015). Culture of Animal Cells: A Manual of Basic Technique. New York, NY: Wiley.

Google Scholar

Glazier, D. (2015). Body-mass scaling of metabolic rate: what are the relative roles of cellular versus systemic effects? Biology 4, 187–199. doi: 10.3390/biology4010187

PubMed Abstract | CrossRef Full Text | Google Scholar

Helmlinger, G., Yuan, F., Dellian, M., and Jain, R. K. (1997). Interstitial pH and pO2 gradients in solid tumors in vivo: high-resolution measurements reveal a lack of correlation. Nat. Med. 3, 177–182. doi: 10.1038/nm0297-177

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, M., Dunn, J. C. Y., and Wu, B. M. (2005). Scaffold fabrication by indirect three-dimensional printing. Biomaterials 26, 4281–4289. doi: 10.1016/j.biomaterials.2004.10.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewis, M. C., MacArthur, B. D., Malda, J., Pettet, G., and Please, C. P. (2005). Heterogeneous proliferation within engineered cartilaginous tissue: the role of oxygen tension. Biotechnol. Bioeng. 91, 607–615. doi: 10.1002/bit.20508

PubMed Abstract | CrossRef Full Text | Google Scholar

Mac Gabhann, F., Ji, J. W., and Popel, A. S. (2007). VEGF gradients, receptor activation, and sprout guidance in resting and exercising skeletal muscle. J. Appl. Physiol. 102, 722–734. doi: 10.1152/japplphysiol.00800.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Magliaro, C., Rinaldo, A., and Ahluwalia, A. (2019). Allometric Scaling of physiologically-relevant organoids. Sci. Rep. 4:559682. doi: 10.1101/559682

CrossRef Full Text | Google Scholar

Malda, J., Rouwkema, J., Martens, D. E., Le Comte, E. P., Kooy, F. K., Tramper, J., et al. (2004). Oxygen gradients in tissue-engineered PEGT/PBT cartilaginous constructs: measurement and modeling. Biotechnol. Bioeng. 86, 9–18. doi: 10.1002/bit.20038

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, Y., and Vermette, P. (2005). Bioreactors for tissue mass culture: design, characterization, and recent advances. Biomaterials 26, 7481–7503. doi: 10.1016/j.biomaterials.2005.05.057

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattei, G., Giusti, S., and Ahluwalia, A. (2014). Design criteria for generating physiologically relevant in vitro models in bioreactors. Processes 2, 548–569. doi: 10.3390/pr2030548

CrossRef Full Text | Google Scholar

Mattei, G., Magliaro, C., Giusti, S., Ramachandran, S. D., Heinz, S., Braspenning, J., et al. (2017a). On the adhesion-cohesion balance and oxygen consumption characteristics of liver organoids. PLoS ONE 12:e0173206. doi: 10.1371/journal.pone.0173206

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattei, G., Magliaro, C., Pirone, A., and Ahluwalia, A. (2017b). Decellularized human liver is too heterogeneous for designing a generic extracellular matrix mimic hepatic scaffold. Artif. Organs 41, E347–E355. doi: 10.1111/aor.12925

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattei, G., Magliaro, C., Pirone, A., and Ahluwalia, A. (2018). Bioinspired liver scaffold design criteria. Organogenesis 14, 129–146. doi: 10.1080/15476278.2018.1505137

PubMed Abstract | CrossRef Full Text | Google Scholar

Nyberg, S. L., Remmel, R. P., Mann, H. J., Peshwa, M. V., Hu, W. S., and Cerra, F. B. (1994). Primary hepatocytes outperform Hep G2 cells as the source of biotransformation functions in a bioartificial liver. Ann. Surg. 220, 59–67.

PubMed Abstract | Google Scholar

Nyberg, S. L., Shirabe, K., Peshwa, M. V., Sielaff, T. D., Crotty, P. L., Mann, H. J., et al. (1993). Extracorporeal application of a gel-entrapment, bioartificial liver: demonstration of drug metabolism and other biochemical functions. Cell Transplant. 2, 441–452. doi: 10.1177/096368979300200602

PubMed Abstract | CrossRef Full Text | Google Scholar

Patzer, J. F. (2004). Oxygen consumption in a hollow fiber bioartificial liver–revisited. Artif. Organs 28, 83–98. doi: 10.1111/j.1525-1594.2004.07150.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Radisic, M., Malda, J., Epping, E., Geng, W., Langer, R., and Vunjak-Novakovic, G. (2006). Oxygen gradients correlate with cell density and cell viability in engineered cardiac tissue. Biotechnol. Bioeng. 93, 332–343. doi: 10.1002/bit.20722

PubMed Abstract | CrossRef Full Text | Google Scholar

Rivron, N. C., Liu, J, J., Rouwkema, J., de Boer, J., and van Blitterswijk, C. A. (2008). Engineering vascularised tissues in vitro. Eur. Cell Mater. 15, 27–40. doi: 10.22203/eCM.v015a03

PubMed Abstract | CrossRef Full Text | Google Scholar

Russell, S., Wojtkowiak, J., Neilson, A., and Gillies, R. J. (2017). Metabolic Profiling of healthy and cancerous tissues in 2D and 3D. Sci. Rep. 7:15285. doi: 10.1038/s41598-017-15325-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Shatford, R. A., Nyberg, S. L., Meier, S. J., White, J. G., Payne, W. D., Hu, W. S., et al. (1992). Hepatocyte function in a hollow fiber bioreactor: a potential bioartificial liver. J. Surg. Res. 53, 549–557. doi: 10.1016/0022-4804(92)90253-V

PubMed Abstract | CrossRef Full Text | Google Scholar

Sielaff, T. D., Nyberg, S. L., Rollins, M. D., Hu, M. Y., Amiot, B., Lee, A., et al. (1997). Characterization of the three-compartment gel-entrapment porcine hepatocyte bioartificial liver. Cell Biol. Toxicol. 13, 357–364. doi: 10.1023/A:1007499727772

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, M. D., Smirthwaite, A. D., Cairns, D. E., Cousins, R. B., and Gaylor, J. D. (1996). Techniques for measurement of oxygen consumption rates of hepatocytes during attachment and post-attachment. Int. J. Artif. Organs 19, 36–44. doi: 10.1177/039139889601900106

PubMed Abstract | CrossRef Full Text | Google Scholar

Tremper, K. K., Barker, S. J., Blatt, D. H., and Wender, R. H. (1986). Effects of anesthetic agents on the drift of a transcutaneous oxygen tension sensor. J. Clin. Monit. 2, 234–236. doi: 10.1007/BF02851171

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, B. A., Venkataraman, S., and Buettner, G. R. (2011). The rate of oxygen utilization by cells. Free Radic. Biol. Med. 51, 700–712. doi: 10.1016/j.freeradbiomed.2011.05.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, S., Cui, Z., and Urban, J. P. G. (2008). Nutrient gradients in engineered cartilage: metabolic kinetics measurement and mass transfer modeling. Biotechnol. Bioeng. 101, 408–421. doi: 10.1002/bit.21887

PubMed Abstract | CrossRef Full Text | Google Scholar