A Causal Relation between Bioluminescence and Oxygen to Quantify the Cell Niche

Bioluminescence imaging assays have become a widely integrated technique to quantify effectiveness of cell-based therapies by monitoring fate and survival of transplanted cells. To date these assays are still largely qualitative and often erroneous due to the complexity and dynamics of local micro-environments (niches) in which the cells reside. Here, we report, using a combined experimental and computational approach, on oxygen that besides being a critical niche component responsible for cellular energy metabolism and cell-fate commitment, also serves a primary role in regulating bioluminescent light kinetics. We demonstrate the potential of an oxygen dependent Michaelis-Menten relation in quantifying intrinsic bioluminescence intensities by resolving cell-associated oxygen gradients from bioluminescent light that is emitted from three-dimensional (3D) cell-seeded hydrogels. Furthermore, the experimental and computational data indicate a strong causal relation of oxygen concentration with emitted bioluminescence intensities. Altogether our approach demonstrates the importance of oxygen to evolve towards quantitative bioluminescence and holds great potential for future microscale measurement of oxygen tension in an easily accessible manner.


Introduction
In situ studies on the mechanisms of cell fate regulation in local microenvironments (niches) has gained considerable interest in the development of cell based therapies for disease and regeneration. These studies are very often complemented with bioluminescence imaging assays that yield valuable information on cell fate and behavior in a dynamic microenvironment. Likewise, specific niche components have been screened for their contribution to therapy outcome, including the matrix elasticity [1], presence of soluble and matrix-bound chemical agents in biomaterials [2], targeted and sustained release of cytokines from transplanted cells [3], and cell-adhesion mediated resistance against therapeutic agents [4].
Bioluminescence imaging relies on the activity of luciferase enzymes that act as catalyst for the conversion of luciferin to oxyluciferin, which is accompanied by the release of a photon [5]. Even though bioluminescence reporter imaging is a well-integrated technique for probing the biological function of living cells in vitro as well as in small-animal models [6,7], the complexity of 3D cellular microenvironments precludes a quantitative interpretation of bioluminescent light [8]. Here we report on three major sources for the ambiguity in interpreting bioluminescence data obtained from cell-seeded hydrogels. At first the availability of bioluminescence substrate molecules (luciferins) to the luciferase enzymes which is dependent on active and/or passive transport from the site of application, secondly the accessibility to oxygen that is required for substrate oxidation, and finally the need for dynamic time point measurements of luciferase activity are all crucial determinants for making robust quantitative analyses. Although previous studies have elucidated the effects of oxygen concentration on the emitted bioluminescence intensity [9,10], we show here how a mathematically validated model aids in resolving the oxygen dependent influences (changes in intensity, Michaelis-Menten kinetics, and decay rates) and how this model can be used to obtain quantitative measurements of the intrinsic bioluminescence intensity. Furthermore, we show that by careful analysis of the bioluminescence signal information can be obtained on local oxygen concentrations, evidencing a causal link between bioluminescence and oxygen.

Bioluminescence Intensities are Influenced by the Available Oxygen Concentration
The importance of the oxygen availability for the bioluminescence reaction is best illustrated by the light flashing mechanism in the adult firefly Photinus pyralis. Tracheolar fluid levels in this system are acting as a diffusive barrier for oxygen supply and effectively control the light flashes emitted from photocytes in the firefly's abdomen [11]. Oxygen-dependent luciferase activity is also observed in vitro in cell monolayer assays and has been used as reporter for cellular oxygen availability upon incubation with nitric oxide [12]. Here, we show that hypoxic conditions applied to firefly luciferase solutions result in a ,3.4 fold difference in total photon flux as compared to normoxia (Fig. 1A). With the resolving power of currently available imaging equipment this fold difference should enable the quantification of oxygen-dependent luciferase activity. Conditions of normoxia (21% O 2 ) and hypoxia (near 0% O 2 ), representing the boundaries of oxygen availability in the artificial (atmospheric) and physiological (in vivo) cell environment [13], were applied in our experimental studies. Intermediate levels of oxygen concentration, which are difficult to maintain using regular gas exchange setups, were not directly applied but their influence was adapted from previously reported relations (which will be described below).
Other cofactors of the bioluminescence reaction (e.g. ATP) also account for a decrease in measured photon flux, but can often be directly or indirectly related to the influence of oxygen. Gradual accumulation of the inactive dehydroluciferyl-adenylate (L-AMP) complex within the cytoplasm of intact cells results in a lower photon flux compared to the free luciferase solution [14,15]. The availability of free diffusing luciferin may be further decreased as membrane-bound ABC transporters not only control the net influx of luciferin into the cytoplasm but also require luciferin as a substrate for their activity [16]. In addition to the lower luciferase activity, ABC transporter-luciferin binding is also reflected in a slower apparent diffusivity of luciferin within cell-seeded hydrogels (Fig. S2). Comparison of photon fluxes from luciferase solutions relative to intact cells for corresponding luciferase concentrations yielded a difference of ,3.5 fold under normoxic conditions, while a much larger difference of ,16.5 fold was observed for the hypoxic environment (Fig. 1A). This effect mainly originates from a decrease in intracellular ATP content under hypoxia, concomitant with a reduction in mitochondrial membrane potential [10,17]. We support this reasoning by visualization of mitochondria with a MitoTracker dye, clearly indicating a strong reduction of stained mitochondria in case of hypoxic incubation (Fig. 1B,C). Lower intracellular ATP concentrations potentially also influence the activity of ABC transporters, leading to a reduced influx of luciferin into the cytoplasm and hence a reduced photon flux [16].

Reduced Oxygen Concentrations Induce Changes in the Bioluminescence Reaction Kinetics
Oxygen dependent changes in initial bioluminescence reaction kinetics were determined from dynamic time point measurements of luciferase activity with varying luciferin concentrations. As peak intensities are reached within less than 1 second after reagent addition, fast operation and manipulation would be required to monitor initial light emission [18]. To circumvent these practical issues specifically for cell lysates, and avoid potential signal interference from measurement equipment, we extrapolated bioluminescent data that were obtained at later time points ( Fig. 2A,B). Typical exponential decay of luciferase activity was observed under both normoxic and hypoxic conditions. Initial reaction velocities were subsequently transformed into Lineweaver-Burk plots in order to retrieve Michaelis-Menten kinetic parameters (Fig. 2C). Under normoxic conditions a higher substrate affinity was found as compared to hypoxia, with corresponding average Michaelis-Menten parameter values (K m,21 and K m,0 ) of 42.6 and 938.4 mM respectively (Table S1). To estimate initial reaction velocities at intermediate oxygen tensions a square root dependency was adopted, similar to what has been reported for ATP-dependent kinetics [19]. This approach resulted in an oxygen dependent luciferase activity that matched with data from Moriyama et al. [10]. In this work, a previously published ordinary differential equation (ODE) model was extended to describe the average photon flux from cell lysates or intact human embryonic kidney 293T cells [20]. The extended model describes the kinetics of luciferin (extracellular as well as intracellular), firefly luciferase, and emitted photon flux, taking into account the oxygen dependency of photon generation (see Text S1 for further details). For intact cells we implemented an apparent Michaelis-Menten parameter (K m ) which is 1.5 fold higher than in cell lysates. This value was obtained by model fitting and reflects the influence of a highly crowded cytoplasm on the luciferase activity [21,22]. These results show great discrepancy with available literature data, where ,7 fold changes in substrate affinity have been reported for cell lysates versus intact cells [23]. Most probably this literature data discrepancy can be attributed to the influence of luciferin membrane transport on the luciferase activity in living cells. In our analysis this influence was included explicitly, by introducing an additional transport term in the model equations (equation (1) and (2) in Text S1).

Bioluminescence Decay Rates are Modified at Reduced Oxygen Concentrations
Intact cells show a significantly slower bioluminescence signal decay compared to the decay in luciferase solution, which is orchestrated by the concentration of cytoplasmic inorganic pyrophosphate (PPi) [14]. Also coenzyme A (CoA) concentrations can influence the decay dynamics. CoA stabilizes the photon flux by thiolysis of L-AMP into dehydroluciferyl-coenzyme A (L-CoA), which is a less powerful inhibitor than L-AMP on the bioluminescence reaction [24]. Under normoxic conditions the average exponential decay rates changed from 3948 s 21 for firefly luciferase extracted from cell lysates to 180 s 21 for intact cells. Hypoxic conditions have a remarkable influence on these decay rates with a decay rate of 1452 s 21 for cell lysates and 3 s 21 for intact cells (Fig. 2D). Cells exposed to an initial luciferin substrate concentration of 470 mM presented a bi-exponential decay in their luciferase activity (Fig. 2E,F) [23]. A fast initial decay (peak with decay rate of 180 s 21 ) which takes place in the first few minutes and which is indicative of a rapid substrate exhaustion, was followed by a slower decay rate of 27 s 21 . Hypoxic environments repressed peak occurrence and transformed the diffusion-limited bioluminescence reaction into a more transition state-limited reaction. These changes were explicitly implemented into the model by a decrease in luciferin substrate affinity under hypoxia (K M , Table S1) and the introduction of a transition threshold parameter (k, Text S1).

Substrate Transport can be Monitored by Fluorescence Recovery after Photobleaching (FRAP)
As indicated before, active cellular membrane transport is not the sole mechanism which controls the availability of luciferin to the cytoplasmic luciferase. When cells are embedded in an extracellular matrix (ECM), luciferin molecules will be transported through this matrix prior to intracellular uptake (Fig. 3A). This necessitated the introduction of diffusion in the model and therefore resulted in a system of partial differential equations (PDEs). Using fluorescence recovery after photobleaching (FRAP) [25], we were able to obtain quantitative measures of substrate diffusion rates at different positions in space and time, in order to rule out any changes in diffusivity. Measurements were performed with an optimized concentration of fluorescein tracer (see Fig. S2). FRAP-based ratios of tracer diffusivity in culture medium versus agarose hydrogel were validated by comparing diffusion rates obtained from Fluorescence Correlation Spectroscopy (FCS) (Fig.  S3).

Oxygen independent Bioluminescence Intensities are a Good Measurement of the Active Bioluminescent Cell Population
Oxygen-and transport-related parameters influencing luciferase activity were subsequently implemented into the PDE-based model. Time-dependent luciferase activity showed a typical exponential decay which was described mathematically by a first order chemical kinetic equation (Equation (3) in Text S1) [26]. Substrate diffusion through cell membranes also accounted for the observed time-dependency of bioluminescent reaction rates, which by virtue of the induced spatial heterogeneities resembled fractallike reaction kinetics [27]. The afore-mentioned model fits the experimental data points with good accuracy, thereby enabling reliable quantitative insights into the bioluminescence reaction kinetics for both cell lysates and intact cells (Fig. 2D,E). The influence of spatial organization on the emitted bioluminescence from luciferase reporter cells was validated for a homogeneously cell-seeded (,1610 6 cells?ml 21 ) cylindrical agarose hydrogel (Fig. 3B). Circular glass plates confined the hydrogel in axial direction, preventing axial diffusion of oxygen and luciferin. This resulted in an axisymmetric setup which simplified further analysis.
Microscale measurements of oxygen concentration were performed by integration of fluorescent oxygen sensitive microbeads that we previously developed [28] (Fig. S4). Stable read-outs of oxygen concentration over a time period of 3 days were obtained for control gels in which no cells were embedded (Fig. S4F). Radial oxygen concentration profiles from cell-seeded hydrogels displayed a gradual decrease towards the center with an averaged availability of oxygen that increased over time (Fig. 3C). Diffusivity measurements obtained from FRAP experiments in cell-seeded hydrogels displayed no significant changes in tracer diffusion rate, with an average value of ,350 mm 2 s 21 (Fig. 3D). As nearly  constant values in space and time for the tracer diffusion rate were measured it can be concluded that the time-dependent increase in radial oxygen concentration profiles was induced by a decrease in initial oxygen consumption rate (OCR) by the cells with time (Fig.  S4G). 293T cells are able to respond to reduced oxygen availability by altering cytochrome c oxidase (COX) subunit composition via a hypoxia-inducible factor 1 (HIF-1) -dependent mechanism [17], that results in a decrease of their OCR. We implemented these time-dependent changes in OCR directly into the model to accurately match the measured radial oxygen profiles ( Fig. 3C and equation (4) in Text S1).  Luciferin is actively transported across the cell membrane (thickness, l c ) and reacts with luciferase (ii) in the cell cytoplasm where this reaction is accompanied by the release of a photon (i). Oxygen availability in the cytoplasm modulates emitted light intensity and kinetics and is described by the Michaelis-Menten kinetics. (B) Setup for validation of the bioluminescence-oxygen model. Oxygen Sensing microBeads (i) and 293T cells (ii) are embedded in an agarose gel that is confined between 2 circular glass plates. Focal volume (v) imaged by combined bioluminescence and fluorescence microscopy reveals luciferase activity in single intact cells (iii) and local oxygen concentrations based on ratiometric intensities obtained from oxygen sensitive and insensitive dyes (iv). Scale bar, 1mm. Scale bar figure insets, 10 mm. Dynamic time point measurements of the average photon flux emitted from the cell-seeded hydrogels revealed the presence of a fast bioluminescence intensity peak that fades away slowly with a decaying signal that was detectable for several hours (Fig. 4A). Peak intensities reached maximum values after 2 days of cultivation and decreased gently with longer cultivation time (Fig. 4B). These results did however not entirely correlate to the quantitative DNA and viability measurements, representing the active bioluminescent cell population and clearly being oxygenindependent readouts (Fig. S5). To decouple bioluminescence signal reshaping caused by oxygen availability from the spatial distribution of luciferase activity, we removed cell respiration in our model, hence maintaining a saturated (and uniform, 21%) oxygen level within the hydrogel and therefore oxygen-independent bioluminescence signal intensities. A prominent role for oxygen in reshaping the bioluminescence signal was clearly observed (Fig. 4C-E). When saturated oxygen conditions were simulated the obtained peak signal intensities matched closely the cell activity measurements (Fig. 4F). Radial profiles of the bioluminescence signal emitted from encapsulated 293T cells indicated high activity near the hydrogel edge at early time points that gradually leveled off towards the center after longer measurement times (Fig. 4G, Video S1). With increasing cultivation times the observed effect of high edge activity became more pronounced and was in good agreement with our simulation results. Combined, these data show a strong influence of oxygen availability on the activity of luciferase reporter cells which should be taken into account when intrinsic (i.e. independent of other environmental factors) measurements and error-free interpretations of bioluminescence intensity are pursued.

Oxygen Concentrations within Cell-seeded Hydrogels can be Derived from Bioluminescence Signal Analysis
With longer cultivation times peak bioluminescence also occurred faster after luciferin substrate addition, in conjunction with a reduced full width at half maximum (FWHM) and a faster decay rate of the bioluminescence signal ( Fig. 4B and Fig. S6A-C). Bioluminescence microscopy of cell activity at the hydrogel center revealed comparable parameter trends, though the change in signal decay rate with culture time was more pronounced (Fig. 4H,I and Fig. S6E-G). Strong discrepancy with the simulated decay rate was observed at 3 days of cultivation. This difference could therefore indicate that a luciferin concentration dependent decay would be required to accurately fit the experimental data, which is supported by a strongly reduced substrate availability in the hydrogel center (Fig. S7) and maintenance of a concentration dependent substrate gradient (partitioning) across the cell membrane (Fig. S8) [14]. A good correspondence in the evolution of bioluminescence peak parameters and available oxygen concentration was clearly observed, evidencing the causal link (Fig.  S6D,H). These results hence show great promise towards quantifying averaged or localized oxygen concentrations in cellbased systems that are based on univocal relations between bioluminescence and oxygen.
Oxygen is a critical component of the cell microenvironment, such as for the determination of cell fate in stem cell niches [29,30]. As such, new non-invasive technologies and methodologies should be developed that are capable of oxygen measurement and control in 3D environments in space and time. In this respect, the methodology as stated above holds great potential as a tool for easily-accessible measurements of oxygen concentration directly from the time-dependent bioluminescence signal emitted within a dynamic microenvironment. Further validation is required to assess practical use of the identified bioluminescence peak parameters in determining oxygen for other experimental setups.
293T cells were transduced with a lentiviral vector (pCH-EF1a-3flag-fLuc-T2A-eGFP-Ires-Bsd, 3.1610 8 TU?ml 21 ) which was a kind donation from Dr. Greetje Vande Velde (MoSAIC, KU Leuven). The day before transduction, cells were seeded in a 96well plate at 1610 4 cells per well. On the day of transduction, medium was replaced by DMEM containing serial dilutions of the vector and incubated for 24 hours. After 24 hours, medium was replaced with DMEM containing 1 mg?ml 21 blastidicin for antibiotic selection of the stably transduced cell population, and was continued for 2-3 weeks. Transduction efficiencies were analyzed by flow cytometric analysis (FACS). Luciferase expression levels of the transduced cells were stable over time (Fig. S1).

In vitro Luciferase Activity Assay
Stable read-outs of luciferase activity were obtained by thoroughly mixing 5 ml solutions of recombinant luciferase (QuantiLum, Promega), firefly luciferase extracted from 293T cells or intact luciferase-reporting 293T cells in suspension in 45 ml of Luciferase Assay Reagent (Promega). Incubation resulted in a glowtype reaction that lasted longer than ,1 min and total photon fluxes were measured using black-walled 96 well plates in the IVIS 100 imaging system (PerkinElmer) with 1 s acquisition time. Image processing and analysis were performed in IGOR Pro software (WaveMetrics). Final concentrations of recombinant luciferase solutions were obtained upon dilution of stock solutions with purified water from a Milli-Q system (Millipore) (Fig. S1). Measurements on cell lysates were performed by lysing 293T cells in Cell Culture Lysis Reagent (Promega) with a final concentration of 1610 5 cells?ml 21 and were compared to an equal concentration of 293T cells in DMEM (Fig. S1). Hypoxic conditions were induced by addition and incubation (.15 min) of 1% Na 2 SO 3 in luciferase solution and assay reagent before mixing. All bioluminescence measurements were performed at 37uC incubation temperature.

Dynamic Time Point Analysis of Free Luciferase Activity
Flash-type bioluminescence reactions were provoked by mixing 10 ml solutions of recombinant luciferase (0.1 nM) with 100 ml luciferin solutions (Beetle luciferin potassium salt (Promega) in 10 mM tricine, 1.07 mM magnesium carbonate pentahydrate, 2.67 mM magnesium sulfate, 0.1 mM ethylenediaminetetraacetic acid (EDTA), 0.5 mM ATP sodium salt, 0.27 mM coenzyme A sodium salt and 33.3 mM dithiothreitol (DTT)). Initial reaction velocities were obtained by fitting exponential decay functions through the measured data points and by extrapolation of the luciferase activity to the initial time point. The initial time point was found by application of a ratiometric criterion between total emitted photon fluxes at normoxia and hypoxia, that we measured via a stable glow-type reaction. This ratio indicated ,3.37 fold difference in initial reaction velocities. Photon fluxes were measured using black-walled 96 well plates in the IVIS 100 imaging system (PerkinElmer) with 1 s acquisition time.

Dynamic Time Point Analysis of Luciferase Activity in Intact Cells
Cell densities (1.09610 5 cells?ml 21 ) used for dynamic time point analysis were calculated from the average intracellular luciferase concentration (0.92 amol?cell 21 ) as the equivalence of a 0.1 nM free recombinant luciferase solution. Cells were suspended in DMEM and at ,20 min before signal measurement, and cells used for the hypoxic condition were resuspended in DMEM containing 1% Na 2 SO 3 . Longer hypoxic incubation times resulted in a modest decrease in luciferase activity. Luciferase activity was measured after mixing 10 ml of cell suspension with 100 ml of luciferin solution (Beetle luciferin potassium salt in DMEM). Photon fluxes were measured using black-walled 96 well plates in the IVIS 100 imaging system (PerkinElmer) with 1 s acquisition time. The stability of luciferase expression levels for stably transduced cells that are exposed to different oxygen concentrations has been shown elsewhere [10].

Single-cell Bioluminescence Microscopy
Glass coverslips were coated with 500 ml Poly-L-Lysine (PLL, Sigma) solution (0.1%, 5 min), rinsed with Milli-Q water and dried overnight. We plated 1610 5 293T cells on PLL-coated coverslips. After overnight cell attachment these plates were placed on the stage of a luminescence microscope (LuminoView 200, Olympus). Cells were incubated with a 470 mM luciferin solution at 37uC (Solent Scientific). Bioluminescence was imaged with a UPLSAPO 606 water objective (NA: 1.2) and transmitted to a cooled CCD camera (ImagEM512, Hamamatsu Photonics) mounted on the bottom port of the microscope. Time-lapse images were collected with 5 s acquisition times and EM gain was set at 12006 (photon imaging mode, Hamamatsu) (Fig. 2F).

Fluorescence Recovery after Photobleaching (FRAP) Analysis of Self-diffusivity
Low melting point agarose gels (Invitrogen) with a final concentration of 2% were prepared on glass coverslips (thickness, 2 mm and diameter, 8 mm). Gels were incubated with DMEM containing various concentrations of fluorescein (Sigma). After overnight incubation slides and gels were transferred to the stage of a confocal fluorescence microscope (FluoView 1000, Olympus) equipped with a UPLSAPO 106 air objective (NA: 0.40) used for observation. Measurements were performed at 37uC. Fluorescence images were collected with the systems PMT. A 488 nm Ar laser was used to bleach and also monitor the recovery of fluorescein tracer molecules. The pinhole size was set to 50 mm. Images were acquired with 72 ms intervals during a 3 s scanning period, and consisted of 2566256 pixels with a pixel size of 0.49760.497 mm (zoom factor: 106). Prebleaching images were acquired to compensate for non-uniform light illumination. Circular regions (radius, 10 mm) were bleached (total bleaching time, 13 ms) with the laser in tornado-scan mode (SIM Scanner, Olympus) to obtain centered, fast, and uniformly bleached spots with an approximate Gaussian shape. Image analysis was performed with a program written in MATLAB (The MathWorks, Natick, MA). The method implemented into this program is based on a spatial frequency analysis of circularly averaged radial intensities of each image [25]. In brief, the recovery of fluorescent tracer was modeled according to Fick's second law. An analytical solution to this equation was obtained via the Hankel transform. Circular averaging on the radial intensities was performed to reduce the noise in the intensity profiles. Finally the analytical solution was fitted to the experimental curves using a nonlinear curve fitting algorithm in MATLAB. We also verified the independence of diffusion rates on spatial frequency, as our setup is characterized by Brownian diffusion. Analysis was performed for a single diffusing component with the fraction of immobile molecules set to zero. The fluorescein tracer concentration was optimized to detect changes in fluorescein diffusivity, caused by local matrix degradation, with maximal sensitivity (Fig. S2B). A fluorescein concentration of 25 mM was found to be optimal. Spatial heterogeneities in tracer diffusivity were detected by FRAP imaging at different circumferential positions along 3 radial positions in the gel (imaging depth, 150 mm). Measurements were performed in duplicate on each gel with a small shift in spatial position.

Fluorescence Correlation Spectroscopy (FCS)
Agarose gels (2% in DMEM supplemented with various concentrations of Rhodamine B) were prepared between two glass coverslips separated by ,450 mm thick spacers. FCS measurements were performed on a custom-built fluorescence spectro/microscopy setup using a 532 nm laser line (Spectra-Physics CW DPSS Excelsior laser; the maximal laser power of 500 mW was attenuated by means of neutral density filters (New Focus)) for fluorophore excitation. The excitation light was circularly polarized by means of a l/4 waveplate (Thorlabs), expanded using a telescope arrangement (to a collimated beam of about 7 mm diameter (1/e 2 intensity), and directed via a dichroic beamsplitter (z532rdc, Chroma Technology, Rochingham, USA) into the oil-immersion objective (UPLSAPO 1006, NA = 1.4) of an inverted optical microscope (Olympus IX71). The center of the confocal volume was positioned 5 mm above the coverslip surface. Fluorescence light was collected by the same objective and was guided, after passing through the dichroic, through a 50 mm pinhole (Linos/Qioptic) in order to reject light from out-of-focus regions. The emitted light was filtered by a longpass filter (HQ545LP, Chroma Technology, Rockingham, USA) focused on a t-SPAD photon counting module (PicoQuant, active area 150 mm in diameter), connected to a HydraHarp 400 TC-SPC module. Measurements were performed until at least 3 million photon-detection events had been recorded, typically requiring a few minutes. Calibration of the one-focus FCS setup was performed using a solution of Rhodamine 6G in Milli-Q water (D < 372 mm 2 s 21 ) [31] using the same excitation laser power as for the samples under investigation (1.69 kW cm 22 at the sample). Temporal autocorrelation curves were calculated from the measured fluorescence intensity fluctuations (photon count time traces) within the SymPhoTime software (PicoQuant). These autocorrelation curves were fitted in IGOR Pro (WaveMetrics) to an expression for autocorrelation decay due to anomalous diffusion [32], In this expression, AENae is the average number of fluorescent dye molecules in the confocal volume element, t D is the diffusive time of dye molecules, a is the degree of anomalous subdiffusion, and v considers the extension of the confocal volume along the optical axis ( = AEw z ae/AEw xy ae). The diffusive time is related to the translational diffusion coefficient D as t D~S w 2 xy T .
4D for a = 1. Anomalous subdiffusion can occur due to macromolecular crowding [32]. As fluorescent tracer molecules applied in this study are small, we did not expect strong deviations in the anomality parameter, which is supported by Figure S3. All FCS measurements were performed at RT (21uC).

Time-lapse Diffusion Measurements
Intrinsic diffusion rates were measured with an optically transparent diffusion setup in which fluorescent tracer movement was imaged using a confocal fluorescence microscope (FV1000, Olympus). Beetle luciferin (MW: 318 g?mol 21 , Promega) was used as a tracer molecule. Luciferin is fluorescent with an emission maximum at 537 nm and absorption maximum near 328 nm in acidic solutions and near 384 nm in basic solutions [33]. The emission maximum of the oxyluciferin complex is 523 nm [33,34]. Tracers were excited with a 375 nm laser line. Imaging was performed with a DM375-405/515/635 primary beam splitter in combination with a BA535-565 emission filter for visualization of luciferin diffusion and a BA505-540 emission filter for detection of accumulated oxyluciferin complex. Tracer motion was dependent on a concentration gradient induced between a saturated tracer concentration (100 mM) in DMEM and a tracer-free 2% agarose gel. Agarose gels were produced in glass Pasteur pipette tips (ID, 1.05 mm; OD, 1.7 mm) and at the onset of the diffusion experiment were connected to transparent silicon tubing (ID, 1.57 mm; Cole-Parmer) containing the saturated tracer solution. Image sequences were acquired with 5 min intervals during a 3 hour scanning period. An area of 19206640 pixels was visualized with a pixel size of 1.4261.42 mm. Imaging was performed with a UPLSAPO 106 air objective (NA: 0.40), focused on the middle plane of the agarose gel, and with a pinhole size of 400 mm. Measurements were performed at 37uC. Empty agarose gels were imaged to correct for background intensities. Intensity profiles acquired during the tracer diffusion experiment were normalized to the average intensity measurements obtained from a tracersaturated agarose gel. Image sequences were processed in Fiji (NIH, Bethesda, MD, USA). Luciferin diffusion rates were obtained by least squares fitting an analytical solution of Fick's second diffusion law (diffusion in semi-infinite media [35]) to the resulting averaged axial intensity profiles, Where, C is the concentration of tracer diffusing into the agarose gel, C 0 is the initial tracer concentration in the tracer-saturated agarose gel, and D is the diffusion coefficient obtained by profile fitting. The profile fitting algorithm was implemented in MATLAB. The initial saturated tracer concentration (in DMEM) contained in the silicon tubing was assumed to remain constant during the diffusion experiment. Validity of this assumption was verified by comparing the analytical solution of luciferin diffusion with a numerical solution of the diffusion problem implemented in COMSOL (COMSOL Multiphysics, Burlington, MA, USA).

OSB Production and Calibration
OSB are microscale oxygen sensors that consist of a core material, on which an oxygen sensitive and an oxygen insensitive dye are deposited, and a shell material, that prevents the dyes from leaching and protects the cells against potential dye toxicity. The protocol for OSB production is based on a strategy we previously developed to produce avidin-coated OSB [28] and is composed of two steps.
The core material of the OSB consists of silica gel (spherical silica gel, 45809, Alfa Aesar) and has a diameter of 5 mm. An initial amount of 50 mg silica gel spheres were brought into suspension with a 1 ml aqueous NaOH (0.01 N) solution and were mechanically stirred at 250 rpm for 30 min. Oxygen measurements relied on the incorporation of two fluorescent dyes; an oxygen sensitive Tris (4,7-diphenyl-1,10-phenanthroline) ruthenium (II) dichloride (76886, Sigma) and an oxygen insensitive reference fluorophore, Rhodamine 6G (56226, Fluka). Two 250 ml solutions were prepared containing these molecules in concentrations of 500 and 50 mM, respectively. These solutions were poured into the silica gel solutions and mechanically stirred (1 hour, 250 rpm). Afterwards the beads were diluted with 4 ml Milli-Q water and washed and centrifuged (13006g, 20 min) three times with Milli-Q water and once in ethanol. Beads were resuspended in 6.5 ml ethanol. Coating of silica gel spheres was performed via a Stöber seed growth method of silica shells described by the van Blaaderen group [36]. In brief, a solution of ammonia in water (25 w/v % NH 4 OH in Milli-Q water, 250 ml) and tetraethyl orthosilicate (99% TEOS, 32.5 ml, 86578, Sigma) was added to the silica spheres suspended in ethanol. Next this solution was mechanically stirred (250 rpm) for 18 hours. Suspended beads were collected by centrifugation (13006g, 20 min) and washed three times in PBS solution (5 ml). The bead solution was autoclaved and supernatant removed after centrifugation (13006g, 10 min). Beads were incubated at low temperature (4uC) in 2 ml DMEM and stored in the dark at a density of ,3610 8 beads?ml 21 . The amount TEOS solution to be added to the Stöber seed growth reaction, was determined from the relation [37,38]: where L s is the initial size of particles before shell growth, L f is the final size after shell growth (to obtain ,400 nm silica shell thickness), W s is the weight of SiO 2 initially present as seeds, W is the amount of SiO 2 added to the reaction mixture in the form of TEOS, and k r = (r 0 /r) 1/3 is a correction factor accounting for the differences in density between seed and final particles [39], with r 0 the density of the initial particle and r the density of the final particle. OSB were calibrated by exposure to hypoxic (0% O 2 via mixing with 1% Na 2 SO 3 ) and normoxic (21% O 2 ) oxygen conditions in DMEM. Microbeads were imaged on an inverted confocal laser scanning microscope (FV1000, Olympus) with a UPLSAPO 206 air objective (NA: 0.75). Excitation of the fluorescent beads was performed with a 488 nm Ar laser. Fluorescence emission was detected with a DM405/488 primary beam splitter in combination with a BA505-540 emission filter for detection of the Rhodamine 6G and a BA575-675 emission filter for detection of the ruthenium complex. Image stacks were acquired along the radial direction and consisted of 6406640 pixels with a pixel size of 0.16560.165 mm and images (n = 5) acquired in z-direction were separated by a distance of 50 mm. Experiments were performed at 37uC incubation temperature (Solent Scientific).
Regions of interest (ROI) were defined around individual beads for calibration and measurement of OSB. Images were processed in Fiji with a threshold macro. Average signal intensities were subsequently inserted into the Stern-Volmer equation. Further details on the calibration procedure can be found elsewhere [28].

Focused Ion Beam Milling and Scanning Electron Microscopy
OSB were placed on SEM specimen holders and coated with a thin conducting platinum layer by a sputter coater system (Quorum Q150 TS, Quorum Technologies Ltd, Laughton, UK). FIB ablation and SEM imaging were performed on a Nova 600 Nanolab dual-beam system (DualBeam TM -SEM/FIB, FEI, Hillsboro, USA). Beams are separated from each other by an angle of 52u. Therefore, vertical cutting of the OSB with the FIB proceeded by tilting the stage with this angle and vertical crosssections with the electron beam were also observed at this angular position. The SEM images were generated in high vacuum mode using an acceleration voltage of 5 keV and an Everhart Thornley detector.

Assembly of Setup for Validation of PDE-based Model
Low melting point agarose gels (2% in DMEM, Invitrogen) containing 1610 6 cells?ml 21 and supplemented with 2.5610 6 OSB?ml 21 , were prepared between two glass coverslips. The slides were separated by a 2 mm thick Teflon spacer, which consisted of two parts and at the central position contained a cylindrical hole (diameter, 8 mm) that served as a mold for the gel production. After the gel was poured into this hole, the second coverslip was quickly placed on top of the spacer and gelation was continued for 5 min. Upon completion of the gelation process, both parts of the mold were carefully removed and the cell-seeded hydrogel was transferred to a tissue culture dish (diameter, 35 mm, Greiner). Hydrogels were submerged in 2 ml DMEM and incubated at 37uC in a humidified atmosphere containing 5% CO 2 .
Bioluminescence intensity measurements. Cultivation media were replaced with DMEM supplemented with 47 mM luciferin. The validation setup was transferred to the imaging platform of an in vivo imaging system (IVIS 100, Perkin-Elmer, USA). Images were taken with a 1-inch CCD camera cooled to 2105uC. The field of view (FOV) was set to 10610 cm. Image sequences (acquisition time, 1 min) were acquired with 5 min intervals during a scanning period of 3 hours, 10 min intervals during 5 hours, and 30 min intervals during 4 hours. Images were processed in the LivingImage software (Perkin-Elmer) and radiance units refer to the number of photons per second that are leaving a square centimeter of the h and radiating into a solid angle of one steradian (V, fraction of the isotropic radiation field which can be thought of as a 3D cone of light emitted from the surface). Measurements were performed at 37uC.
Single-cell bioluminescence microscopy. After luciferin substrate addition, the validation setup was transferred to the stage of luminescence microscope (LV 200, Olympus). Bioluminescence was imaged with a UPLSAPO 606 oil objective (NA: 1.35) and transmitted to a cooled CCD camera. Time-lapse images were collected with 5 s acquisition times and EM gain was set at 12006 (photon imaging mode, Hamamatsu) Active cells were imaged at the central hydrogel position within a confocal volume at the maximal working distance (i.e. 150 mm) from the glass coverslip.

Quantification of Active Cell Numbers in Agarose Gels
Viable cell numbers were quantified using the Live/Dead viability/cytotoxicity kit (Invitrogen). Agarose hydrogels were rinsed with PBS solution, covered with Live/Dead staining solution containing 2 mM calcein AM and 4 mM ethidium homodimer-1 in PBS, and incubated for 1 hour in the dark. The dye solution was discarded afterwards and background due to residual stain was washed away with PBS. Live and dead cell numbers were counted via a sphere fitting algorithm in Imaris Bitplane (Zurich, Switzerland). DNA content was evaluated according to a protocol described by Grayson et al. [40] Briefly, agarose hydrogels were washed in PBS, transferred to 400 ml of digestion buffer (10 mM Tris, 1 mM EDTA, and 0.1% Triton X-100) with 0.1 mg/ml proteinase K in centrifugation tubes, and incubated at 56uC for 3 hours. Supernatant was collected after removal of debris by centrifugation (13,000 rpm, 1 min) and measured with a Qubit system (Invitrogen).

Visualization of Cell Mitochondria
293T cells were seeded on PLL coated glass coverslips (200,000 cells per slide). Cells were incubated overnight prior to staining. Cell mitochondria were stained with a 200 nM MitoTracker dye concentration (MitoTracker Red CM-H2Xros, Invitrogen) in DMEM for 20 min in the dark. The staining solution was supplemented with Na 2 SO 3 (1%) to expose the cells to hypoxia. Cell nuclei were stained by addition of the cell permeable nucleic acid stain, Hoechst (Invitrogen). After the staining was completed, adherent cells were washed with pre-warmed DMEM. Cells were subsequently fixed by medium replacement with DMEM containing 2% formaldehyde at 37uC (30 min). After fixation, the cells were rinsed with PBS and stored in PBS for visualization. Cell mitochondria were visualized on an inverted confocal microscope (FV1000, Olympus) with a UPLSAPO 1006 oil objective (NA: 1.4). Image stacks were acquired with a pixel size of 73673 nm and consisted of 102461024 pixels. . Reference values for Rhodamine 6G tracer diffusion in DMEM are also shown. The experimentally obtained autocorrelation curves were best fitted with a function that accounts for anomalous subdiffusion, presumably induced by macromolecular crowding [32,48]. (B) Anomality parameter for tracer diffusion in agarose or in DMEM. As expected from the small size of fluorescent tracer molecules and the low agarose concentration, we did not observe significant changes in anomality. All measurements were performed at RT (21uC). Error bars, 61 s.d. unit; n = 5.  Figure S8 (A) Fluorescence imaging of a luciferin concentration gradient established between a saturated solution (left) and a tracer-free 293T cell-seeded agarose gel (right). Positions occupied by the cells appear as dark spots (red arrowheads). (B) Multichannel fluorescence image shows cells embedded in the agarose (green) in combination with diffusing luciferin (red) and reacted oxyluciferin (blue). Scale bar, 400 mm.

(TIF)
Table S1 Overview of the parameter values implemented in the bioluminescence-oxygen model.

(DOCX)
Video S1 2D bioluminescence profiles of cell-seeded agarose gels imaged from top position (hydrogel diameter, 8 mm). Profiles are imaged with an IVIS 100 imaging system. Radiance units are in p s 21 cm 22 . (AVI)