2.1.

Patient-Derived Cancer Organoid Culture

PCOs were generated from two patient-derived metastatic colorectal cancer lines (P1 and P2), according to a previously published protocol.⁴⁷ These lines were both derived from colorectal cancer metastases to the liver. Five wells (one well for each of the five treatment groups) were generated in each of the 24-well plates. The base culture medium used was DMEM/F12 (ThermoFisher) supplemented with 10% FBS (Sigma), and 1% penicillin–streptomycin (Sigma). Briefly, previously grown PCOs were singularized with 0.25% trypsin, resuspended in base medium, and mixed with Matrigel (corning # 354,234) at a 1:1 ratio. The cell–Matrigel mixture was pipetted onto the glass surface of each well in a 24-well glass bottom plate (black frame, #0 cover glass, Cellvis, P24-0-N). The plate was incubated at 37°C for 2 to 3 min to allow the mixture to solidify and then inverted to ensure PCOs were suspended in the Matrigel. The mixture was left to solidify for at least 20 min, then base medium supplemented with Wnt3a-conditioned medium (1:1 ratio) in each well [Fig. 1(a), organoid culture]. Medium was refreshed every two days. Aberrant Wnt/$\beta$-catenin signaling is a hallmark of colorectal cancer, and previous studies show that medium conditioned with Wnt3a supports colorectal PCO growth.⁴⁸ Wnt3a-conditioned medium is generated by harvesting medium, in which murine L Wnt3a cells (ATCC, CRL-2647) have been cultured.

Fig. 1

Redox imaging of PCOs. An overview of the protocol for redox imaging and quantitative image analysis. (a) Graphical protocol showing the culturing of PCOs, pretreatment redox imaging, drug treatment, and posttreatment redox imaging time course. Redox imaging uses pairs of NAD(P)H and FAD fluorescence images from the same field-of-view to calculate the ORR image [NAD(P)H/(FAD + NAD(P)H)]. “+” in first well indicates registration mark. (b) Graphical representation of the image analysis pipeline, which includes registration of all frames in each image time series, organoid segmentation, single-organoid tracking, and quantification of metabolic and morphological variables to capture PCO drug response at each time point. Scale bar: $500\mu \mathrm{m}$.

2.2.

Treatment Protocol

Five wells for each PCO line were grown for 14 days before treatment. At day 14, media in the wells were refreshed with normal growth medium or medium containing sodium cyanide, 2-deoxy-glucose (2DG), cisplatin, or paclitaxel [Fig. 1(a), drug treatment]. Specifics of each treatment group, including the mechanism of action, concentration used, and number of PCOs per group, are listed in Table 1.

Table 1

Treatment groups, drug concentrations, and number of PCOs (n) assessed for P1 and P2 lines.

2.3.

Redox Imaging

Redox imaging was performed using an inverted epifluorescence microscope (Nikon Ti-U). System specifications included a $4\times$ air objective (Nikon CFI Plan Fluor, NA 0.13, FOV: $3.33\mathrm{mm}\times 3.33\mathrm{mm}$, lateral resolution: $2.16\mu \mathrm{m}$ at $460\mathrm{nm}/2.46\mu \mathrm{m}$ at 525 nm), a white light LED source (SOLA FISH, Lumencor), a scientific grade CMOS camera (Flash4, Hamamatsu, image pixel size: $1.625\mu \mathrm{m}$), and an $XYZ$ automated stage (MS-2000, ASI). Microscope control was performed with $\mu \text{Manager}$.⁵³ NAD(P)H fluorescence was excited with a DAPI filter cube (Nikon, ex: 361 to $389\mathrm{nm}/\mathrm{em}$: 435 to 485 nm) and integrated over 3 s. FAD fluorescence was excited with a GFP filter cube (Nikon, ex: 450 to $490\mathrm{nm}/\mathrm{em}$: 500 to 550 nm) and integrated over 5 s.

An NAD(P)H fluorescence image and FAD fluorescence image were acquired for each well pretreatment and at 20 min, 12, 24, 36, and 48 h posttreatment [Fig. 1(a), redox imaging]. The NAD(P)H and FAD images were used to calculate the 2D ORR for each pixel position $(i,j)$ using the following equation:

The entire imaging protocol including $XYZ$ stage repositioning and image acquisition was automated with $\mu \text{Manager}$. For each plate, a total of five wells were imaged. Total imaging time per plate was $<2\min$, capturing 48 PCOs across five wells for P1 and 35 PCOs across five wells for P2. Note that the first well of the 24-well plate was used to assist in positioning the plate before imaging, which enabled reproducible measurements of the same field-of-view over multiple time points [Fig. 2(a)].

Fig. 2

Single-organoid tracking. Organoid tracking over multiple time points was achieved with a registration mark and single-particle tracking. (a) A registration mark was inscribed on the bottom of the first well of the 24-well plate and referenced at each time point to ensure accurate $XYZ$ positioning. (b) Organoid tracking was achieved by registering the NAD(P)H image time series, segmenting the organoids in each frame, and detecting the centroids of each segmented organoid. These centroids were tracked using an open-source single-particle tracking tool (TrackMate). (c) A segmented organoid mask at 48-h posttreatment is overlaid with the centroid tracking results at each imaging time point. Scale bar: $500\mu \mathrm{m}$.

2.4.

Quantitative Image Analysis

A framework for redox imaging and automated image analysis was developed to track each PCO across multiple time points and quantify organoid metabolism and morphology [Fig. 1(b)]. Linkage of each organoid over time was performed via a combined approach of NAD(P)H image time series registration and a 2D centroid matching-based approach with single-particle tracking. This pipeline was implemented using MATLAB (2019b, Mathworks), except when specifically noted.

Registration of each image time series was achieved using both a registration mark and rigid image registration [Figs. 2(a) and 2(b)]. A registration mark was etched on the glass bottom of the first well of each plate using a diamond-tip scribe and was referenced at every time point to precisely position the plate before imaging to minimize $XYZ$ drift [Fig. 2(a)]. After acquisition, $XY$ drift between frames was corrected through rigid image registration, which registers two frames by maximizing the cross correlation of pixel values.^54,55 The NAD(P)H image time series was used to calculate the $XY$ shift values needed to align the frames in the NAD(P)H image time series, and these shift values were applied to the corresponding FAD and ORR image time series. These shift values were also used to determine areas that are common across all frames. Image areas that did not overlap across the image time series were excluded.

Segmentation of the PCOs was independently performed on the NAD(P)H image at each time point to capture changes in organoid size and shape [Fig. 2(b)]. However, segmentation based on intensity alone is challenging because autofluorescence images have a low signal-to-background ratio (typically between 2 and 3). A segmentation algorithm based on edge enhancement was developed using MATLAB that reliably segmented all PCOs imaged in this study. Algorithm 1 is as follows, below.

Algorithm 1

Organoid Segmentation

Specific values for algorithm parameters were chosen empirically based on visual inspection. Validation of the segmentation algorithm was performed by calculating the Sørenson–Dice similarity coefficient for pairs of automatically segmented images and manual segmented images (Fig. S1 in the Supplemental Materials, $n=5$, $\text{mean}\pm \mathrm{SEM}=0.93\pm 0.05$).

Tracking of each PCO was achieved using an open-source single-particle tracking tool (TrackMate, ImageJ/FIJI).⁵⁷ TrackMate achieves single-particle tracking by detecting particles in each frame, linking detected particles from frame-to-frame, and combining all links into the most likely tracks.⁵⁸ Image time series were generated using MATLAB with 2D Gaussian particles (kernel: $10\times 10\text{pixels}$) at the centroid of each segmented PCO [Fig. 2(b), centroid detection]. These centroids were then identified in each frame using the difference-of-Gaussian blob detector (options: 10 pixel diameter, no spot filtering). The linear assignment problem algorithm was chosen to perform particle tracking because it is sensitive to events such as particle drift, particle merging, or over-/under-segmentation (options: 200 pixel maximum linking distance, 200 pixel max gap-closing distance, and 2 frame max gap-closing distance). All key parameters were initially chosen by visual inspection of the tracking results for the P1 control time series and then applied to all other image time series [Fig. 2(c)]. Once completed for an image time series, the single-organoid tracking results were used to link the metabolic and morphological variables for each organoid over time.

Quantification of variables that describe organoid morphology and metabolism was performed for each PCO at every time point ($n=83$ PCOs total, $n=48$ PCOs for P1, and $n=35$ PCOs for P2). In total, 24 variables (12 metabolic and 12 morphological) were quantified for each segmented organoid: mean, minimum, maximum, and standard deviation of the ORR, NAD(P)H intensity, and FAD intensity values; organoid area, perimeter, solidity, extent, eccentricity, circularity, minimum Feret’s diameter, maximum Feret’s diameter, minor axis, major axis, convex area, and equivalent diameter (regionprops). To limit the effects of variation in the culture medium autofluorescence, no background subtraction was performed on the images before quantification. Detailed definitions of these variables are provided in Table S1 in the Supplemental Materials. PCOs connected to the image border were excluded from statistical analysis.

2.5.

Statistics

All statistical analyses were performed using MATLAB or R. Open-source MATLAB toolboxes used were Gramm, a data visualization and statistics package, and drtoolbox, a dimensionality reduction package.^59,60 Open-source R toolboxes included nlme, a mixed effect modeling package, and lsmeans, a least-squares means estimation package.^61,62

Linear regression was also used to assess the bivariate relationship between organoid area and ORR, NAD(P)H intensity, and FAD intensity pretreatment (stat_glm and Gramm).⁶⁰ The linear pairwise correlations were also calculated for all 24 variables pretreatment (corr).

Organoid-level-normalized time series were calculated for each organoid as $X_{\text{posttreatment}}/X_{\text{pretreatment}\text{organoid}}$, where $X$ is the mean ORR or organoid area. Linear mixed-effect models were used to analyze these organoid-level-normalized time series for each condition and PCO line (nlme, R).⁶¹ Linear mixed-effect models can account for organoid-level differences in measurements using a random-intercept model, which provides a way to assess treatment-induced effects in the presence of organoid-level heterogeneity.^63,64 A random-intercept model was specified as $Y\sim \text{time}+\text{treatment}+\text{treatment}*\text{time}+(1|\text{organoid})+\text{error}$, where $Y$ is either the pretreatment-normalized ORR or organoid area, time is a categorical variable encoding time, treatment is a categorical variable encoding the five treatment groups, organoid is the organoid-level-normalized value at 20-min posttreatment, and error is the random error in the model. Time is defined as a categorical variable to avoid the assumption of linearity between time and the ORR or organoid area. Organoid here is encoded as a random effect (“1|organoid” in the statistical model notation) to account for variability at the organoid level. Pretreatment values were excluded from the analysis because all time points were normalized to pretreatment values (i.e., pretreatment variance is one). Thus the first time point included in the analysis is 20-min posttreatment. A first-order autoregressive [i.e., AR(1)] covariance structure was used to account for the correlations found in these time-course data, where the correlation between time points decreases with increasing separation in time. Least-squares means was then used to calculate the pairwise differences between all treatment groups at each time point, compute the $p$-value for each pairwise difference using Tukey’s “honest significant difference” method, and account for multiple comparisons (lsmeans, R).⁶²

Comparison between single-organoid tracking and pooled analysis was performed. The pooled analysis was the same as the single-organoid tracking analysis, except for the normalization approach and statistical model. To mimic pooled organoid data without organoid tracking, data for each organoid was normalized to the well-level pretreatment mean as $X_{\text{posttreatment}}/X_{\text{pretreatment well}}$, where $X$ is the mean ORR or organoid area. A fixed effects model was specified as $Y\sim \text{time}+\text{treatment}+\text{treatment}*\text{time}+\text{error}$, where $Y$ is either the well-level-normalized ORR or organoid area, time is a categorical variable encoding time, treatment is a categorical variable encoding the five treatment groups, and error is the random error in the model. Pretreatment values were also not included in the analysis. All other parts of the pooled analysis were the same as the organoid-level analysis including the AR(1) covariance structure, least-squares mean for pairwise comparisons, and correction for multiple comparisons.

Principal component analysis (PCA) was used to explore the multivariate relationships between the 24 variables measured from each PCO. PCA was performed using an open-source MATLAB toolbox for dimensionality reduction (compute_mapping, out_of_sample, and drtoolbox). Data from all PCOs pretreatment, 24-h posttreatment, and 48-h posttreatment were standardized within each time point before PCA to account for differences in the scale of each variable ($z\text{score}$). The loadings that projected the data into the principal component (PC) basis were computed using the data from all PCOs pretreatment. These loadings were then applied to data at 24-h and 48-h posttreatment to visualize the time-dependent effects of treatment on the variables from each PCO. Loading vectors were plotted to visualize the contribution of each variable to the PCs. For clarity, only the top 12 variables that loaded on PC1 and PC2 were displayed as they are the variables that contributed the most variability. The projected data were scaled to fit in the loadings interval (biplot). The scaled scores of the data projected on the first two PCs were plotted and fitted with ellipses (i.e., bivariate normal distribution) that contained 95% of the points for each group (geom_point, stat_ellipse, and Gramm).

3.1.

Optical Redox Ratio is Independent of Organoid Morphology

The goal of this study was to segment and track individual PCOs with wide-field one-photon redox imaging to extract morphological and metabolic variables of drug response. A quantitative image analysis framework was developed to automatically quantify 12 metabolic variables that summarize NAD(P)H intensity, FAD intensity, and ORR, along with 12 morphological variables that quantify organoid size and shape (Table S1 in the Supplemental Materials), for each organoid over time. Relationships between the area of a PCO and its mean NAD(P)H intensity, FAD intensity, and ORR are shown in Figs. 3(a)–3(c). Linear models were fit to paired measurements from all PCO pretreatment ($n=48$ PCOs for P1, $n=35$ PCOs for P2). There is no correlation between the mean ORR of each PCO and its area [Fig. 3(a), $r=0.04$, $p=0.70$]. However, there is a moderate positive correlation between organoid area and mean NAD(P)H intensity [Fig. 3(b), $r=0.48$, $p<0.001$] and between organoid area and mean FAD intensity [Fig. 3(c), $r=0.36$, $p<0.001$].

Fig. 3

Analysis of organoid-level redox imaging variables. Analysis of variables extracted from all PCOs pretreatment. (a) Linear regression was used to quantify the relationship between organoid area and its mean ORR. The fitted model indicates the mean ORR value of each PCO is independent of its area ($r=0.04$, $p=0.70$) for all PCOs from patient 1 (P1) and patient (P2) pretreatment. Statistical significance of the model coefficients is indicated as *** for $p<0.001$. (b) Linear regression was used to quantify the relationship between organoid area and its mean NAD(P)H intensity. The fitted model indicates a moderate positive correlation ($r=0.48$, $p<0.001$) between organoid area and mean NAD(P)H intensity. (c) Linear regression was used to quantify the relationship between organoid area and its mean FAD intensity. The fitted model indicates a moderate positive correlation ($r=0.36$, $p<0.001$) between organoid area and mean FAD intensity. (d) Linear pairwise correlation (Pearson’s $r$) matrix of the variables extracted from each PCOs. Heatmap shows only statistically significant correlations ($p<0.05$); n.s., not significant ($p>0.05$, black boxes). Table 1 shows the number of PCOs in each treatment group for each patient.

The linear pairwise correlations of all 24 variables were computed to understand the relationships between the metabolic and morphological variables [Fig. 3(d)]. As expected, morphological variables were significantly correlated ($p<0.05$) with most other morphological variables apart from eccentricity, which was only significantly correlated with the minor axis length, solidity, and circularity. Likewise, metabolic variables were significantly correlated ($p<0.05$) with many other metabolic variables except for the minimum and maximum values, which were uncorrelated with many of the other metabolic variables. Statistically significant correlations between morphological and metabolic variables were also observed, especially for standard deviations in FAD intensity and ORR that positively correlate with greater organoid area/diameter/length. Consistent with the linear model [Fig. 3(a)], the mean ORR had no significant correlations with any morphological parameters. Previous studies have shown that the ORR of cancer organoids predicts treatment response,^9–12 so these data indicate that the mean ORR may provide information about PCO drug response that is independent of morphology.

3.2.

Drug Response Lowers Optical Redox Ratio

Changes in the ORR were calculated in two PCO lines over a 48-h treatment time course. Previous studies have found that a treatment-induced decrease in the ORR of cancer organoids predicts in vivo treatment response.^25,40 Treatment-induced changes in mean NAD(P)H intensity and mean FAD intensity can be found in Fig. S2(A)–(D) in the Supplemental Materials, which reflects the trends in the ORR time series. Both PCO lines were treated with cisplatin and paclitaxel, two chemotherapies, and sodium cyanide and 2DG, two metabolic inhibitors. There was one well for each treatment and an additional untreated well served as a negative control (five wells total). Three representative ORR image time series are shown in Fig. 4(a), which show the change in ORR with cyanide treatment, large change in ORR with paclitaxel treatment, and the stability of the ORR in untreated organoids. Patient 1 PCOs treated with cyanide and paclitaxel were representative of a minimal and large response at 48-h posttreatment, respectively. Arrows in Fig. 4(a) indicate specific PCOs that are representative of changes observed for that condition. Single-organoid tracking was used to normalize each PCO to its pretreatment ORR, so that relative changes in ORR could be compared across PCOs. Least-squares means was used to estimate the mean of each group and to calculate the pairwise differences at each time point between treatment groups [Figs. 4(b) and 4(c)].

Fig. 4

Treatment-induced changes in PCO ORR. Quantitative image analysis tracked changes in the ORR of each PCO over time. (a) Representative ORR image time series for patient 1 (P1) and patient 2 (P2). Arrows indicate specific organoids representative of changes observed for that condition. B.G, background. (b) Heatmaps show the pairwise percent differences between pretreatment-normalized ORR between all treatment groups for P1. Top: pairwise percent differences between each treatment and control for P1. Bottom: pairwise percent differences between drug treatments (excluding control) for P1. (c) Heatmaps show the pairwise percent differences between pretreatment-normalized ORR between all treatment groups for P2. Top: pairwise percent differences between each treatment and control for P2. Bottom: pairwise percent differences between drug treatments for P2. All pairwise differences were calculated from the linear mixed-effect models via least-squares means. n.s., not significant ($p>0.05$, black boxes). Table 1 shows the number of PCOs in each treatment group for each patient. Scale bar: $500\mu \mathrm{m}$.

For P1 PCOs, the ORR decreased compared to control at 36 h for the cyanide ($p<0.01$) and 2DG ($p<0.01$) groups [Fig. 4(b)]. Cisplatin treatment in P1 resulted in a decreased ORR compared to control only at 48 h [Fig. 4(b), $p<0.05$]. Paclitaxel treatment in P1 resulted in a decreased ORR compared to control starting at 12 h [Fig. 4(b), $p<0.05$]. For P2 PCOs, the ORR decreased compared to control at 12 h for the 2DG ($p<0.01$), cisplatin ($p<0.05$), and paclitaxel ($p<0.05$) groups [Fig. 4(c)]. Cyanide treatment in P2 resulted in a decreased ORR compared to control starting at 36 h [Fig. 4(c), $p<0.01$]. Additionally, cyanide in both P1 and P2 caused a sharp rise in the ORR ($p<0.001$) at 20 min, characteristic of electron transport chain inhibition [Figs. 4(b) and 4(c)]. Heatmaps [Figs. 4(b) and 4(c)] also show the pairwise differences between treatment groups.

3.3.

Drug Response Slows Organoid Growth

Changes in organoid size or growth rate can occur due to effective drug response, so the relative change in organoid area over time was also quantified. All other morphological variables remained constant or had a limited change over 48 h compared to pretreatment values [a subset can be found in Fig. S2(E)–(J) in the Supplemental Materials]. Three representative image time series in Fig. 5(a) show segmented PCOs with the color-coded value mapped to the pretreatment normalized organoid area. These image time series illustrate increased organoid area over time in control conditions and relatively modest increases in organoid area over time in cyanide and paclitaxel treatment conditions. Arrows in Fig. 5(a) indicate specific PCOs representative of changes observed for that condition. Heterogeneity in organoid area was observed pretreatment, so single-organoid tracking was used to normalize each organoid to its pretreatment area. This enabled comparisons of relative changes in area across organoids. Least-squares means was used to estimate the mean of each group and to calculate the pairwise differences at each time point between control and treatment groups [Figs. 5(b) and 5(c)].

Fig. 5

Treatment-induced changes in PCO areas. Quantitative image analysis tracked changes in the area of each PCO over time. (a) Representative segmented image time series for patient 1 (P1) and patient 2 (P2). The color-coded value of each segmented organoid corresponds to value at each time point divided by its pretreatment value. Arrows indicate specific organoids that are representative of changes observed for that condition. (b) Heatmaps show the pairwise percent differences between pretreatment-normalized area between all treatment groups for P1. Top: pairwise percent differences between each treatment and control for P1. Bottom: pairwise percent differences between drug treatments (excluding control) for P1. (C) Heatmaps show the pairwise percent differences between pretreatment-normalized area between all treatment groups for P2. Top: pairwise percent differences between each treatment and control for P2. Bottom: pairwise percent differences between drug treatments for P2. All pairwise differences were calculated from the linear mixed-effect models via least-squares means. n.s., not significant ($p>0.05$, black boxes). Table 1 shows the number of PCO in each treatment group for each patient. Scale bar: $500\mu \mathrm{m}$.

Two separate patterns of response emerged for metabolic inhibitors and chemotherapies. For P1 groups treated with metabolic inhibitors, lower organoid areas compared to control were observed starting at 12 h for the cyanide ($p<0.01$) and 2DG ($p<0.01$) [Fig. 5(b)]. However, paclitaxel-treated P1 PCOs continued to grow until 24 h and were only significantly different from control beginning at 36 h ($p<0.05$) [Fig. 5(b)]. For P1 PCOs, cisplatin treatment resulted in lower organoid areas compared to control only at 48 h ($p<0.05$) [Fig. 5(b)]. For P2 PCOs treated with metabolic inhibitors, lower organoid areas compared to control were observed starting at 12 h for the cyanide ($p<0.01$) and 2DG ($p<0.01$) groups [Fig. 5(c)]. For P2 PCOs treated with chemotherapies, continued growth was observed until 24 h, and organoid area was only lower than control beginning at 36 h for paclitaxel ($p<0.01$) and cisplatin ($p<0.01$) groups [Fig. 5(c)]. Heatmaps [Figs. 5(b) and 5(c)] also show the pairwise differences between treatment groups.

3.4.

Single-Organoid Tracking Provides Improved Sensitivity

PCOs are heterogeneous and exhibited a range of ORR and organoid area values. Single-organoid tracking can monitor dynamic changes in organoid-level variables over time by normalizing each PCO to its own pretreatment value. This enabled the observation of the ORR and organoid area on an organoid-level in the P1 and P2 control groups. In both P1 and P2 control groups, the ORR was stable over 48 h [Figs. 6(a) and 6(c), $p>0.05$ pretreatment versus each time-point posttreatment using organoid-level normalization]. In both P1 and P2 control groups, PCOs continue to grow over 48 h, with a statistically significant increase in organoid area of 10% and 11%, respectively, at 48 h compared to pretreatment [Figs. 6(e) and 6(g), $p<0.05$ pretreatment versus each time-point posttreatment using organoid-level normalization].

Fig. 6

Comparison of single-organoid tracking and pooled analysis. Comparison of single-organoid tracking and pooled analysis using time series data from patient 1 (P1) and patient 2 (P2). Normalization refers to dividing the value of each organoid at each time point by either its own pretreatment value (organoid-level normalization) or the well-level mean (well-level normalization). (a), (b) Pretreatment normalized ORR using organoid-level normalization (left) and well-level normalization (right) for each treatment group for P1. (c), (d) Pretreatment normalized ORR using organoid-level normalization (left) and well-level normalization (right) for each treatment group for P2. (e), (f) Pretreatment normalized organoid area using organoid-level normalization (left) and well-level normalization (right) for each treatment group for P1. (g), (h) Pretreatment normalized organoid area using organoid-level normalization (left) and well-level normalization (right) for each treatment group for P2. Data in (a), (c), (e), (g) were analyzed using a linear mixed-effect model that uses single-organoid tracking to account for organoid-level variability. Data in (b), (d), (f), (h) were analyzed using a linear mixed-effect model that does not account for organoid-level variability due to the lack of single-organoid tracking. Data plotted are mean (lines) ± standard error of the mean (shaded area). Significant differences between a treatment and control ($p<0.05$) are indicated with a circle color-coded for the treatment group (see legend). Table 1 shows the number of PCOs in each group and for each patient.

To demonstrate the utility of single-organoid tracking, a pooled analysis was performed. In addition to the previous analysis using organoid-level normalization, data were analyzed using well-level normalization where each PCO was normalized to the well-level mean. Figure 6 shows the ORR and organoid area over time for P1 and P2 analyzed using single-organoid tracking (organoid-level normalization) or pooled analysis (well-level normalization). Single-organoid tracking [Figs. 6(a), 6(c), 6(e), and 6(g)] provided data with lower variability compared to pooled analysis [Figs. 6(b), 6(d), 6(f), and 6(h)], as the well-level mean does not fully account for organoid-level heterogeneity. In addition, pooled analysis does not provide organoid-level time series for use in a statistical model. Single-organoid tracking found more differences between control and treatment groups when compared to pooled analysis, which demonstrates the improved sensitivity to treatment response afforded by single-organoid tracking.

3.5.

Multivariate Analysis Reveals Patient-Derived Cancer Organoid Subpopulations

Two subpopulations of organoids with distinct morphology (solid and hollow) are present in cultures from P1, whereas all P2 organoids are morphologically similar [Fig. 7(a)]. The first P1 subpopulation, P1 solid, were generally dense (highly scattering), irregular in morphology, and had core regions with lower ORR values than the periphery. The second P1 subpopulation, P1 hollow, were less dense (lower scattering), had a spherical morphology, and had core regions with higher ORR values than the outer regions. In addition, the ORR values of P1 hollow PCOs were generally higher than P1 solid PCOs, driven by higher NAD(P)H intensity. P2 appeared to have one type of organoid morphology that was qualitatively similar to P1 solid: dense, irregular in morphology, and with a core region with a low ORR compared to the periphery. The differences in the internal structure of P1 solid and P1 hollow PCOs were confirmed by confocal microscopy of P1 solid and P1 hollow PCOs stained for epithelial cell adhesion molecule (EpCAM, Fig. S3 in the Supplemental Materials). By qualitative examination, 39 P1 solid, 9 P1 hollow, and 35 P2 PCOs were identified.

Fig. 7

Phenotypic identification of PCO subpopulations. Qualitative and quantitative (PCA) evaluation of PCOs with distinct phenotypes from patient 1 (P1) and patient 2 (P2). (a) Paired bright-field (top) and ORR (bottom) images of representative PCO from each of the three phenotypic subpopulations: P1 solid, P1 hollow, and P2. Scale bar: $250\mu \mathrm{m}$. B.G., background. (b) Loading vectors for the 12 variables with the highest loading on PC1 and PC2, which were defined with pretreatment data only. Green are metabolic variables and blue are morphological variables. (c) Data from all organoids projected onto PC1 and PC2. Colors correspond to the three qualitatively identified phenotypic organoid subpopulations: P1 solid, P1 hollow, and P2. Ellipse containing 95% of the points for each group were also plotted. (d) The pretreatment PCs applied to data from the PCOs at 24-h posttreatment. (e) The pretreatment PCs applied to data from PCOs at 48-h posttreatment.

PCA was applied to the 24 morphological and metabolic variables (Table S1 in the Supplemental Materials) for all pretreatment PCOs to determine if these variables could cluster PCOs by P1 solid, P1 hollow, and P2. A plot of the loading vectors [Fig. 7(b)] shows the contribution of each of the pretreatment variables on the first two PCs. The length of the vector corresponds to the magnitude of the contribution of a specific variable to the first two PCs. The angle between vectors as well as the angle between each vector and the PC axes indicates the level of correlation. The colors indicate which are metabolic or morphological variables. Only the 12 variables that have the largest contributions to the first two PCs were displayed for clarity. Figure 7(c) shows the organoid-level data projected onto the PC1–PC2 basis, with clear clustering of subpopulations based on P1 solid, P1 hollow, and P2. Notably, P1 hollow was correlated with metabolic variables, clustering in the lower left quadrant apart from the P1 and P2 clusters [Figs. 7(b) and 7(c)]. The same pretreatment PC1–PC2 basis was used to plot PCOs at 24- and 48-h posttreatment. At 24-h posttreatment, smaller variability is observed within and between clusters compared to pretreatment variability [Fig. 7(d)], and cluster variability increased at 48 h [Fig. 7(e)].