A sense of proximity: Cell packing modulates oxygen consumption

Accurately modeling oxygen transport and consumption is crucial to predict metabolic dynamics in cell cultures and optimize the design of tissue and organ models. We present a methodology to characterize the Michaelis–Menten oxygen consumption parameters in vitro, integrating novel experimental techniques and computational tools. The parameters were derived for hepatic cell cultures with different dimensionality (i.e., 2D and 3D) and with different surface and volumetric densities. To quantify cell packing regardless of the dimensionality of cultures, we devised an image-based metric, referred to as the proximity index. The Michaelis–Menten parameters were related to the proximity index through an uptake coefficient, analogous to a diffusion constant, enabling the quantitative analysis of oxygen dynamics across dimensions. Our results show that Michaelis–Menten parameters are not constant for a given cell type but change with dimensionality and cell density. The maximum consumption rate per cell decreases significantly with cell surface and volumetric density, while the Michaelis–Menten constant tends to increase. In addition, the dependency of the uptake coefficient on the proximity index suggests that the oxygen consumption rate of hepatic cells is superadaptive, as they modulate their oxygen utilization according to its local availability and to the proximity of other cells. We describe, for the first time, how cells consume oxygen as a function of cell proximity, through a quantitative index, which combines cell density and dimensionality. This study enhances our understanding of how cell–cell interaction affects oxygen dynamics and enables better prediction of aerobic metabolism in tissue models, improving their translational value.


INTRODUCTION
Oxygen (O 2 ) is a key player in cellular respiration.Its consumption rate (i.e., metabolism) is one of the major determinants of metabolic activity, which in turn strongly influences cell growth and fate. 1 O 2 's low solubility in aqueous media results in an inadequate supply of this primary energy resource to cells, particularly in high-density,three-dimensional (3D) tissue constructs, which are known to form hypoxic cores.][4][5][6] Investigating O 2 metabolism in cellular models-be they in silico or in vitro-is not straightforward because it depends on both intrinsic and extrinsic factors.][9][10][11] The size or extent of cell cultures, their density, and spatial arrangement in two or three dimensions (i.e., the dimensionality), as well as structural and functional characteristics of the substrate are also known to have a major impact on cellular activity. 12,13][16] A widely accepted analytical formulation for cellular O 2 consumption kinetics is the Michaelis-Menten (MM) model, 4,17,18 which is fully described by two parameters-the maximum single-cell O 2 consumption rate (sOCR, in mol s À1 ) and the MM constant (k M , in mol m À3 ).The MM can be used to express the O 2 consumption rate (R cell , in mol s À1 ) for a single cell as follows: where C (mol m À3 ) is the O 2 concentration perceived by the cell.Thus, the cellular consumption rate is adaptive with respect to ambient O 2 : at high O 2 levels (C ) k M ) a cell consumes at its maximal rate (i.e., R cell ffi ÀsOCR), while when resource availability is low (C ( k M ), the rate of consumption decreases proportionately (i.e., R cell ¼ À sOCR kM C).One of the best-known variables that conditions O 2 availability and, hence, consumption rate is the height of cell culture media in dishes, plates, or wells.Some theoretical studies and, as far as we know, only one experimental investigation 19 have sought to demonstrate the dependence of O 2 metabolism on the level of medium.It should also be noted that, although the MM kinetic parameters are often intended as intrinsic features of the specific cell phenotype, recent reports show that they are likely to vary with experimental conditions (e.g., cell density), 20 with direct implications on aerobic metabolism and, thus, cell and tissue function, as well as on related biophysical phenomena, such as size-dependent metabolic scaling. 21Nonetheless, most investigations on in vitro cultures disregard MM kinetics and consider cellular O 2 consumption rates to be constant and independent of differing ambient O 2 levels or other factors. 10,22he quest to engineer physiologically relevant models-that is, models which better recapitulate the behavior and function of in vivo tissues and organs-has propelled the use of 3D culture systems, such as scaffolds, spheroids, organoids, and microtissues.This in turn necessitates the design of constructs with adequate levels of O 2 supply to all cells, requiring more accurate tools and methods to quantify and therefore predict O 2 consumption in vitro.Do cells change their kinetic parameters when interacting with neighboring cells?Does the number of neighbors and their vicinity modulate their metabolic behavior?In short, can cells sense their local dimension, and can we exploit this knowledge to develop more physiologically relevant and, thus, translatable in vitro systems?
To address these questions, we designed an integrated approach to study the O 2 consumption characteristics of a human hepatic cell line, HepG2, 23 using the MM equation as a starting point.More in detail, using the experimental set-up in Fig. 1, O 2 concentration measurements recorded over time in two-dimensional (2D) monolayer cultures and 3D cell-laden spheroids were analyzed to extract reliable values of sOCR and k M through a purposely developed multiparameter identification algorithm.The MM parameters were then related to the cell density and dimensionality, quantified by means of imaging-based metrics.Using this approach, we show that kinetic parameters are not pre-determined but context-dependent, and this superadaptive behavior may underlie the emergence of cooperative and, hence, physiologically relevant behavior in specific culture conditions.

MM parameters and cell density
The values of sOCR and k M identified for different surface cell densities of monolayers are shown in Figs.2(A) and 2(B), respectively.For both parameters, the Kruskal-Wallis test highlighted that medians are not constant across densities (p < 0.0001 and p ¼ 0.002 for sOCR and k M , respectively).sOCR decreases with the surface cell density [Fig.2(A), Spearman coefficient r ¼ À0.8], while k M increases significantly [Fig.2(B), r ¼ 1].Note that k M values were not identifiable for the lowest cell density considered in the study, since the number of seeded cells was too low to guarantee a stable hypoxic condition, implying that k M is below the limit of sensitivity of our measurement setup.Hence, it was assumed equal to zero [Fig.2(B)], that is, cells in the monolayer consume O 2 according to a zero-order kinetics independent of ambient O 2 levels.][26] At higher densities-higher than those typically used in most experiments-we observe a notable increase in k M , which might indicate metabolic cooperation whereby the cells self-limit their consumption rate.Magliaro et al. report a similar trend in cylindrical cell-laden hydrogels. 20he dependency of sOCR on the volumetric cell density of spheroids reported in Fig. 2(C) could be interpreted in the same light.In fact, statistical differences emerging through the Kruskal-Wallis test (p < 0.0001) and the associated negative correlation with respect to volumetric cell densities (r ¼ À0.8) suggest cooperative behavior in terms of O 2 consumption, analogous to that noticed for 2D cultures [Fig.2(D)].Neither significant differences among medians (p ¼ 0.9030) nor correlations (r ¼ À0.2) were observed for k M .However, the median values of k M tend to increase for cell densities up to 5 Â 10 12 cell m À3 [Fig.2(D)], as might be expected for cooperative behavior.
It is worth noting that the dataset obtained for k M is characterized by a large dispersion, probably associated with the relatively small number of tested samples as well as unavoidable experimental errors.The estimation of k M is particularly sensitive to the latter because its value is often close to sensing limits. 27,28Such factors may mask differences, if any, related to cell density.We should also point out that the low k M values identified for the densest spheroids are likely to be an artifact due to the reduced cell viability that we noted (the fraction of viable cells was consistently $90% compared with !98% for the other FIG. 2 .MM parameters identified for monolayers (A) and (B) and cell-laden spheroids (C) and (D).In (B), k M !0 for q cell;s ¼ 3.1 Â 10 8 cells m À2 as it was not identifiable.Pairwise statistical differences through post hoc Dunn's multiple comparisons ( Ã p < 0.05, ÃÃÃ p < 0.0005, and ÃÃÃÃ p < 0.0001).5][26] All data are expressed as median 6 range.densities).Taken together, our data show a modulation of resource uptake driven by cooperation when O 2 depletion occurs in densely populated spheroids-statistically confirmed only for sOCR.][26] Concerning the mono-parametric estimation of MM parameters, reproducible estimations of k M were returned for some literature values of sOCR (data not shown).On the other hand, the algorithm generated widely dispersed and inconsistent estimations of sOCR or was even unable to identify this parameter within the threshold of tolerance if k M was fixed.This confirms that a global optimization approach like ours is necessary to properly estimate both MM parameters and further underlines the challenge of measuring k M .

Dimensionality and sphericity
To identify a common index to describe the extent of cell packing independent of the culture dimensionality, we extracted a proximity index, PI, through image processing of confocal stacks [Eq.( 6)]. Figure 3(A) shows PI values for the different surface and volumetric cell densities.Pairwise Dunn's comparisons resulted in a basically constant PI for all surface cell densities in the monolayers, while it increases notably with volumetric density for 3D cultures.Furthermore, the PI value of the monolayers was significantly lower than those obtained for spheroids.Irrespective of the slicing orientation, no relevant variations were observed when computing the PI on 10 lm cross-sectional slices in 3D, suggesting that the extent of packing can be reasonably assumed homogeneous within the same cell culture.
Nucleus sphericity was significantly different for the two dimensionalities investigated but did not depend on their density.Thus, data from all monolayers and spheroids were (respectively) lumped together.In coherence with the PI, the sphericity determined for the nuclei of HepG2 cells was significantly lower in 2D than in 3D [Fig.3(B)].While it is well-known that cells plated in monolayers tend to spread adopting a flatter shape than in 3D, our results suggest that the extent of flattening is irrespective of surface density.Similarly, cells encapsulated in the hydrogel matrix maintain a rounder shape and this roundness does not depend on the number of cells per unit volume.

Finding a common metric
We further introduced an uptake coefficient u ¼ sOCR PI = kM (lm 2 s -1 ), which can be interpreted as the area around a cell through which the O 2 it consumes passes per unit time.It expresses an analog of D for consumption and depends on the extent of cell packing (i.e., the combination of dimensionality and cell density) of the culture; the higher its value, the greater the area available to a cell.Thus, a dimensionless group expressing the local balance between u and D, which we refer to as the local supply-to-uptake ratio (LSU ¼ u/D) can be defined.The LSU is conceptually similar to the Thiele modulus, since it represents the ratio between the O 2 consumption and diffusion behavior of the culture, irrespective of its size, density, and dimensionality.The dependence of O 2 consumption kinetics on the PI, and the derived parameters, u and LSU, was investigated.Although statistical differences among medians for both sOCR and k M with respect to PI were observed, there were no significant trends or correlations.Refer to the supplementary material (Sec.SM5 and Fig. S3) for further details in this regard.
However, combining sOCR and k M to compute the uptake coefficient u, we noticed a significant correlation between u and corresponding PI values (r ¼ 0.7857, Fig. 4).After a steep raise for cell cultures with the lowest values of PI, the coefficient u saturates.An arbitrary exponential equation reasonably fits the u vs PI data (R 2 ¼ 0.5698), as reported in Fig. 4. Since u represents a measure of the surface per second related to O 2 uptake, this trend indicates that the ability of cells to recruit O 2 from the surroundings increases All data are reported as median 6 range.Ã p < 0.05, ÃÃ p < 0.005, ÃÃÃ p < 0.0005, and ÃÃÃÃ p < 0.0001.
rapidly when the cells are far apart (irrespective of dimensionality) but tends to a constant when they are tightly packed.Sigmoidal and Langmuir-type equations were also tested but resulted in weaker statistical significance (see supplementary material, Sec.SM6).
LSU is < 0.1 across all PI values, indicating that the diffusion locally remains faster than consumption even in highly populated cell cultures.In coherence with previous observations, this result suggests that when cells are close to their neighbors, they modulate their O 2 uptake to maintain a balance with respect to the diffusive supply.

DISCUSSION
Here, we demonstrate a novel method for systematically characterizing the O 2 consumption kinetics of cells across different dimensionalities and densities, starting from the well-known MM formulation.Our results underline the fundamental role of cooperation in shaping the O 2 metabolism of cell cultures. 29,30In fact, for all cultures, the values of sOCR and k M identified as well as the trends observed indicate that O 2 uptake is modulated as a function of (surface or volumetric) cell density (Fig. 2).
However, the contribution of dimensionality is impossible to infer unless a common metric that combines cell density and spatial arrangement can be identified.Here, we introduce a new imaging-based index-the so-called PI-for quantifying the closeness of cells regardless of their dimensionality [Eqs.( 5) and ( 6)].Based on this unique index, an uptake coefficient u, which combines the experimentally derived MM parameters and PI, was defined.u synthetically describes the ability of cells to "sense" their neighbors and, thus, modulate their O 2 uptake.In particular, u increases with PI for poorly packed cultures, up to a saturation level for the highest PI values (Fig. 4).Consequently, a proper supply-to-uptake equilibrium is locally guaranteed.We also introduced the LSU, a dimensionless number, which represents the ratio between u and D. Similar to Thiele's modulus, it quantifies the competing contributions of O 2 consumption and diffusion dynamics at the cellular level.We observed that LSU < 0.1 (Fig. 4), even in conditions of maximal packing, indicating that diffusion is faster than consumption.
In this study, we used HepG2 monolayers and HepG2-laden spheroids, measuring local O 2 consumption rates with microprobes or patches.While our data are in agreement with literature values for similar systems, the sOCR reported here differs significantly from maximal O 2 consumption rates per cell measured for perfused in vitro hepatic microtissues 22 and ex vivo rat and mouse livers [31][32][33][34] (our sOCR is about three orders of magnitude lower).However, these studies refer to hepatic tissues-be they engineered or from animal models-containing a variety of cell types present in the liver and rely on the measurement of O 2 concentration differences between the inlet and outlet of the perfusion circuit rather than pointwise monitoring of the O 2 level over time.Moreover, most reports 22,31,32 model O 2 consumption using zero-order kinetics, neglecting adaptation to local availability as described by MM kinetics.Despite these differences, where available, the half saturation constants estimated for perfused rat and mouse livers ($0.02 mol m À3 ) match our k M values, particularly for the higher volumetric densities.
It is well known that cells behave differently in 2-and 3D and that cell-cell interaction is fundamental to their function.In this study, we sought to identify a single metric capable of describing the volumetric O 2 consumption rate of hepatic cells as a function of their dimensionality and extent of cell-cell vicinity.The dependence of parameter u-and the corresponding dimensionless LSU-on the proximity index PI suggests that the metabolic behavior of the cells depends only on consumption parameters and the extent of cell packing rather than on culture conditions (e.g., O 2 availability, dimensionality).Although further investigationsmore experimental repeats covering a wider range of densities and the use of different cell phenotypes-are necessary to verify the statistical significance and generalizability of our data, this study demonstrates the importance of considering the O 2 consumption rate of a cell culture in terms of phenomena such as cell adaptation, cell-cell cooperation, and coexistence. 35These phenomena require a superadaptive kinetic formulation and cannot be explicitly described by means of the MM model, which only considers the dependence of consumption rate on the local O 2 concentration.consumption kinetics.In particular, measuring maximal rates of consumption (i.e., sOCR) in O 2 -saturated conditions and estimating the PI of cell cultures using well-established imaging techniques would be more feasible and reproducible than attempting to monitor O 2 metabolism in stable hypoxia.Further studies to elucidate the biochemical mechanisms underlying the observed superadaptive behavior could help better understand why sOCR and k M depend on the average proximity of cells.In fact, linking the modulation of MM parameters to the underlying biochemical pathways (e.g., hypoxia inducible factorrelated mechanisms) could contribute to the development of a more accurate mathematical formulation of O 2 consumption kinetics in cell cultures, which implicitly includes the dependency on cell proximity.Finally, a more exhaustive description of cell superadaptivity to both local O 2 concentration and cell packing might have implications beyond tissue engineering, from life sciences to metabolic regulation.

METHODS The mathematical model
Analytically speaking, measuring kinetic parameters requires monitoring the O 2 concentration field C ¼ C x; y; x; t ð Þover time in a region of the space with a fixed volume, where the generalized second Fick's law [Eq.( 2)] holds, where ) is the O 2 production/consumption rate per unit volume of interest, and D (m 2 s À1 ) is the diffusion constant of O 2 in water.In the case of consumption, R is negative and, as outlined in the introduction, typically described via the MM model.Thus, where q cell;v (cell m À3 ) is the cell density in the volume, and sOCR and k M are the MM parameters to be identified.Note that Eq. (3) assumes that sOCR and k M are constants, that is, each cell in the volume considered has the same metabolic sensitivity to O 2 levels and consumes at the same maximal rate, and that the cells are homogenously distributed in the volume.
In the case of monolayers, assuming that cells are homogeneously distributed over a plane (e.g., the bottom of a culture dish), O 2 dynamics is only dependent on the axial direction [i.e., C ¼ C z; t ð Þ, see Fig. 1(A)].Then, neglecting cell thickness, O 2 consumption can be modeled as a surface reaction, physically equivalent to an inward O 2 flux (Table I).Thus, Eq. ( 2) can be written as The challenge of measuring sOCR and k M Several practical considerations need to be considered to estimate the MM parameters for O 2 metabolism.Obviously, the sensor must have a response time, which is much faster than the kinetics being monitored and possess adequate sensitivity to changes in O 2 levels.Furthermore, to determine sOCR and k M , the cells must be probed over a range of O 2 concentrations.Specifically, to accurately estimate k M , they need to be exposed to 100% O 2 saturation conditions down to concentrations much less than k M .As k M is often much lower than ambient O 2 , this requires the cells to be brought to near anoxia.To achieve this condition, great care needs to be taken to ensure that the system is impermeable to atmospheric O 2 .We found that the standard method for achieving this (i.e., a layer of silicone oil atop the medium) is not sufficient, particularly when R is high.Thus, additional precautions for reducing O 2 diffusion at the air-media interface were purposely implemented as described herein.
Additionally, in 3D constructs, the rate of change of O 2 concentration varies in space and time, as expressed in Eq. ( 2).Hence, it must be locally monitored over time with high spatial resolution.Specifically, measuring changes in concentration over time requires a sensor with dimensions that are negligible in comparison to the volume of interest, and its positioning must be accurately controlled and known during experiments.On the other hand, 2D cultures require sensors able to measure O 2 levels in the vicinity of a surface.Finally, to make appropriate comparisons as a function of dimensionality, the sensors should ideally exploit the same physical principle.

Cell preparation
HepG2 cells are a human hepatocellular carcinoma cell line (ATCC V R , Manassas, Virginia, USA).While they lack drug-metabolizing enzymes, their endogenous metabolism is largely intact. 36They were cultured as described in the supplementary material (Sec.SM1); then, 2D and 3D cell aggregates were prepared as follows.
For monolayers, before cell seeding, the surface of the chosen flat sensor (see the sub-section "Oxygen concentration measurements") was coated with a 1% w/v alginic acid (Sigma-Aldrich V R , St Louis, Missouri, USA) solution containing 0.004% type I collagen (Life Technologies Corporation V R , New York, USA), followed by an overlaying volume of 20 ll of 0.1 M calcium chloride (CaCl 2 ) for gelation.After 15 min of incubation, the CaCl 2 solution was removed, and the surface of the hydrogel layer rinsed with phosphate buffered saline (PBS, Lonza V R ) 1Â.Then, HepG2 cells were suspended in 200 ll of Dulbecco's modified Eagle medium (DMEM, Sigma-Aldrich V R ), dispensed onto the prepared sensor and incubated overnight (37 C, 5% CO 2 , 95% RH), so as to obtain surface densities ranging from 3.10 Â 10 8 to 3.30 Â 10 9 cell m À2 , as reported in the supplementary material (Table S1).
Cell-laden spheroids were fabricated through the SpHyGa (spherical hydrogel generator) bioprinting platform (see supplementary material, Fig. S1). 37The cells were suspended at different volumetric densities, ranging from 1 Â 10 12 to 7.5 Â 10 12 cell m À3 (Table S2) in DMEM containing 1% w/v alginic acid and 0.004% type I collagen.The suspension was extruded as discrete droplets from the SpHyGa platform in an electrolytic bath of CaCl 2 0.1 M and incubated during alginate cross-linking for 15 min.The spheroids were moved to fresh DMEM and incubated overnight before testing.Details of the extrusion process are reported in the supplementary material (Sec.SM2), while the spheroid fabrication is depicted in Fig. 5 (and the associated video in the supplementary material).

Cell viability and density assessment
Before starting each batch of experiments, at least three samples were used to evaluate cell viability and the effective cell density.In the monolayers, cells were detached from the culture surface with trypsin-EDTA (Lonza V R , Basel, Switzerland); spheroids were left in a solution of 55 mM citric acid (Alfa Aesar V R , Haverhill, Massachusetts, USA) and incubated until complete dissolution of the hydrogel matrix.Then, the cells were centrifuged for 5 min at 900 rpm (Centrifuge 5702, Eppendorf V R , Hamburg, Germany) and resuspended in 200 ll of DMEM.
A sample of 30 ll was then mixed with the same volume of trypan blue (Thermo Fisher Scientific V R , Walthan, Massachusetts, USA), and an automated cell counter (Countess II FL, Thermo Fisher Scientific V R ) was used to determine cell number and viability.O 2 measurements were only performed on batches in which the percentage of viable cells was > 90%.
The cell surface density was obtained knowing the size of a 96well bottom.The volumetric cell density was instead determined from images of the spheroids acquired to compute their volume.More specifically, construct equivalent radii were estimated as the mean of those returned by using three different measurement tools in ImageJ. 38ygen concentration measurements O 2 concentration measurements were performed using commercial sensors based on O 2 quenching (Pyroscience V R GmbH, Aachen, Germany).A wide variety of sensing elements with the same functioning principle are available.
For cell monolayers, a contactless O 2 sensor spot (OXP5) attached to the bottom of a 96-well plate [Fig.1(A)] was used to acquire the average O 2 concentration at the cell level as a function of time.The well size was chosen to match spot dimensions, so that O 2 concentration is averaged over the whole surface of the cell culture.Sensor spots were fixed onto the plate bottom using a silicone glue (RS Components V R , London, UK), and optical signals were collected and transmitted to the processing unit through an optical fiber.For cellladen spheroids, a needle probe microsensor (OXR50), consisting of a needle housing a retractable optic fiber with a 50 lm sensing tip [Fig.1(B)], was directly connected to the processing unit.For both kinds of probe, a two-point sensor calibration was carried out in exactly the same conditions as those of the corresponding experiments for FIG. 5. Fabrication of cell-laden spheroids using the SpHyGa bioprinting platform. 37he video shows the whole process performed under a laminar flow hood: first, the filling of the syringe with the cell suspension; subsequently, the drop-by-drop extrusion and release of the suspension in the cross-linking bath to form spheroidal cellladen hydrogels (video available in the supplementary material).TABLE I. Initial and boundary conditions implemented for simulating O 2 dynamics for both dimensionalities (cell monolayers ¼ 2D, cell-laden spheroids ¼ 3D).n denotes the pointwise normal direction to the considered boundary.

Dimensionality
Boundary Condition Description Surface consumption 3D Air-medium interface, cuvette walls n Á J 6 L 2 ; y; z; t establishing the level of 100% O 2 saturation and the totally anoxic condition, using the protocol reported in Ref. 20.
For both 2D and 3D cultures, the air-medium interface was rendered O 2 -impermeable by combining the use of a thick layer of silicone oil [Agilent Technologies V R , St Clara, California, USA-Figs.1(A) and 1(B)] and an O 2 -consuming chemical filter (AnaeroGen, Thermo Fisher Scientific V R ).This combination was necessary because silicone oil alone-despite being advertised as O 2 impermeable-was not sufficient to guarantee the absence of O 2 diffusion into the culture medium in the timescales of interest.All experiments were carried out at 37 C using a purposely designed heating system composed of a transparent chamber made up of poly-methyl-methacrylate (PMMA, Plexiglass V R ) and an electronic module implementing a PID control.The chemical filter was left in the chamber to continuously scavenge O 2 and so maintain the boundary concentration at the air-medium interface close to zero.A water reservoir was also placed in the chamber to avoid evaporation of culture medium.Given the limited duration of measurements, the CO 2 partial pressure was not controlled.
For monolayers, the initial condition was set to 100% O 2 saturation by dispensing a minimum volume (20 ll) of fresh DMEM onto the cells, and the acquisition started immediately.The measurement terminated after a hypoxic steady state (i.e., O 2 saturation lower than 20%, corresponding to a concentration of 0.04 mol m À3 ) 5 was reached, which took from half an hour to a few hours depending on the cell surface density.Sensor spots were re-used for a few cycles, after cell detachment with trypsin-EDTA and rinsing with a 70% v/v ethanol solution.
Cell-laden spheroids were placed in a transparent cuvette containing 200 ll of fresh DMEM.A micro-manipulator and the camera of an optical tensiometer (Attension Theta Lite, Biolin Scientific V R , Goteborg, Sweden) were used to accurately position the sensor tip at the center of the spheroid [Fig.1(B)] until the stationary hypoxic condition was reached (implying measurements lasting from 10 min to an hour).Note that a single pointwise measurement over time of the local O 2 concentration at a controlled spatial position within the spheroid was enough to identify MM parameters, because the setup geometry (i.e., culture medium height, construct diameter) is known and-as shown in Fig. 6 (and associated video in the supplementary material)-the insertion of the probe did not affect the spheroidal shape of the construct.Between consecutive acquisitions, the sensing tip of the probe was rinsed with a 70% v/v ethanol solution.
For all experiments, several replicates were performed for each cell density and dimensionality (see supplementary material, Tables S1  and S2).

Kinetic parameter identification
The concentration profiles obtained were post-processed to estimate numerical values of the MM parameters through a multiparameter identification pipeline shown in Fig. 7. Basically, sOCR and k M were identified for each tested configuration by iteratively minimizing the mean squared error (MSE) computed between the measured O 2 concentration profile and the one obtained solving the nondimensional form of Eqs. ( 2) and (4) using the partial differential equations (PDE) toolbox in MATLAB (R2021b).All parameters as well as domain and boundary conditions used for modeling O 2 measurements in both cases are listed in Tables I-III, respectively.For each combination of sOCR and k M , the O 2 concentration profile over time was extracted and rescaled from the nondimensional concentration field returned by the PDE solver to evaluate the MSE with respect to the corresponding experimental dataset.Starting from reasonable initial ranges of the two MM parameters (Table II), the minimization problem was iteratively solved, and the ranges of parametric sweep correspondingly refined around the last estimated values.The algorithm stops and returns the final estimation of sOCR and k M when the MSE is below a predefined threshold of tolerance (set as 0.005 mol m À3 , $ an order of magnitude lower than the assumed hypoxic threshold) or stabilizes over two consecutive iterations (Fig. 7).
To assess whether a global fitting for the simultaneous determination of both MM parameters improves the quality of results, a (mono)parameter identification was also run; k M was set at literature values to estimate sOCR and vice versa.Section SM4 in the supplementary material reports further details and data on the multiparameter identification procedure.

Staining and imaging
For each cell density, three samples were fixed with a 4% w/v paraformaldehyde (PFA) solution and then stained with phalloidin red (Phalloidin-iFluor 647, Abcam V R , Cambridge, UK) and 4 0 , 6-diamidino-2-phenylindole (DAPI, Sigma-Aldrich V R ) using the manufacturers' protocols (supplementary material, Fig. S2).Constructs were left in PBS 1Â and stored in the dark at 4 C until image acquisition.Imaging was performed using a Nikon A1 confocal laser microscope (Nikon V R , Tokyo, Japan), equipped with a 10Â objective (pixel size: 0.31 lm).FIG. 6. Positioning of the optical microsensor at the center of the cell-laden spheroid.The video-directly recorded by the frontal camera of the tensiometer-shows the spheroid placed within the cuvette during the insertion of the sensor tip, which is accurately moved to the desired spatial point using the micromanipulator.The shape retention of the construct after puncturing can be also noticed.

Proximity index estimation
To identify a common metric to describe the extent of cell packing independent of the culture dimensionality (i.e., for both 2D and 3D cultures), confocal images were post-processed using ad hoc routines developed in MATLAB.Specifically, the blue channel was segmented to highlight cell nuclei (see the supplementary material, Fig. S2) and identify their centroid coordinates.Then, the weighted mean p i (lm À1 ) of the inverse distance of each centroid from all the other nuclei within the imaged region was computed.Sholl analysis 40,41 allowed the estimation of the number of cell nuclei n ij at a distance d j from the ith cell, representing the jth weight for averaging inverse distances from that cell.Thus, p i is given by the following equation:  where N d is the number of distances considered for the Sholl analysis to encompass the range [5.0; 179.2] lm with a step of 10 lm.The lower limit of the range corresponds to the radius of the smallest nucleus detected, while the higher limit is set by the size of the stack.Note that, for consistency, the total number of cells in the stack is N cell ¼ P N d j¼1 n ij .The proximity index PI (lm À1 ) was then defined as in the following equation, which is the average value of p i over all cells detected in the sample: For 3D spheroids, Eqs. ( 5) and ( 6) were used to compute PI values of both the whole construct and 10 lm-thick cross-sectional regions across its volume, obtained using different slicing orientations.The slice thickness corresponds to the diameter of the biggest nucleus detected.The PI accounts for both cell density and dimensionality and can thus be considered as a unique metric for comparing 2D and 3D cultures in terms of cell packing.We also estimated the average sphericity S (a.u.) of cell nuclei within each image stack using the following equation: 42 where A i (m 2 ) and V i (m 3 ) are, respectively, the surface area and the volume of the ith cell detected.They were automatically determined using the MATLAB's built-in connected component analysis.

Statistical analysis
Statistical analyses were performed in GraphPad Prism (v7) by means of a non-parametric Kruskal-Wallis test and pairwise post hoc Dunn's multiple comparisons for PI and sphericity.The Kruskal-Wallis test and related multiple comparisons, as well as nonlinear correlation (i.e., non-parametric Spearman coefficient) with respect to (i) the surface or volumetric cell density and/or (ii) the sample PI were used to analyze the dataset of identified sOCR and k M values.Nonlinear correlation with PI was also tested for u.
Note that, given the limited number of samples used, all quantities were expressed as median 6 range, and only non-parametric tests were performed.

SUPPLEMENTARY MATERIAL
See the supplementary material that contains supporting technical information on the cell maintenance, fabrication of cell-laden-spheroids, multiparameter identification and statistical analysis."Figures 1 and 3" of the supplementary material are videos of, respectively, the bioprinting process of spheroids and the microsensor positioning recorded via the frontal camera of the tensiometer.

FIG. 1 .
FIG. 1. O 2 sensing setup for the monolayer and 3D configurations.(A) The 96-well plate containing the O 2 sensor spot where HepG2 monolayers were cultured.(B) HepG2laden spheroid punctured with the O 2 microsensor, frontal view through a transparent cuvette.A schematic representation highlighting the main components and O 2 transport directions is reported for both configurations.

CONCLUSIONA 2 FIG. 4 .
FIG.4.Outcome of the multiparameter identification procedure expressed through the uptake coefficient (u, left vertical axis) and the local supply-to-uptake ratio (LSU, right vertical axis) as a function of the PI.An arbitrary exponential fit of u as a function of PI (black dashed curve, R 2 ¼ 0.5698) is reported.All data are expressed as median 6 range, computed from measured values of sOCR, k M , and PI by applying standard rules for error propagation.Datapoints were fitted weighting by the corresponding extents of dispersion.For the sake of visualization, the vertical axes are in logarithmic scale.Other fittings are reported in the supplementary material (Sec.SM6).

FIG. 7 .
FIG.7.Pipeline of the multiparameter identification workflow for sOCR and k M .For each sample, the optimal values of sOCR and k M are determined by iteratively simulating the O 2 transport and consumption dynamics of the system using the PDE toolbox in MATLAB, until the MSE (mean squared error) between predicted and measured O 2 concentration profiles is minimized.
a Median of values identified in Ref.20, in which gels and cell densities were similar to those used here.