Oxygen Consumption Characteristics in 3D Constructs Depend on Cell Density

Oxygen is not only crucial for cell survival but also a determinant for cell fate and function. However, the supply of oxygen and other nutrients as well as the removal of toxic waste products often limit cell viability in 3-dimensional (3D) engineered tissues. The aim of this study was to determine the oxygen consumption characteristics of 3D constructs as a function of their cell density. The oxygen concentration was measured at the base of hepatocyte laden constructs and a tightly controlled experimental and analytical framework was used to reduce the system geometry to a single coordinate and enable the precise identification of initial and boundary conditions. Then dynamic process modeling was used to fit the measured oxygen vs. time profiles to a reaction and diffusion model. We show that oxygen consumption rates are well-described by Michaelis-Menten kinetics. However, the reaction parameters are not literature constants but depend on the cell density. Moreover, the average cellular oxygen consumption rate (or OCR) also varies with density. We discuss why the OCR of cells is often misinterpreted and erroneously reported, particularly in the case of 3D tissues and scaffolds.


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. 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 2 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 reoxygenated 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 −1 · s −1 ], typically estimated by measuring the overall oxygen consumption rate (in mol·s −1 ) in the system and dividing by the total number of cells present. The range of OCR is usually between 1 × 10 −16 and 1 × 10 −18 mol·cell −1 · s −1 (or 100 to 1 amol·cell −1 · s −1 ) (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., 1993Nyberg et al., , 1994Sielaff 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 −3 · s −1 ]) is described by the well-known Michaelis-Menten law (Equation 1): 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 −3 ]). 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 −3 ]) 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.

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 R 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 noninvasive/non-destructive measurements through a transparent material without consuming oxygen (measurement range: 0-21% O 2 ; accuracy: ±0.01% O 2 ; 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 −3 (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 2 partial pressure, i.e., 0 mol·m −3 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 2 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 12 , 1·10 12 , and 5·10 12 cells·m −3 ) 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 R 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.

Oxygen Consumption Measurements
Oxygen measurements were performed by placing the setup inside a 37 • C/5% CO 2 humidified incubator. At the beginning of the experiment (t = 0), hydrogel constructs polymerised and maintained in the 20% O 2 incubator (thus containing 20% O 2 concentration, in equilibrium with the overlying culture medium) were exposed to a step change in the O 2 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 −2 mol·m −3 , instead of the expected 0 mol·m −3 , 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. 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.

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: where C [mol . m −3 ] is the concentration of oxygen as a function of space and time, D [m 2 · s −1 ] is the diffusion coefficient of oxygen in the gel (assumed to be spatially homogeneous), and R [mol O 2 consumed·m −3 · s −1 ] 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 V max = ρ cell ·sOCR and considering only oxygen concentration variations in time and height, i.e., C(z, t), obtaining: In Equation (3), the sOCR (mol·cell −1 · s −1 ) 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) = V max = ρ 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 −1 · s −1 ]), 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): (4) Given the cylindrical symmetry of the problem investigated in this work, Equation (4) simplifies to: 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 2 in cell culture incubators due to the presence of 5% CO 2 . Because all the hydrogels were prepared under ambient conditions and then polymerised and maintained in a 20% O 2 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 = 0.2 mol·m −3 . Conversely, the O 2 concentration in the overlying culture medium was experimentally found to be C b = 1% O 2 (section Oxygen Consumption Measurements), thus giving the boundary condition C (z = H, t) = 1·10 −2 mol·m −3 . 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:

Non-dimensionalisation of the Governing Equation
The governing diffusion-reaction equation (Equation 2) was rendered non-dimensional by introducing the following group of variables: 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): Equation (3) can then be rewritten in the following dimensionless form: which is subjected to the initial and boundary conditions below: A summary of all the model parameters is presented in Table 1.

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.

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).
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).
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 2 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 deathreasonable 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.
The assumption that ∂c m ∂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 2 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 2 concentration decreases and returns a time-dependent parameter, differently than the classical Thiele modulus which is a constant (Equation 10) (Fink et al., 1973): 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: 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).

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.