Hypoxia-induced drug-resistance bias 3D cancer spheroid drug screens

TL;DR

This study shows that hypoxia in 3D cancer spheroids causes drug resistance, and offers a framework to select effective drugs for all cell types in spheroids.

Contribution

The study identifies hypoxia-induced quiescence as a key barrier to drug efficacy in 3D spheroids and proposes a method to overcome it.

Findings

Oxygen decay in spheroids causes hypoxia beyond 20 μm from the surface.

Spheroids larger than 150 μm cannot achieve IC50 for drugs targeting proliferating cells.

Targeting viability barriers enables effective drug selection for all cell populations in spheroids.

Abstract

Cellular 3D cancer spheroid technologies are novel tools that facilitate large-scale drug screening to bridge the in vitro–in vivo gap, without the cross-species effects of animal models. However, many spheroid studies fail to achieve
IC50 (dosage for 50% inhibition) even for unreasonably high applied drug concentrations (up to 1000× 2D
IC50). By mapping oxygen transport in patient-derived pancreatic cancer spheroids, this limiting viability is attributed to a near-universal oxygen decay gradient that renders cells deeper than 20 μm from the spheroid surface hypoxically quiescent and resistant to many chemotherapeutic drugs. The dose-independent viability barrier prevents
IC50 from being achieved for spheroids larger than 150 μm in diameter if the applied drug is dependent on the proliferating cell behavior. By examining three cancer cell types and five chemotherapeutic drugs, targeting this limiting viability barrier allows the selection of drugs and adjuvants that are effective in treating all cell populations within a spheroid. The reported analysis provides a framework for the accurate assessment of drug efficacy to target both well-oxygenated proliferating cells and hypoxically quiescent cells in biologically relevant and realistic 3D spheroid systems.

Figures

Dual drug–oxygen reaction–diffusion gradient model. (a) HepG2 spheroids cultured in increasing concentrations of PTX (top to bottom) stained with CellTracker™ Blue (blue) and BOBO-3 (red). Scale bars represent 50 μm. (b) PANC-1 spheroids of various sizes stained with Acridine Orange (green) and Propidium Iodide (red). Scale bars represent 50 μm. (c) Schematic of regions present in the dual reaction–diffusion gradient model featuring resultant layers caused by transport of applied therapeutic drugs (left) and oxygen (right). The center displays the two key dimensionless Thiele Moduli ( ΦDrug2 and ΦO22) and their key characteristic reaction–diffusion constants [ (k/D)Drug and (k/D)O2].

Reaction–diffusion model explains dose–response behaviors at low overdosage values. (a)–(c) Representative profiles of HepG2 spheroids cultured in increasing concentrations of PTX (left to right) stained with CellTracker™ Blue (blue) and BOBO-3 (red). (a) Multiphoton confocal images of spheroids with the top row showing the full z-stack averaged composite micrograph and second row showing the central cross section z-stack slice of the spheroid. Scale bars represent 50 μm. (b) Radial profiles of CellTracker™ Blue (blue) and BOBO-3 (red) for each spheroid. All profiles are plotted against normalized radial coordinates with fluorescent intensities in absorbance units. (c) Calculated radial viability distributions (black) for each spheroid plotted against expected local overdosage (purple) based on fitted (k/D)Drug=0.15 μm−2. Dotted purple line represents ODr=r*=1 (or Cr=r*=CTox), and shaded purple region represents region where local drug concentrations exceed the value of 2D IC50. Blue shaded region represents fully oxygenated region. (d) Dose–response curve plotted against fold-increase overdosages. Boxes and whiskers display quartile ranges for each treatment group (n = 25 for each group). Black lines represent theoretical viability fractions given the reaction–diffusion model per Eqs. (2) and (4), and spheroid size distribution. (e) Dose–response curve plotted against drug overdosage values normalized to drug Thiele moduli per Eq. (6). Theoretical curve represents expected dose–response behavior given high Φ2 reaction–diffusion models. Error bars represent one standard deviation of the mean.

Spatial oxygen concentration profiles are estimated using projection profile fitting. (a) Procedure for mapping measured fluorescent images of PANC-1 spheroids stained with BioTracker™ 520 Green Hypoxia Dye to a dataset of prepared projections for an exhaustive set of expected radial profiles of spheroids of variable size and ΦO22 accounting for light attenuation through tissue. (b) Representative images of select overlay micrographs of PANC-1 spheroids of various sizes in the dark-field (white) and stained with BioTracker 520 Green Hypoxia Dye (green). Scale bars represent 100 μm. (c) Fitted (k/D)O2 values for PANC-1 cell line spheroids (n = 26) and patient-derived pancreatic cancer spheroids (n = 56), displayed with literature-derived (k/D)O2 for HepG2 in hydrogel scaled to cell density comparable to spheroids (1 cell per 1000 μm3).34,35

Oxygen gradient toxicity limit. (a) Diagram of radial regions within spheroid produced by local oxygen concentrations. Starting from the surface of the spheroid, these are a fully viable well-oxygenated outer shell ( r∈(RQ,R]), a mantle region of viable cells with subpar oxygen rendered quiescent ( r∈(R2,RQ]), a secondary shell partially dying from lack of oxygen ( r∈(R1,R2]), and an inner core of fully dead, anoxic cells ( r∈[0,R1]). (b) Fluorescent micrographs of PANC-1 cell line spheroids for three representative spheroids. Scale bars for micrographs represent 50 μm. (c) Simulated radial profiles (left) and light-attenuated projections (right) of live stain for the three imaged spheroids. (d) Radial fluorescent profiles (black line) and optimally fitted projection profiles (red line) for the three imaged spheroids.

Fitted oxygen profiles predict hypoxic core cell death in large spheroids. Radial distribution plot for PANC-1 spheroids ( n=37) with diameters ranging from 61 to 436 μm. The surface of the spheroid is represented by the line of identity with unit-slope, with the regions below representing the space within the spheroids. Theoretical local oxygen concentrations [calculated from (k/D)O2,PANC−1 in Fig. 3(c)] delineate regions of full oxygenation (blue), quiescently low oxygenation (green), and lethally low oxygenation (red). The green and red lines represent the threshold radii where local CO2(r) = 100 and 5 μM (10% and 0.5%), respectively, given (k/D)O2=1.82×10−3 μm−2. (R1,R2)-tuple pairs are plotted for each spheroid in accordance with their sizes as crosses and open circles, respectively. The dotted lines represent best fit for non-zero values of the R1 and R2 series, and the gray region represents the expected range of the dying layer. The red arrow at the top right corner represents the layer of well-oxygenated proliferating cells that are not resistant to any applied drug.

Many cell–drug spheroid systems do not experience a drug dose–response below a thin oxygenated layer (red arrow). (a) ΦDrug2-normalized dose–response curves for (i) large LNCaP (blue) and HepG2 (red) spheroids of similar sizes and (ii) small LNCaP (green) spheroids. LNCaP spheroids were treated with Doc, as adapted from results featured Ref. 13. HepG2 spheroids were treated with PTX. (b) Adapted Φ2-normalized dose–response curves from Ref. 15 for W1 ovarian cancer cell line spheroids treated using (i) DNA inhibitors (blue), microtubule inhibitors (red), and (ii) topoisomerase inhibitors (green). Error bars represent one standard deviation. The red arrow in each figure corresponds to the same red arrow in Fig. 5 representing the population of proliferating cells in the well-oxygenated surface of the spheroid. These cells are not resistant to any applied drugs. As a result, the arrow also corresponds to the limiting (asymptotic) net viability of a spheroid if the drug can be resisted through hypoxic quiescence. This limiting viability is demonstrated at high drug dosages for all cell–drug combinations in each system within a(i, ii) and for the four drugs in b(i). Yet, DOX and TOP treated W1 cells in b(ii) surpass this limit.

Universal viability limit of drug dosage efficacy bounded by spheroid oxygenation. (a) and (b) Shaded zones represent the fraction of cells that are well oxygenated (blue) and hypoxically quiescent (red) for a given spheroid size. (a) The universal curve (red dashed curve) denoting expected viability plateaus for a given spheroid size represented by the dimensionless value of ρ50 given (k/D)O2=1.82±0.67×10−3 μm−2. Experimental data for nine cell–drug combinations using spheroids of various sizes are shown.13,15 Error bars represent one standard deviation. (b) Predicted oxygenated shell thickness for each spheroid system calculated using dose–response plateaus. Error bars represent the standard errors of the mean. In the large W1 spheroids with R>RCrit, only DOX and TOP could achieve IC50 and yield FPlat<0.50 at OD>100 [green data points in (a)].

Full text

INTRODUCTION

I.

In 2025, over 2 million new cancer cases are expected to occur.1 In order to combat this, a broader tool-belt of therapeutic options must be explored.2–4 Yet, development of novel therapeutics is stymied by the slow, often ineffective drug screening methodologies currently in place, rendering less than 1% of drug candidates to pass clinical testing.5,6 For classical small-molecule drug screening, standardized efficacy is reported using the half-maximal inhibitory concentration ( $[eqn]$ ).7–9 Over the years, various complications such as drug resistance have called into question the use of $[eqn]$ as a standard metric since it often fails to accurately capture response to dose escalation.10–12 Particularly, the impetus to introduce monoculture or coculture cellular 3D spheroids to attain more physiologically relevant responses has resulted in systems where previously effective treatments can only achieve at most 50% of their prior efficacy even using elevated dosages several folds above toxicity.13–15 Earlier studies have demonstrated that transport limitations of small-molecule drugs reduce their efficacy in 3D systems.16,17 This effect is further exacerbated for applications of highly anticipated antibody-based immunotherapies, leading to reduced clinical implementation in the field.18 To resolve these limitations, reaction–diffusion models have been introduced to simulate transport phenomena within spheroid systems.19,20 Many of these studies focus on ameliorating nutrient and waste gradients in order to maintain long-term 3D cultures as drug screening tools. However, studies that scrutinize only drug transport into tumors often do not capture the dose–response behaviors seen in screening studies.21

Hypoxic conditions within the tumor microenvironment have been shown to exhibit resistance to applied drug and radiation therapies.22–24 Qualitative models of multi-shell spheroids based on local oxygen concentrations have been proposed as far back as 1979, detailing a proliferating outer shell, a non-cycling intermediary region, and a necrotic core.15,25 Local oxygen concentrations from 62.5 to 100 μM (6.25%–10% atmospheric $[eqn]$ ) have been reported to produce hypoxia-inducible factor (HIF-1α) and multi-drug resistance (MDR1) leading to downstream escape from the cell cycle and subsequent drug resistance.26–29 Severe hypoxia at levels of below 5 μM (0.5%) oxygen (known as anoxia) has been shown to be lethal to cells, causing apoptosis.30 Some have reported that hypoxia directly induces intrinsic drug resistance through altered gene expression, where cells from hypoxic and anoxic layers more closely reflect the gene expression landscape of in vivo tumors.31 Despite mounting interest in drug–tumor interactions, limited work has been done to rigorously interrogate the implications of drug transport gradients in concert with other gradient systems. In this work, we establish a dual-gradient model that pairs reaction–diffusion gradients for the transport of applied small-molecule therapeutics and oxygen (Fig. 1). Using image processing of fluorescence micrographs, we estimate and report several transport coefficients for various drug–cell combinations. This work proposes a framework for the development of new metrics in order to better predict dose–response behavior for spheroids of various length scales and drug doses in the presence of a dominant oxygen gradient. We highlight the importance of selecting and maintaining consistent spheroid sizes, as dual-gradient behavior varies with high sensitivity at different length scales. However, if size uniformity is not achievable, we offer up a normalization protocol to produce a metric that can be used to compare spheroids of all length scales.

Dual drug–oxygen reaction–diffusion gradient model. (a) HepG2 spheroids cultured in increasing concentrations of PTX (top to bottom) stained with CellTracker™ Blue (blue) and BOBO-3 (red). Scale bars represent 50 μm. (b) PANC-1 spheroids of various sizes stained with Acridine Orange (green) and Propidium Iodide (red). Scale bars represent 50 μm. (c) Schematic of regions present in the dual reaction–diffusion gradient model featuring resultant layers caused by transport of applied therapeutic drugs (left) and oxygen (right). The center displays the two key dimensionless Thiele Moduli ( ΦDrug2 and ΦO22) and their key characteristic reaction–diffusion constants [ (k/D)Drug and (k/D)O2].

RESULTS

II.

Dual-gradient reaction–diffusion spheroid model

A.

To characterize concentration gradients, we model the cellular spheroid as a homogeneous, steady-state spherical catalytic particle of radius (R) subjected to mass transfer via diffusion and consumption via cellular uptake.32 Scaling the radial coordinate by the spheroid radius (R) and the concentration distribution by the applied bulk concentration ( $[eqn]$ ), one obtains the dimensionless equation

$$[eqn]$$

where $[eqn]$ is the Thiele modulus of the species i in the model spheroid that defines the ratio between the species uptake rate (k) and mass-transfer rate ( $[eqn]$ ). In the cases of the drug and oxygen transport gradients, the key factors to $[eqn]$ are highlighted in Fig. 1(c). The solution to Eq. (1) has the explicit form

$$[eqn]$$

For large $[eqn]$ , the majority of species i within the spheroid would be confined to the thin boundary layer near the spheroid surface ( $[eqn]$ ) and Eq. (2) can be approximated by

$$[eqn]$$

We note this asymptote is independent of the spheroid radius (R) and corresponds to a flat tissue layer. To highlight this behavior, consider the reaction–diffusion system of oxygen into spheroids. For large spheroids ( $[eqn]$ ), it suggests an oxygen decay length of $[eqn]$ over which the oxygen concentration decays by a factor of e per Eq. (3). If the applied bulk concentration is fixed at the standard environmental value of 200 μM at atmospheric 20% partial pressure conditions ( $[eqn]$ ) and the critical oxygen level for quiescence is a fixed fraction of that bulk value, a universal thickness can be estimated as a factor of $[eqn]$ for the fully oxygenated cells outside the hypoxic layer. If the drug kills only the fully oxygenated cells, an asymptotic viability fraction exists that cannot be reduced regardless of applied bulk drug concentrations ( $[eqn]$ ). If this asymptotic viability fraction is above 50%, it is by definition impossible to kill half of the spheroid cells and realize an $[eqn]$ .

This prediction is consistent with micrographs in Fig. 1(a). Increasing dosages of paclitaxel (PTX) treatment does not change the thickness of the outer shell of dead HepG2 cells as stained by BOBO-3 (red), which is not readily able to permeate intact cell membranes. The inner dead cells due to hypoxia in a spheroid of comparable size [bottom image of Fig. 1(b)] do not overlap with the outer dead HepG2 cells due to PTX treatment even under the highest dosage. This suggests an intermediary layer of hypoxic cells that are viable and drug-resistant. Note that the cores of Fig. 1(a) (stained with CellTracker™ Blue) do not contain significant amounts of dead cells.

For smaller, well-oxygenated spheroids ( $[eqn]$ ), hypoxically quiescent cells do not exist, allowing for the drug gradient to solely dictate cell death in the spheroid. To model the dose–response of a small spheroid system, we define a dimensionless radius $[eqn]$ within a spheroid representing the radial coordinate above which the local drug concentration is greater than the known $[eqn]$ for the same cells within a 2D culture $[eqn]$ . The overdosage (OD) is a dimensionless quantity representing the ratio of bulk drug concentration to the $[eqn]$ of the drug applied to the 2D culture. Note that, for PTX on HepG2, 2D $[eqn]$ was set as 0.856 μM from 2D cell-culture dose–response studies.33 The value of $[eqn]$ can be estimated from Eq. (2), where the transported species is PTX. Without hypoxically induced drug resistance, the viability fraction (F) of the spheroid is only a function of $[eqn]$ as the ratio of the volumes of viable core to the total spheroid,

$$[eqn]$$

This result can be made explicit for the case of $[eqn]$ by using the analog of Eq. (3) where the species is the applied drug,

$$[eqn]$$ $$[eqn]$$

In Fig. 2(d), expected drug response profiles are fitted to experimental dose–response data for HepG2 cells treated with PTX. In Fig. 6(a), the same is done for dose–response data for LNCaP cells treated with dosages of Docetaxel (Doc) from Ref. 13. The $[eqn]$ for 2D cultures of LNCaP cells treated with Doc was set as 2.45 nM.13 Curves were fitted using reported spheroid size distributions, contributing to the families of expected curves. Each cell–drug testing pair was fitted with curves yielding estimated $[eqn]$ values of 0.15 μm^−2^ for HepG2-PTX and 0.05 μm^−2^ for LNCaP-Doc. These correspond to decay lengths $[eqn]$ of 2.6 and 4.5 μm, respectively. Over these lengths, the local drug concentration decays by a factor of e of the applied bulk concentration ( $[eqn]$ ). This small decay length suggests that the OD ratio must exceed 100 to raise the core concentration beyond the toxic level even for 20 μm diameter spheroids, the smallest size acceptable to be considered a fully formed cellular spheroid. For reasonably sized spheroids, $[eqn]$ is required to achieve an $[eqn]$ , if at all, as seen in Fig. 2. As expected, Eq. (6) quantitatively captures all viability data in Fig. 6 for larger spheroids ( $[eqn]$ ) for the 20 μm region nearest to the surface.

Reaction–diffusion model explains dose–response behaviors at low overdosage values. (a)–(c) Representative profiles of HepG2 spheroids cultured in increasing concentrations of PTX (left to right) stained with CellTracker™ Blue (blue) and BOBO-3 (red). (a) Multiphoton confocal images of spheroids with the top row showing the full z-stack averaged composite micrograph and second row showing the central cross section z-stack slice of the spheroid. Scale bars represent 50 μm. (b) Radial profiles of CellTracker™ Blue (blue) and BOBO-3 (red) for each spheroid. All profiles are plotted against normalized radial coordinates with fluorescent intensities in absorbance units. (c) Calculated radial viability distributions (black) for each spheroid plotted against expected local overdosage (purple) based on fitted (k/D)Drug=0.15 μm−2. Dotted purple line represents ODr=r=1 (or Cr=r*=CTox), and shaded purple region represents region where local drug concentrations exceed the value of 2D IC50. Blue shaded region represents fully oxygenated region. (d) Dose–response curve plotted against fold-increase overdosages. Boxes and whiskers display quartile ranges for each treatment group (n = 25 for each group). Black lines represent theoretical viability fractions given the reaction–diffusion model per Eqs. (2) and (4), and spheroid size distribution. (e) Dose–response curve plotted against drug overdosage values normalized to drug Thiele moduli per Eq. (6). Theoretical curve represents expected dose–response behavior given high Φ2 reaction–diffusion models. Error bars represent one standard deviation of the mean.*

Larger spheroids ( $[eqn]$ ) require exponentially higher OD values of $[eqn]$ and $[eqn]$ to penetrate and inhibit the outer one to two layers of cells, respectively. Astronomically higher dosages relative to those required to achieve $[eqn]$ for 2D cell cultures are hence required to achieve similar inhibition effects in 3D cellular spheroids. These dosages are impossible to administer to patients safely in practical applications. The implication is that large spheroids ( $[eqn]$ ) can never achieve the 50% inhibition necessary to derive a traditional $[eqn]$ . Moreover, from our analysis of the dose–response data in Fig. 6(a), $[eqn]$ values cannot be calculated for spheroids larger than 150 μm in diameter ( $[eqn]$ ) even for OD values exceeding 100. As seen in Fig. 6(a), the relationship presented in Eq. (5) begins to break down at reasonable OD values even for larger spheroids ( $[eqn]$ ). This behavior suggests that, beyond exponential decay of the radial drug concentration profiles, the coexisting oxygen transport profile is causing drug resistance in subsurface cell layers due to oxygen deprivation. In order to determine the thickness of the fully viable, drug-resistant, quiescent cell region, we must resolve the radial oxygen profile and estimate the associated $[eqn]$ , corresponding to the radius R in Fig. 3 when $[eqn]$ .

Spatial oxygen concentration profiles are estimated using projection profile fitting. (a) Procedure for mapping measured fluorescent images of PANC-1 spheroids stained with BioTracker™ 520 Green Hypoxia Dye to a dataset of prepared projections for an exhaustive set of expected radial profiles of spheroids of variable size and ΦO22 accounting for light attenuation through tissue. (b) Representative images of select overlay micrographs of PANC-1 spheroids of various sizes in the dark-field (white) and stained with BioTracker 520 Green Hypoxia Dye (green). Scale bars represent 100 μm. (c) Fitted (k/D)O2 values for PANC-1 cell line spheroids (n = 26) and patient-derived pancreatic cancer spheroids (n = 56), displayed with literature-derived (k/D)O2 for HepG2 in hydrogel scaled to cell density comparable to spheroids (1 cell per 1000 μm3).34,35

Oxygen transport behavior

B.

To estimate $[eqn]$ for different cell lines, cellular spheroids generated from the PANC-1 pancreatic cancer cell line and 3R-4/3R-5-1 patient-derived pancreatic cancer cells were stained using BioTracker™ 520 Green Hypoxia Dye and measured using fluorescence microscopy. The dye increases in fluorescence inversely to local oxygen concentrations. Assuming the concentration distribution described in Eq. (3), we would expect to see a radial profile where fluorescence is highest in the spheroid center and lowest near the surface.

The images captured via microscopy are 2D circular projections of the original radial-symmetric 3D profile. A workflow [Fig. 3(a)] was developed to generate expected projection profiles of variable $[eqn]$ for given sizes of spheroids. Light attenuation through the spheroid structure was accounted for using an adapted form of Beer–Lambert law [Eq. (7)]. The absorption and scattering coefficients for excitation and emission wavelengths of these spheroids are approximated to that of human dermis,36

$$[eqn]$$

The radial hypoxia intensity profiles and effective diameters were measured in the spheroid micrographs [Fig. 3(b)]. The profiles were subsequently normalized to both spheroid size and maximum intensities while retaining their curvature. By comparing the curvature of the normalized intensity profiles, one can approximate a Thiele modulus ( $[eqn]$ ) [see Fig. 3(a) for typical profiles at different Thiele moduli]. These experimental radial intensity profiles were fitted to a database of simulated projection curves to approximate $[eqn]$ for a given spheroid size. Using the definition of the Thiele modulus, we can approximate $[eqn]$ -values by noting that the $[eqn]$ is directly related to $[eqn]$ as seen in Eq. (8). By fitting curves of slope ( $[eqn]$ ), those of best fit would yield intercept values approximating this constant term,

$$[eqn]$$

Using the above approach, $[eqn]$ were approximated for the PANC-1 cell line spheroids (n = 26) and the patient-derived pancreatic cancer (3R-4/3R-5-1) spheroids (n = 56) [Fig. 3(c)]. These values correspond to an oxygen decay length $[eqn]$ of about 10–23 μm, approximating to 2× or 5× deeper penetration as compared to small-molecule drugs and lining up well with the expected range of $[eqn]$ in typical cancer cell lines and extracellular matrices (ECM) compositions. In naturally occurring concentrations, the diffusivities within these matrices do not exceed ranges of $[eqn]$ for Matrigel to $[eqn]$ for low-density collagen,37 with edge cases of $[eqn]$ for small intestinal submucosa and $[eqn]$ for human dermis.38 For cellular oxygen uptake rates, overall oxygen consumption across 19 human and murine cancer cell lines was found to be $[eqn]$  attomoles O_2_/s per cell with edge cases of $[eqn]$  attomoles O_2_/s per cell for HL60 (leukemia cell line) and $[eqn]$  attomoles O_2_/s per cell for HeLa (cervical cancer).39 Taking these ranges, we calculate possible oxygenation depths ranging from 3.5 to 18 μm, with the average landing around 8.5 μm. As such, the oxygenation behavior appears consistent between most cell lines with a fairly universal depth that is always about a single layer of cells (10–15 μm). Therefore, the oxygen decay length is near universal for various cell types and spheroid sizes. It remains to determine the critical oxygen level for cellular quiescence and necrosis.

Delineating hypoxia-derived regions

C.

The establishment of oxygen gradients within spheroids results in an explicit demarcation of spatial cell viability based on local oxygen concentrations. In Fig. 4(a), we define an untreated spheroid at an oxygen concentration of 20% (200 μM). Within this spheroid, four potential shell layers may form based on local oxygen concentrations: (i) fully oxygenated and proliferative cells, (ii) viable but hypoxically quiescent cells subject to cell cycle arrest or gene expression-linked drug resistance, (iii) cells subject to severely hypoxic conditions with reduced local viability, and (iv) a core of fully necrotic cells. These layers have been observed and documented in earlier studies of large cancer spheroids.24,40

Oxygen gradient toxicity limit. (a) Diagram of radial regions within spheroid produced by local oxygen concentrations. Starting from the surface of the spheroid, these are a fully viable well-oxygenated outer shell ( r∈(RQ,R]), a mantle region of viable cells with subpar oxygen rendered quiescent ( r∈(R2,RQ]), a secondary shell partially dying from lack of oxygen ( r∈(R1,R2]), and an inner core of fully dead, anoxic cells ( r∈[0,R1]). (b) Fluorescent micrographs of PANC-1 cell line spheroids for three representative spheroids. Scale bars for micrographs represent 50 μm. (c) Simulated radial profiles (left) and light-attenuated projections (right) of live stain for the three imaged spheroids. (d) Radial fluorescent profiles (black line) and optimally fitted projection profiles (red line) for the three imaged spheroids.

Using the PANC-1 cell line, spheroids were generated with diameters ranging from 50 to 450 μm. These spheroids were stained using Acridine Orange (AO) as a live cell dye. Based on our model, we would expect to see layers (i) and (ii) exhibit full fluorescence, while layer (iii) would exhibit reduced fluorescence, and layer (iv) would be dark. A stepwise function modeled this radial intensity distribution per Eq. (9), where the intermediary layer of (iii) is represented by a simple linear scaling between zero and one,

$$[eqn]$$

Applying the same approach from Fig. 3(a), a family of intensity projection curves was simulated for varying spheroid diameter, $[eqn]$ and $[eqn]$ . These two radial parameters establish the transitional region of cells in the process of dying between the hypoxic quiescent cells and necrotic core. By measuring and normalizing the radial intensity profiles of the AO-stained live cells within each spheroid [Fig. 4(b)], each spheroid was fitted to an associated simulated projection for a given spheroid size and $[eqn]$ -tuple [Figs. 4(c) and 4(d)]. We compiled the fitted $[eqn]$ -tuples of the PANC-1 spheroid catalog (n = 37) and plotted them against their measured sizes in Fig. 5. In accordance with the theoretical profile defined by Eq. (9), we visualize the expected transitional dying shell boundaries for all experimentally measured spheroids. We find that spheroids smaller than 290 μm in diameter (145 μm in radius) do not contain necrotic cores.

Fitted oxygen profiles predict hypoxic core cell death in large spheroids. Radial distribution plot for PANC-1 spheroids ( n=37) with diameters ranging from 61 to 436 μm. The surface of the spheroid is represented by the line of identity with unit-slope, with the regions below representing the space within the spheroids. Theoretical local oxygen concentrations [calculated from (k/D)O2,PANC−1 in Fig. 3(c)] delineate regions of full oxygenation (blue), quiescently low oxygenation (green), and lethally low oxygenation (red). The green and red lines represent the threshold radii where local CO2(r) = 100 and 5 μM (10% and 0.5%), respectively, given (k/D)O2=1.82×10−3 μm−2. (R1,R2)-tuple pairs are plotted for each spheroid in accordance with their sizes as crosses and open circles, respectively. The dotted lines represent best fit for non-zero values of the R1 and R2 series, and the gray region represents the expected range of the dying layer. The red arrow at the top right corner represents the layer of well-oxygenated proliferating cells that are not resistant to any applied drug.

Using the $[eqn]$ -value derived for PANC-1 cell line spheroids in Fig. 3(c), we can overlay curve traces representing the reported threshold concentrations for quiescence (green line) and lethality (red line) at 100 and 5 μM $[eqn]$ representing a reduction of the bulk oxygen concentration ( $[eqn]$ ) by a factor of 2 and 40, respectively. These traces establish three defined shaded regions that arise in different-sized spheroids: (i) a fully-oxygenated outer shell (blue zone), (ii) a region of hypoxically quiescent yet viable cells (green zone), and (iii) a hypoxically necrotic core (red zone). Likewise, based on these demarcations, three distinct sizes of spheroids can be defined: (i) small, well-oxygenated spheroids (<90 μm in diameter), (ii) mid-sized, partially quiescent spheroids (>90 and <290 μm in diameter), and (iii) large spheroids with necrotic cores (>290 μm in diameter). Per the simplified continuum model for large spheroids ( $[eqn]$ ) from Eq. (3) and the fitted $[eqn]$ for PANC-1, the approximate depth of the fully-oxygenated cell layer is a universal value of $[eqn]$ , or 15 μm, from the surface of the spheroid. Cell death is expected to start approximately $[eqn]$ , or 86 μm, from the surface.

Overlaying these regions in Fig. 5, we note that the zone of lethal oxygen concentrations lands well upon the fitted dying cell layer, where we expected to see spheroids exceeding 290 μm in diameter to exhibit a necrotic core. We do not expect spheroids smaller than 90 μm in diameter to exhibit quiescent cells in their composition. Key among all observations is that, for spheroids larger than 90 μm in diameter, only the initial layer of 15 μm from the surface of the spheroid is fully-oxygenated, normally-proliferating cells. As most chemotherapeutics treat cancer by targeting highly proliferative cells by intercepting cell division, these therapies only work effectively on spheroids smaller than 90 μm in diameter. For mid-sized and large spheroids, only the thin outer shell is susceptible to such therapeutic techniques.

Oxygenation-limited dose–response behavior

D.

Our establishment of spheroid oxygenation regions allows us to estimate the asymptotic viability fraction at arbitrarily large OD values. To verify that only the thin oxygenated surface of spheroids is susceptible to applied therapeutic treatments, we review the dose–response experiments once more. In Fig. 2(b), we radially measure the intensities of the CellTracker™ Blue and BOBO-3 stains for each HepG2 spheroid to attain the radial viability profiles of spheroids under different applied PTX dosages. In Fig. 2(c), with each fold-increase in PTX concentration, more cells within the given spheroid die as expected. However, only cells located within the 15-μm surface layer of the spheroid (indicated by the blue region and red arrows) are inhibited. Moreover, we observe that the final tenfold increase from 100- to 1000-μM yields no significant change in additional cell death.

Assuming a uniform drug diffusion from its surroundings, the net viability of the spheroid given a radial viability distribution can be calculated via Eq. (10). Plotting net viability for all 100 spheroids ( $[eqn]$ per treatment cohort), we derive the dose–response curve found in Fig. 2(d). The dose–response data follow the drug reaction–diffusion model (dark curves) at lower overdosages. The trend deviates when the expected viability fraction lands within a region representing the expected viability fraction at which only the well-oxygenated layer is inhibited (red arrow),

$$[eqn]$$

This minimum viability plateau is sensitively dependent on spheroid size, minimally dependent on the spheroid composition, and variably dependent on the drug application. Large (168 μm in diameter) LNCaP spheroids in Fig. 6(a, i) have a lower fraction of their total composition consisting of well-oxygenated cells as compared to small spheroids (84 μm in diameter) in Fig. 6(a, ii). This yields a smaller fraction of cells being susceptible to drug effects despite being exposed to typically lethal dosages of therapeutics. The size dependence is the primary determinant of the fraction of inhibitable cells within a spheroid. In Fig. 6(a, i), two different spheroid–drug systems with similar spheroid sizing converge to a similar minimum viability plateau. In terms of drug dependence, we note that the majority of drugs, such as DNA inhibitors (Methotrexate (MET) and Cisplatin (CIS)) and microtubule inhibitors (Vincristine (VIN) and Paclitaxel (PTX)) that target the disruption of cell division activities, are heavily susceptible to the hypoxia-induced oxygenation layer limit [Fig. 6(b, i)]. This occurs when hypoxically quiescent cells escape the cell cycle and fall dormant to preserve resources and maintain basal level metabolism.26–29 However, this is untrue for topoisomerase inhibitors (Doxorubicin (DOX) and Topotecan (TOP)) that are able to penetrate the spheroid past the well-oxygenated layer and show a consistent dose escalation response [Fig. 6(b, ii)] unlike all other cases with high significance (supplementary material Table 1).

Many cell–drug spheroid systems do not experience a drug dose–response below a thin oxygenated layer (red arrow). (a) ΦDrug2-normalized dose–response curves for (i) large LNCaP (blue) and HepG2 (red) spheroids of similar sizes and (ii) small LNCaP (green) spheroids. LNCaP spheroids were treated with Doc, as adapted from results featured Ref. 13. HepG2 spheroids were treated with PTX. (b) Adapted Φ2-normalized dose–response curves from Ref. 15 for W1 ovarian cancer cell line spheroids treated using (i) DNA inhibitors (blue), microtubule inhibitors (red), and (ii) topoisomerase inhibitors (green). Error bars represent one standard deviation. The red arrow in each figure corresponds to the same red arrow in Fig. 5 representing the population of proliferating cells in the well-oxygenated surface of the spheroid. These cells are not resistant to any applied drugs. As a result, the arrow also corresponds to the limiting (asymptotic) net viability of a spheroid if the drug can be resisted through hypoxic quiescence. This limiting viability is demonstrated at high drug dosages for all cell–drug combinations in each system within a(i, ii) and for the four drugs in b(i). Yet, DOX and TOP treated W1 cells in b(ii) surpass this limit.

As high OD viability plateaus develop, these asymptotic values depend only on the size-normalized oxygen penetration depth. Simply, viability plateaus become fundamentally independent of the drug profile. To predict this limit for all spheroid sizes, we define a new term $[eqn]$ as the ratio between the necessary oxygenation depth required to achieve 50% inhibition ( $[eqn]$ ) and the actual penetration depth of the oxygenated layer ( $[eqn]$ ). This oxygenated depth is defined as the spheroid radius at which the rates of oxygen uptake and mass transport are balanced ( $[eqn]$ ). In Sec. II C, we derived that in most cellular spheroids the $[eqn]$ corresponds to an oxygen decay length of 23 μm and a quiescent layer starting at a depth that is a factor $[eqn]$ of this length. Hence,

$$[eqn]$$

should be unity when $[eqn]$ and the asymptotic viability fraction are exactly 50% inhibited. All spheroids larger than this universal critical size ( $[eqn]$ ) will never achieve a dose–response at or below 50%. Here, the critical spheroid radius $[eqn]$ is 72 μm, corresponding to a spheroid diameter of about 150 μm. This value is reasonable as it falls into the expected range of possible $[eqn]$ values, i.e., 17–85 μm calculated from literature-derived k and D values for oxygen across several cell lines and ECM compositions.37–39 In Fig. 7(a), the estimated variability given this range is represented by the orange region surrounding the red theoretical curve. To find the theoretical relationship between $[eqn]$ and the viability plateau at high OD values ( $[eqn]$ ), we apply the geometric definition of dead-shell/live-core viability fraction employed earlier in Eq. (4). By defining $[eqn]$ and using our definitions of $[eqn]$ and $[eqn]$ , we derive the final relation of the theoretical viability plateau curve ( $[eqn]$ ) in relation to the $[eqn]$ -coordinate

$$[eqn]$$

Universal viability limit of drug dosage efficacy bounded by spheroid oxygenation. (a) and (b) Shaded zones represent the fraction of cells that are well oxygenated (blue) and hypoxically quiescent (red) for a given spheroid size. (a) The universal curve (red dashed curve) denoting expected viability plateaus for a given spheroid size represented by the dimensionless value of ρ50 given (k/D)O2=1.82±0.67×10−3 μm−2. Experimental data for nine cell–drug combinations using spheroids of various sizes are shown.13,15 Error bars represent one standard deviation. (b) Predicted oxygenated shell thickness for each spheroid system calculated using dose–response plateaus. Error bars represent the standard errors of the mean. In the large W1 spheroids with R>RCrit, only DOX and TOP could achieve IC50 and yield FPlat<0.50 at OD>100 [green data points in (a)].

In Fig. 7, the red curve illustrates the expected viability plateau for the range of possible $[eqn]$ -values plotted alongside the viability plateaus of several cellular spheroid systems (LNCaP, HepG2, and W1). These sphere systems are tested against a variety of chemotherapeutic drugs, including DNA inhibitors (CIS, DOX, and MTX), microtubule inhibitors (PTX and VIN), and a topoisomerase inhibitor (TOP). Given the universality of $[eqn]$ , the asymptotic viability fraction is only a function of spheroid size R for large spheroids ( $[eqn]$ ). This simple geometric relationship is consistent with data for all three cell lines and five drug conditions when $[eqn]$ cannot be realized ( $[eqn]$ ) where both $[eqn]$ and $[eqn]$ .

DISCUSSION

III.

This study demonstrates that looking exclusively at drug transport is insufficient to capture expected drug responses within 3D spheroid models, due to competing gradients such as oxygen transport. Our results highlight the importance of tight spheroid size control in dose–response studies, as discrepancies in spheroid length scales significantly alter local drug and oxygen concentrations available to cells throughout the spheroid. These gradients yield drastic incongruities in local cellular environments, leading to noticeably different viability profiles across spheroids of various sizes. Spheroids larger than 290 μm in diameter exhibit anoxic, necrotic cores leading to bias from drug-independent cell death. Spheroids larger than 90 μm in diameter exhibit hypoxic, quiescent cores leading to bias from hypoxia-induced drug resistance. Our analysis underscores how large cancer spheroids fail to fully recapitulate in situ tumors, which compensate in response to hypoxic conditions by activating angiogenesis and migration.41,42 Invasion by tumor-associated immune cells43,44 and blood/lymph angiogenesis into the tumor45–48 strongly impacts tumor microenvironment (TME) composition and behavior. The added complexities would greatly alter the material gradients and drug efficacy within the tumor in ways beyond the scope of this project. New models, both experimental and theoretical, must be developed with intentionality and control to further improve the drug discovery workflow. These behaviors cannot be fully reproduced by in vitro spheroid models. Recent works49,50 to develop physiologically representative vasculature in tumor organoid studies provide promising steps forward to bypass these limitations.

If the scope of a drug screening project focuses solely on analyzing drug transport behaviors, we suggest limiting spheroids to below 90 μm in diameter. For applications where hypoxia-induced drug resistance is part of the investigation, spheroids between 90 and 290 μm should be chosen. Our model provides the framework by which drug response may be evaluated under variable oxygen conditions. For instance, in Fig. 7, all cell–drug combination series exhibit asymptotic viability plateaus except for the W1 ovarian cancer spheroids treated with the topoisomerase inhibitors (DOX and TOP). While the other drugs primarily target cell division, these drugs disrupt DNA transcription, crippling both proliferating and quiescent cells, thus breaking the oxygenation barrier. An alternate strategy would be to oxygenate the entire tumor. Oxygen delivery to hypoxic regions in tumors has become a focus for improving treatment response in oxygen-deprived tissue.51,52 Furthermore, our normalized reaction–diffusion model captures the short-term (<24 h) hypoxia-induced drug resistance for each of our examined cell lines and drugs. However, we have noted a long-term (7 day) phenomenon of shell shedding that we have only observed in SP-2 (PDAC) spheroids treated with Gemcitabine (see supplementary material Figs. 3 and 4). Tumor spheroids have been reported to exhibit continuous tumor cell and fragment shedding as a part of the behavior of the TME niche.53 This process can be exacerbated in the presence of chemotherapeutic drug application. Treatment with drugs such as PTX has been shown to modify the TME in a pro-invasion fashion by altering ECM composition and signaling through surface receptors.54 Furthermore, chemotherapy drugs induce tumor release of enzymatic EVs that degrade the structural integrity of the tumor ECM.55 Gemcitabine specifically has been shown to directly cause epithelial–mesenchymal transition (EMT)-like changes in pancreatic cancer models.56 These behaviors are likely linked to the findings that show cell apoptosis events cause the release of pro-metastatic, pro-invasion signals in the form of active extracellular vesicles and metabolites that promote TME escape and migration.57 Such behaviors relate heavily to more realistic treatment regimens and necessitate a detailed transient model that will be scrutinized in a future study. With large-scale spheroid drug screening platforms such as the inverted colloidal crystal (iCC) platform by Ref. 14 becoming popularized, there is a growing need for standardized, physiologically relevant metrics to assess treatment efficacy within and between studies.

Despite the promise of powerful universal threshold values like those proposed in this work, it is key to remember the limitations of mathematical modeling in biological systems. In real hands-on cell-culture and drug-screening experiments, the systems are usually non-ideal, and errors caused by simplifications from presupposed assumptions may amplify. Two key examples when interpreting the work presented in this article are (1) irregularity of spheroid shape leading to breakdown of the assumption of spherical radial symmetry and (2) internal granularity caused by variable cell shape and distribution within the spheroid leading to the breakdown of the assumption of internal uniformity. We address the first by analyzing the spheroids individually rather than as a population, which allows us to apply circularity cutoffs (>0.85) and limit our analysis to systems that may be reasonably represented by our model. We address the second by setting a minimum spheroid size threshold of 50 μm in diameter, where any spheroid below would have a non-negligible granularity, with each spheroid being consistent of fewer than five cell layers. Even with all this considered, our work details physically informed relationships that go deeper than traditional $[eqn]$ evaluations and provides a framework for evaluating drug performance in complex cancer models.

METHODS

IV.

Simulation of fluorescence image projection libraries for radial hypoxia and viability profiles

A.

Spherical radial distribution profiles were generated using Python for cases of hypoxia, and viability was generated using Eqs. (2) and (9), respectively. Dimensionless radial coordinates ranged from 0 to 1 with step sizes of 0.01. For the hypoxia distributions, Thiele moduli ( $[eqn]$ ) ranging from $[eqn]$ to $[eqn]$ were simulated exponentially with step sizes of $[eqn]$ . For the viability distributions, dimensionless $[eqn]$ -tuples ranging from 0 to 1 were simulated with step sizes of 0.05. Radial profiles were simulated for a differential hemispheric slice of the spheroid parallel to the direction of incident illumination. Decay profiles were generated for simulated diameters of 0–500 μm with step sizes of 2 μm. The simplified decay function of $[eqn]$ derived from Eq. (7)—where z is the depth from the incident surface of the spheroid in micrometers and A is the simplified constant yielding $[eqn]$ and $[eqn]$ —was applied to the appropriate unattenuated profile.36 Each light-attenuated profile was summed for each r-coordinate along the axis of incident illumination and normalized to the maximum summed intensity value, producing the resultant radial projection profile.

Hypoxia assessment of 3D cell cultures of PANC-1, 3R-4, and 3R-5-1 spheroids

B.

Spheroids derived from the pancreatic cancer cell line PANC-1 and pancreatic cancer patient-ascites-derived cancer cells (3R-4 and 3R-5-1; IRB-TPEVGH No. 2023-04-005BC#2) were generated using the GEcoll ACD 3D Culture Kit according to the manufacturer's protocols.58 Briefly, cells were seeded at a density of 1000 cells per 50 μl of hydrogel per well in 24-well plates. The cell–hydrogel mixture was incubated at 37 °C in a 5% $[eqn]$ atmosphere for 14 days. The culture medium was replaced every 3 days to sustain cell growth and maintain appropriate nutrient levels.

After 14 days of spheroid culture, hypoxia was assessed using the BioTracker™ 520 Green Hypoxia Dye (Merck). The dye was prepared according to the manufacturer's protocols to achieve a stock concentration of 1 mM. The stock solution was further diluted with DMEM to a final staining concentration of 5 μM. Before staining, the culture medium was removed, and the spheroids were rinsed twice with phosphate-buffered saline (PBS). The prepared stain solution was added directly to the spheroids and incubated at 37 °C for 1 h. Following incubation, the dye was removed, and the spheroids were washed twice with PBS. Fresh DMEM was added, and spheroids were incubated for an additional 3 h prior to imaging.

Spheroid viability was evaluated using the Cyto3D® Live-Dead Assay Kit (TheWell Bioscience) following the manufacturer's protocol. The Cyto3D reagent was brought to room temperature and then directly added at a ratio of 20 μl per 1000 μl of total volume in each well. The spheroids were then incubated at 37 °C for 15 min to allow proper staining before analysis.

Imaging of stained spheroids was performed using an Olympus IX83 Motorized Inverted Fluorescence Microscope. Hypoxia staining was visualized at excitation/emission wavelengths of approximately 498/520 nm. Live-dead assays were visualized using fluorescence microscopy to distinguish viable (AO, 494/517 nm, green fluorescence) from non-viable (PI, 535/617 nm, red fluorescence) cells. All images were captured and processed using the microscope's associated imaging software, ensuring consistent magnifications and exposure settings for accurate comparative analyses.

Micrographs were processed using ImageJ to derive radial distribution profiles. The Analyze Particles function was used to determine spheroid centroids and approximate sizes. Radial Profile Plot and Specify ROI plugins were used to reproducibly extract radially normalized intensity profiles.59,60 Intensity profiles were normalized to maximum intensities and spheroid radius in Python. Least squares regression fitting was used to fit curves given the approximate spheroid size using the appropriate simulated hypoxia and viability curve libraries. Best fit values of $[eqn]$ and $[eqn]$ were selected from these fittings.

Drug testing of 3D cell cultures of LNCaP spheroids in degradable hydrogel droplets

C.

LNCaP cells were obtained from ATCC and cultured in RPMI-1640 medium modified to contain 2 mM L-glutamine, 10 mM HEPES, 1 mM sodium pyruvate, 4500 mg/l glucose, 1500 mg/l sodium bicarbonate (ATCC, Manassas, VA, USA), and 1 vol. % of antibiotic–antimycotic (Life Technologies Corporation, Carlsbad, CA, USA). Cell lines were incubated at 37 °C with 5% $[eqn]$ and routinely tested for mycoplasma contamination and were negative throughout the present study.

Norbornene-modified hyaluronic acid (NorHA) with a 22% degree of substitution (DS) was synthesized following established procedures.50,61 To quantify the DS, lyophilized polymer was dissolved in deuterium oxide ( $[eqn]$ ) and analyzed via ^1^H-NMR (Bruker AVANCE III HD 400 Nanobay).61 Polymer solutions were generated by dissolving NorHA (2% wt./vol.) in 1× Dulbecco's phosphate-buffered saline (DPBS) (Corning Life Sciences, MA, USA) with 2.1 mM DL-dithiothreitol (DTT) (Sigma-Aldrich, LO, USA) and 0.5 mM thiolated MMP-sensitive cross-linker (GCRDGPQG↓IWGQDRCG, MW: 1754.0 Da; down arrow indicates the site of proteolytic cleavage, GenScript). The water-soluble photoinitiator 2-Hydroxy-4′-(2-hydroxyethoxy)-2-methylpropiophenone (Irgacure 2959, Sigma-Aldrich) was used at 0.05% (wt./vol.). LNCaP cells were passaged and resuspended in the NorHA solutions at $[eqn]$ .

Cell-laden microgels were produced through emulsification via pipetting using a 008-FluoroSurfactant (1 wt. %) dissolved in HFE7500 (RAN Biotechnologies, Inc., MA, USA). A homemade pipette microfluidic device was used for microdroplet generation. Elliptical pipette tips were fabricated following our previous reports.62,63 Tips were deformed with a torque screwdriver using 20 in.-pounds. The microdroplets were exposed to UV light (Omnicure S2000, Excelitas Canada Inc., Canada) with a power of 20 mW cm^2^ for 2 min to prompt their gelation via thiol-ene reaction.64 Microgels were collected using a cell strainer and washed with 30% 1H,1H,2H,2H-Perfluoro-1-decanol (Sigma-Aldrich, LO, USA) in HFE7500 for breaking the emulsion. Further washes were made with PBS to purify the microgels.

Cell-laden microgels were incubated, and cell viability was tested on days 1 and 8. A range of doxorubicin (DOX) solutions (0.01–10 μM) in cell medium was used to culture spheroids grown for 8 days for an additional 24 h. Spheroids were stained using calcein and SYTOX (Thermo Fisher, Waltham, MA, USA) following the procedures recommended by the manufacturer. Fluorescence microscopy (ECHO Revolve, San Diego, CA, USA) was utilized to take micrographs of the spheroids. $[eqn]$ was calculated based on the fluorescence signal using an exponential decay fit.

Drug testing of 3D cell cultured HepG2 spheroids in inverted colloidal crystal (iCC) framework

D.

Alginic acid sodium salt (180947), calcium chloride dihydrate ( $[eqn]$ ; 223506), anhydrous sodium hydroxide (NaOH; S5881), and alginate lyase (A1603) were purchased from Sigma-Aldrich (MO, USA). Cryopres™ dimethyl sulfoxide (DMSO; 092780148) was obtained from MP Biomedicals (CA, USA). Certified™ low-melt agarose (161-3111) was sourced from Bio-Rad (CA, USA). Sylgard™ 184 silicone elastomer (PDMS; 102092-312) was purchased from VWR (PA, USA). Minimum Essential Medium Eagle (MEM; 10-010-CV) and phosphate-buffered saline (PBS; 21-040-CM) were obtained from Corning (NY, USA). Fetal bovine serum (FBS; 1300-500) of USDA-approved origin was sourced from Avantor Seradigm (PA, USA) and stored at $[eqn]$ . Antibiotic–antimycotic solution (15240096) was purchased from Fisher Scientific (MA, USA). Trypsin-EDTA (25200072) and BOBO-3 Iodide (570/602) were obtained from Thermo Fisher Scientific (MA, USA). 4% paraformaldehyde in 0.1 M phosphate buffer (15735-50S) was sourced from Electron Microscopy Sciences (PA, USA). Cell Explorer™ Live Cell Tracking Kit (22620) was purchased from AAT Bioquest (CA, USA). Paclitaxel (103607-762) was purchased from VWR (PA, USA).

The inverted colloidal crystal (iCC) framework was fabricated based on our previously reported method,14 with minor adjustments for enhanced scalability and throughput. Monodisperse alginate microgels were generated via electrospraying a 2.5% wt./vol. alginate solution into a 200 mM $[eqn]$ bath using the Spraybase® electrospraying system. A 15.6% wt./vol. low-melt agarose solution was prepared and thermally processed to form a moldable hydrogel, which was then reshaped into 1.5 mm thick disks using thermal cycling inside a glass-bottom dish.

To assemble the microgel template, $[eqn]$ alginate microgels were suspended in 200 mM $[eqn]$ , deposited into the glass-bottom dish, and covered with a cover glass. Mild sonication removed excess liquid, facilitating microgel array formation. The agarose hydrogel was located over the array, covered with a coverslip, and heated at 70 °C for 30 min to ensure full embedding, followed by cooling at room temperature to solidify infiltrated agarose hydrogel. Embedded microgels were removed by sequential enzymatic and ionic degradation: incubation with alginate lyase (1 mg/ml) at 37 °C for 24 h, followed by extensive washing with de-ionized water. Residual calcium ions were chelated with five rounds of PBS washes (200 ml each) at 1-h intervals, making the iCC framework transparent. The final iCC framework was rinsed and UV-sterilized before cell seeding.

Human hepatocellular carcinoma cells (HepG2; HB-8065, ATCC, VA, USA) were cultured in MEM supplemented with 10% FBS and 1% antibiotic–antimycotic. Cells were maintained at 37 °C with 5% $[eqn]$ in a humidified incubator (MCO-15AC; Sanyo, Japan) and used between passages 6 and 12. Upon reaching >80% confluency, cells were trypsinized and resuspended at $[eqn]$ . For live tracking, cells were incubated with 20 μM Cell Explorer™ dye in 2 ml culture medium for 30 min, followed by PBS wash and resuspension in fresh medium at the same concentration.

To seed the iCC framework, 4.5 ml of fresh culture medium was added to the framework positioned inside a PDMS-molded, perpendicularly oriented 50 ml centrifuge tube. A 250 μl aliquot containing $[eqn]$ cells was gently mixed into the setup and centrifuged at 400 rpm for 3 min to initiate cell loading. Excess cells outside the framework were gently redistributed around the structure. This process was repeated to achieve a total loading of $[eqn]$ cells per framework. The seeded framework was transferred to a non-treated 60 mm Petri dish containing 5 ml of fresh medium and incubated at 37 °C with 5% $[eqn]$ for 6 days. The medium was replaced every 48 h.

Paclitaxel was dissolved in DMSO (50 mg/ml; 58.6 mM) and diluted into complete culture medium to final concentrations of 10, 100, and 1000 μM (2 ml each). Spheroids cultured within the iCC framework were treated with the respective concentrations for 24 h. Post-treatment, spheroids were stained with BOBO-3 Iodide for 15 min to visualize compromised/dead cells. Multiphoton confocal imaging was performed using the Leica Stellaris 8 DIVE system (Wetzlar, Germany), capturing dual-channel images for live (blue) and dead (red) cell fluorescence. The imaging parameters were as follows: 24 z-slices over 240 μm with a voxel dimension of [1.354, 1.354, 10] μm.

Radial intensity profiles for total and dead cells were generated using central cross sections of the Blue CellTracker™ and red BOBO-3 stained micrographs using the Radial Profile Plot and Specify ROI plugins. The radial viability distribution ( $[eqn]$ ) was estimated as the ratio of red dye intensity against blue dye intensity. Approximating each spheroid to a rough sphere, the net viability for the spheroid in each treatment case was derived using Python as

$$[eqn]$$

Statistical analysis

E.

The data presented here are expressed as the mean and standard deviation or standard error of the mean as calculated using Excel, Python, or GraphPad Prism 8 software functions. Errors are propagated through formulas using standard propagation of uncertainties techniques.

SUPPLEMENTARY MATERIAL

See the supplementary material for supplementary Figs. 1–5 and supplementary Table 1.