1. Introduction

Patient-derived cancer organoids (PCOs) are three-dimensional (3D) in vitro cultures isolated from tumor biopsies that recapitulate features of the in vivo tumor, including gene expression, cell morphology, and cell phenotype [1–3]. PCOs reflect in vivo drug sensitivity, offering a platform to perform high-throughput drug screens. Drug screens using PCOs have the potential to optimize treatment selection, improve treatment efficacy, and reduce side effects and toxicity [4,5]. Subclones captured in PCOs have also been shown to mirror the heterogeneity found in the original tumors, which often leads to treatment resistance [1,6,7]. This heterogeneity must be considered for robust predictions of patient drug sensitivity [7,8]. Heterogeneity manifests itself in PCOs both within (cell-level) and between (organoid-level) organoids, which motivates the assessment of large populations of PCOs to determine drug sensitivity. However, current methods to assess in vitro drug response in PCOs are limited because they are generally destructive, prohibiting time-course measurements (e.g., immunohistochemistry, single-cell sequencing/imaging, flow cytometry) or lack the ability to volumetrically image individual PCO size and morphology (e.g., brightfield microscopy, plate reader assays). There is a need to quantify volumetric growth and drug response in PCOs using non-destructive tools that are sensitive to organoid-level heterogeneity.

Optical coherence tomography (OCT) is a label-free imaging technique that detects light backscattered from a sample to reconstruct depth-resolved structure. OCT is based on low-coherence interferometry and provides rapid acquisition times of large 3D volumes (hundreds of mm³). Therefore, OCT can capture micron-scale volumetric information from large populations of PCOs. Additionally, OCT can track the volumetric growth of individual PCOs over time without destructive processing or exogenous labels. Previous studies have used OCT to monitor drug response in PCOs [9–15]. However, these previous studies have considered only a single PCO or pooled response across all PCOs in a well over time, and have not used single-organoid tracking to monitor heterogeneity in response on the individual PCO level.

Heterogeneity at the organoid-level results in sub-populations of organoids with differing treatment sensitivities [8,16]. Therefore, accurate OCT measurements of drug response requires imaging many PCOs over time to capture the response of these sub-populations. To compensate for organoid-level variability, the relative volumetric change from baseline for each PCO must also be captured. Accordingly, quantitative image analysis methods to automatically quantify and track organoids over time are needed for robust OCT-based high-throughput drug screens.

Here, we provide the first demonstration of 3D single-organoid tracking of PCOs using OCT. An image analysis framework was developed to automatically quantify organoid-level changes in volumetric growth. We measured the response to two chemotherapies (cisplatin, paclitaxel) and two metabolic inhibitors (sodium cyanide, 2-deoxy-glucose) across two patient-derived colorectal PCO lines. A quantitative image analysis pipeline was developed to track volumetric growth and drug response in each PCO before treatment and over multiple time points post-treatment up to 48 hours. Linear mixed-effect models were then applied to these organoid-level volumes over time, which provided improved sensitivity to drug response compared to pooling organoid volume over time (i.e., averaging all PCO volumes in a well, treatment group, or time point). Direct volume measurement using OCT provided more accurate quantification of PCO volume compared to a volume estimation technique based on 2D projections, similar to volume estimates commonly used with brightfield microscopy. This 2D volume estimation error was minimized in PCOs that were highly spherical. Lastly, volumetric growth rate was found to be heterogeneous within one patient, and this heterogeneity was attributed to PCO morphological sub-type (hollow or solid). This work demonstrates that OCT and 3D single-organoid tracking are attractive tools to monitor volumetric growth and drug response in PCOs, and that 3D single-organoid tracking provides improved sensitivity to PCO drug response and heterogeneity.

2. Methods

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 [17]. 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 # 354234) 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-3 minutes 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 minutes, then base medium supplemented with Wnt3a-conditioned medium (1:1 ratio) was added to each well. Medium was refreshed every two days. Aberrant Wnt/β-catenin signaling is a hallmark of colorectal cancer, and previous studies show that medium conditioned with Wnt3a supports colorectal PCO growth [18]. Wnt3a-conditioned medium is generated by harvesting medium in which murine L Wnt3a cells (ATCC, CRL-2647) have been cultured.

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 replaced with media containing sodium cyanide, 2-deoxy-glucose (2DG), cisplatin, or paclitaxel (Fig. 1(B), Drug Treatment). Media in the control wells were refreshed with unaltered media (no treatments). Specifics of each treatment group, including the mechanism of action, concentration used, and number of PCOs per group, are listed in Table 1.

 

Fig. 1. Optical Coherence Tomography and 3D Single-Organoid Tracking of Patient-Derived Cancer Organoids. An overview of the custom OCT instrument, experimental protocol, and quantitative image analysis pipeline. (A) A diagram of the custom swept-source OCT instrument used in this study. The 100kHz swept-source laser used for this study had a center wavelength of 1040nm and a bandwidth of 109nm. The off-axis inverted geometry enabled imaging of patient-derived cancer organoids in a standard glass-bottom multiwell plate. BD: Balanced detector, GP: XY galvanometer pair, FBG: fiber Bragg grating, Circ.: Fiber circulator, DC: Dispersion compensation, M: mirror. AI/AO: analog in/analog out. (B) Experimental protocol illustration including drug treatment and imaging time points. (C) Graphical representation of the image analysis pipeline, which includes registration of all frames in each image time series, 3D organoid segmentation, single-organoid tracking, and quantification of organoid volume. Scale bar: 500µm.

Download Full Size | PDF

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

View Table

2.3 Optical coherence tomography imaging

All OCT measurements were performed on a custom inverted swept-source OCT system (Fig. 1(A)) [23]. The swept-wavelength source (Axsun Technologies) had a center wavelength of 1040nm, bandwidth of 109nm, and a line rate of 100kHz. The calculated axial resolution was 6.5µm in culture medium (n = 1.33 at 1040nm). The lateral resolution was 16µm (Thorlabs, LSM03-BB, NA 0.55). This low numerical aperture scan lens was chosen to ensure that the PCOs in culture were within the depth-of-focus of the lens. A line clock was generated using a fiber Bragg grating with narrowband reflectivity at 990nm (AC Photonics) and a fast photodetector (Newport, 1811-AC-FC, 125 MHz bandwidth). The signal from the balanced detector (Newport, 1617-AC-FC, 800MHz bandwidth) and line clock photodetector were digitized using a high-speed digitizer (AlazarTech, ATS9350, 12 bit, 500 MS/s). The 90:10 and 50:50 fiber-coupled splitters were wideband (+100nm) with a center wavelength of 1064nm (Thorlabs, 90:10 TW1064R2A2A, 50:50 TW1064R5A2A). The circulators were wideband and polarization insensitive (AC Photonics). The beam was scanned using an XY galvanometer pair (Cambridge, 6210H). The scan waveform and frame trigger was output using an analog in/analog out data acquisition (DAQ) board (National Instruments, BNC-2110 and PCIe-6353). The computer used for acquisition was custom built with an Intel i7 quad-core processor, 32GB DDR3 RAM, GTX 550 Ti video card, and 240GB SSHDD. To eliminate specular reflection from the air-glass interface of the glass-bottom multiwell plate, the scan lens was mounted on a custom rotational mount to tilt the lens and XY galvanometers off-axis by 9° in the X direction.

Each well within a 24-well plate was imaged with a 6mm x 6mm lateral field-of-view (XY pixel size: 6µm). The volumes contain 1000 A-scans per B-scan, 1000 B-scans per volume, and 4 repeated B-scans per position. The A-scan length was 1337 pixels. Standard OCT processing was used to calculate intensity images [24] and the repeated B-scans were averaged for each position. Data was acquired pre-treatment, immediately post-treatment and at 12 hours, 24 hours, 36 hours, and 48 hours post-treatment.

The XYZ stage repositioning was automated with µManager [25]. OCT acquisition was performed using a custom LabView script. For each of the 2 plates, a total of five wells were imaged. Total imaging time per plate was less than 5 minutes, capturing 124 PCOs across five wells for P1 and 98 PCOs across five wells for P2 (Table 1). 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. 1(B)).

2.4 Quantitative image analysis

A framework for automated volumetric image analysis was developed to process OCT volumes and track each PCO in 3D across multiple time points (Fig. 1(C)). This extends previous work that developed 2D single-organoid tracking for autofluorescence images [26]. Linkage of each organoid over time was performed via a 3D centroid matching-based approach with single-particle tracking. This pipeline was implemented using MATLAB (2019b, Mathworks), except when specifically noted.

Preprocessing of each OCT volume was required before reliable segmentation and single-organoid tracking could be achieved (Fig. 1(C)). The goal here was to take OCT volumes that were captured off-axis and included the glass bottom of the multiwell plate and provide a flattened OCT volume that contains only organoids (Fig. 2(A),(B)). Figure 2(B) shows the image processing workflow for OCT volume image processing. First, the glass bottom was initially detected by calculating the Z gradient of the OCT volume and recording the coordinate of the maximum Z gradient value, which corresponded to the location of the air-glass interface. This provided a two-dimensional (2D) matrix of the location of the air-glass interface in Z for each XY location. A 3D second-order polynomial surface was then fit to the air-glass interface using least absolute residual (LAR) fitting (curvefit, MATLAB). A constant value based on the thickness of a #0 coverslip (n = 1.50 at 1040 nm) was added to the fitted coordinate values, which moved the coordinate values to the location of the glass-medium interface. These coordinate values were then used to shift each A-scan in the Z direction to yield a flattened volume that contained only organoids in the volumetric field-of-view. There was typically a large region in the volume in the Z direction that did not contain organoids, therefore each volume was cropped to contain only 200 pixels in Z. The final volume size was 1000 pixels by 1000 pixels by 200 pixels.

 

Fig. 2. Quantitative Image Analysis of OCT Volumes A quantitative image analysis workflow was developed that enabled volumetric growth tracking of PCOs using OCT. (A) Representative OCT volume of PCOs pre-treatment. (B) Tracking of volumetric growth of organoids was achieved using an image processing workflow shown here with OCT B-scans for clarity. The goal was to correct for the off-axis imaging geometry and provide a flattened OCT volume that contains only organoids. The air-media interface was detected in 3D and used to shift each A-scan in the Z direction to yield a flattened volume that contained only organoids in the volumetric field-of-view. A segmentation algorithm based on OCT intensity was then used to segment PCOs. (C) The automated segmentation algorithm was compared against manual segmentation. The Sørenson-Dice similarity coefficients for n = 3 pairs of automatically and manually segmented OCT volumes are plotted here as individual data points along with the mean and standard error. (D) Representative 3D single-organoid tracking results for an OCT time course. Each organoid was tracked over multiple time points, and the tracked centroid is displayed as an overlay to a 2D maximum intensity projection of pre-treatment OCT volume. Scale bar: 500µm.

Download Full Size | PDF

Registration of each OCT volume time series was achieved using both a registration mark and rigid image registration (Fig. 1(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. After acquisition, XY drift between frames was corrected through rigid image registration, which registers two frames by maximizing the cross-correlation of pixel values (dftregistration, MATLAB Central File Exchange) [27,28]. The maximum intensity projection (MIP) of the preprocessed OCT volumes at each time point was used to calculate the XY shift values needed to align the volumes.

PCO segmentation was performed on each OCT volume independently to capture changes in organoid volume (Fig. 1(C)). A segmentation algorithm based on OCT intensity segmented PCOs after OCT volume preprocessing. First, the PCOs were initially segmented using two-class k-means clustering (kmeans, MATLAB). Postprocessing of the segmented regions involved filling holes in the volumes (imfill, MATLAB). Then any segmented region smaller than 5 pixels by 5 pixels by 5 pixels was removed (bwareaopen, MATLAB) to set a lower bound on the organoid size of ∼1.5 × 10⁻⁵ mm³ (equivalent organoid diameter of 32 µm). Filling holes in the PCO volume was performed for PCOs with either a solid morphology (cellular core) or hollow morphology (non-cellular core) [26,29,30]. The volume and 3D centroid of each segmented PCO was measured for each time point (regionprops3, MATLAB). The voxel size (6 µm x 6 µm x 3.5 µm) was used to scale the volume measurements from pixels to mm³. The performance of the automated segmentation algorithm was compared against manual segmentation. OCT volumes were manually segmented and the Sørenson-Dice similarity coefficient was calculated for pairs of automatically segmented and manual segmented volumes (Fig. 2(C), n = 3, mean ± SEM = 0.89 ± 0.03).

Tracking of each PCO over time was achieved using an open-source single-particle tracking tool (Fig. 1(C), TrackMate, ImageJ/FIJI) [31]. 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 [32]. Volume time series were generated using MATLAB with 3D Gaussian particles (kernel: 10 × 10 × 10 pixels) at the 3D centroid of each segmented PCO. 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 (LAP) 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, 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. Once completed for an OCT volume time series, the single-organoid tracking results were used to link the organoid volume measurements over time.

2.5 Organoid volume estimation and comparison

One approach for measuring organoid volume involves diameter measurements from 2D images of organoids (e.g., brightfield microscopy) that are used in an equation to estimate the organoid volume [10,33,34]. This approach assumes that each PCO is spherical (Eq. (1)).

where d is the diameter of the 2D projection of each PCO. To create 2D images for organoid diameter measurements, first MIPs were calculated for each pre-treatment OCT volume. These provide similar information to brightfield microscopy images as both OCT and brightfield use light scattering as contrast. To automatically measure the diameter of each organoid, 2D PCO segmentation was performed corresponding to the 3D case described above. 2D segmentation was achieved using two-class k-means clustering (kmeans, MATLAB). Postprocessing of the segmented regions involved filling holes (imfill, MATLAB). Then any segmented region smaller than 5 pixels by 5 pixels was removed (bwareaopen, MATLAB) to set a lower bound on the organoid size of approximately 32 µm diameter. The diameter was then quantified for each 2D segmented organoid (regionprops, MATLAB), and these diameters were used to calculate the estimated organoid volume (Eq. (1)). The estimated organoid volume (Eq. (1), “2D Est.”) was then compared to the measured organoid volume (regionprops3, MATLAB, “OCT”) for each organoid at pre-treatment, and the volume estimation error between the two values was quantified by calculating the percent error ($|{({OCT - 2D\; Est.} )/OCT} |\ast 100$).

Sphericity is a measure of how close a 3D volume is to a sphere and is based on variables (volume and surface area) that are easily measured for each segmented PCO [35]. Measuring organoid sphericity for each PCO enables a quantitative investigation of the error between direct volume measurement using OCT and volume estimation from 2D diameter measurements for organoids with different shapes.

where V and SA are the organoid volume and surface area measured using 3D segmentation of each PCO from OCT volumes (regionprops3, MATLAB).

2.6 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 [36,37]. Open-source R toolboxes included nlme, a mixed effect modeling package, and lsmeans, a least-squares means estimation package [38,39].

Organoid-level-normalized time series were calculated for each organoid as X_{post-treatment}/X_{pre-treatment_organoid}, where X is the organoid volume. Linear mixed-effect models were used to analyze these organoid-level-normalized time series for each treatment and PCO line, specifically a longitudinal random-intercept model was used (nlme, R) [38]. Random-intercept models compensate for organoid variability at pre-treatment, which enables a more accurate assessment of individual PCO drug response [40,41]. A random-intercept model was specified as $Y\sim \; Time + Treatment + Treatment\ast Time + ({1\textrm{|}Organoid} )+ Error$, where Y is the pre-treatment-normalized organoid volume, 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 minutes post-treatment, 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 organoid volume. Organoid here is encoded as a random effect to account for variability at the organoid-level. Pre-treatment values were excluded from the analysis because all time points were normalized to pre-treatment values (i.e., pre-treatment variance is one). 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. Pairwise comparisons were performed between all treatment groups at each time point using least-squares means (lsmeans, R) [39]. Tukey’s Honest Significant Difference method was used to account for multiple comparisons.

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 pre-treatment mean as X_{post-treatment}/X_{pre-treatment_well}, where X is the organoid volume. A fixed effects model was specified as $Y\sim \; Time + Treatment + Treatment\ast Time + Error$, where Y is the well-level-normalized organoid volume, 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. Pre-treatment 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.

The difference between OCT-measured organoid volume and 2D-estimated organoid volume was tested using the non-parametric Wilcoxon Signed Rank test. The correlation between volume estimation error and OCT-measured organoid volume was quantified using Spearman’s correlation coefficient. The correlation between volume estimation error and organoid sphericity was quantified using Spearman’s correlation coefficient.

For the growth rate heterogeneity analysis, single-component Gaussian distributions were fit to the pre-treatment normalized organoid volume at 48 hours post-treatment for each of the treatment groups. Single-component Gaussian were chosen from visual inspection of the data, which was unimodal and normally distributed for each treatment group. Then, PCOs with hollow morphologies were identified via visual inspection. The Z-score [($X - {\mu _{Treatment}})/{\sigma _{Treatment}}]\; $for all PCOs were calculated and stratified into hollow or solid morphology groups, where X is the volume of an organoid in a treatment group, µ_Treatment is the mean volume of all organoids within that treatment group, and σ_Treatment is the variance of all organoids within that treatment group. The Z-score enables comparison between treatment groups. The non-parametric Mann-Whitney U test was used to test the difference between the hollow and solid morphology groups.

3. Results

3.1 Drug response reduces organoid volumetric growth rate

Drug response changes PCO volume or volumetric growth rate. Figure 3(A) shows representative MIP of individual PCOs imaged over time, including one each from the P1 and P2 control (no treatment) groups and one from the P1 cyanide treated group. Single-organoid tracking was used to track the volumetric growth of PCOs in response to treatments, which were compared against PCOs grown in a control well. In addition, heterogeneity exists in both organoid morphology and size. Again, singe-organoid tracking enabled normalization of each PCO volume to its value pre-treatment, so that relative changes in organoid volume can be compared across organoids of different sizes. Linear mixed models were fit to the pre-treatment-normalized organoid volume time-series of all organoids from a PCO line (P1 or P2). The change in organoid volume caused by drug treatment was quantified for all five treatment groups. Least-squares means was used to calculate the pairwise differences between each treatment group and control for each PCO line (Fig. 3). Figure 3 also shows the pairwise differences between all treatment groups for P1 and P2.

 

Fig. 3. Treatment-induced Changes in Patient-Derived Cancer Organoid Volume. Single-organoid tracking quantified changes in organoid volume for each treatment group of PCOs over time. (A) Representative MIP of individual PCOs imaged over 48 hours, including one each from the P1 and P2 control groups and one from the P1 cyanide treated group. The 20-minute post-treatment time point was omitted for clarity as there was no change in PCO volume. (B) Heatmaps show the pairwise differences between pre-treatment-normalized organoid volume for all treatment groups for P1. Top: Pairwise differences between each treatment and control for P1. Bottom: Pairwise differences between drug treatments (excluding control) for P1. (C) Heatmaps show the pairwise differences between pre-treatment-normalized organoid volume for all treatment groups for P2. Top: Pairwise differences between each treatment and control for P2. Bottom: Pairwise differences between drug treatments for P2. All pairwise differences were calculated from the linear mixed-effect models via least-squares means and are significant at the p < 0.05 level. Black boxes indicate that these differences are not significant (p > 0.05, black boxes). 2DG: 2-deoxyglucose, Pac: paclitaxel, Cis: cisplatin. Table 1 shows the number of PCOs in each treatment group for each patient. Scale bar: 50µm.

Download Full Size | PDF

PCOs from P1 and P2 had similar drug response patterns. For P1 groups treated with metabolic inhibitors (cyanide, 2-deoxyglucose), lower organoid volumes compared to control were observed starting at 12 hours for the cyanide (p < 0.001) and 2DG (p < 0.001) groups (Fig. 3(B)). P1 PCOs treated with cisplatin or paclitaxel were significantly different from control beginning at 12 hours (p < 0.01, Fig. 3(B)). For P2 PCOs treated with metabolic inhibitors, lower organoid volumes compared to control were observed starting at 12 hours for the cyanide (p < 0.001) and 2DG (p < 0.001) groups (Fig. 3(C)). For P2 PCOs treated with cisplatin or paclitaxel, organoid volume was lower than control starting at 12 hours (p < 0.01, Fig. 3(C)). Heatmaps (Fig. 3(B),(C)) also show the pairwise differences between treatment groups.

In both P1 and P2 control groups analyzed using single-organoid tracking, PCOs continue to grow over 48 hours, with a statistically significant increase in organoid volume of 32% and 30%, respectively, at 48 hours compared to pre-treatment (Fig. 4(A),(C) p < 0.001 pre-treatment vs. each time-point post-treatment using organoid-level normalization).

 

Fig. 4. 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) PCOs. Normalization refers to dividing the value of each organoid at each time-point by either its own pre-treatment value (organoid-level normalization) or the well-level mean (well-level normalization). (A,B) Pre-treatment normalized organoid volume for P1 using organoid-level normalization (A) and well-level normalization (B) for each treatment group. (C,D) Pre-treatment normalized organoid volume for P2 using organoid-level normalization (C) and well-level normalization (D) for each treatment group. Data in A,C were analyzed using a linear mixed-effect model that uses single-organoid tracking to account for organoid-level variability. Data in B,D 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 is mean (lines) ± standard error of the mean (shaded area). Significant differences between a treatment and control (p < 0.05) are indicated with a square color-coded for the treatment group (see legend). Table 1 shows the number of PCOs in each group and for each patient.

Download Full Size | PDF

3.2 Single-organoid tracking increases sensitivity to drug response

The volumetric growth of each PCO was quantified using single-organoid tracking. Pre-treatment normalization of each organoid to its initial volume accounts for heterogeneity before treatment, and thereby provides more robust analysis of treatment effects. A comparison between single-organoid tracking and a pooled analysis was performed. For the single-organoid tracking, the volume for each PCO was normalized to the pre-treatment volume of the same PCO. For the pooled analysis, each PCO was normalized to the well-level mean of all PCOs in that well. Figure 3 shows organoid volume for P1 and P2 analyzed using single-organoid tracking (organoid-level normalization) or pooled analysis (well-level normalization). Single-organoid tracking (Fig. 4(A),(C)) yields lower variability compared to pooled analysis (Fig. 4(B),(D)), as the well-level mean does not account for organoid-level variability. Importantly, the statistical model that can be applied to organoid-level time series provided by single-organoid tracking increases statistically significant differences between treatment. Overall, single-organoid tracking improves sensitivity to drug response compared to pooled analysis.

3.3 Optical coherence tomography directly quantifies organoid volume

The PCO diameter in 2D images, typically brightfield microscopy images, are often used to estimate changes in PCO volume. This estimation procedure assumes that each PCO is spherical, and therefore large errors can occur with irregular shaped organoids. OCT enables direct measurement of PCO volume, which eliminates the need for volume estimation. Figure 4 compares 3D measurements of PCO volumes using OCT to 2D volume estimates using OCT 2D maximum intensity projections (MIPs) that mimic brightfield measurements. A graphical illustration and equation for volume estimation based on 2D diameter measurements in OCT MIPs is shown (Fig. 5(A)). Dot plots of the OCT measured and 2D estimated PCO volumes are shown for both P1 and P2 (Fig. 5(B)). The distributions of OCT volumes are log-normal, where there are more small organoids than large organoids, which agrees with qualitative observations of the sizes of PCOs present in the culture. Note that organoids smaller than 32 µm diameter are excluded from this analysis.

 

Fig. 5. Quantification of Volume Estimation Error. Percent error between direct volume measurement using OCT and volume estimation from diameter measurements of 2D maximum intensity projections (MIPs) of OCT volumes that mimic brightfield measurements. (A) Illustration of the estimated organoid volume (V) based on a diameter measurement (d) in a 2D MIP of an OCT volume. (B) Dot plots of organoid volumes either directly measured from 3D OCT (“OCT”) or estimated from 2D MIPs of OCT volumes (“2D Est.”, as shown in (A)). Statistical significance was assessed with the non-parametric Wilcoxon Signed Rank test (** for p < 0.01; n.s., not significant). Mean (dotted line) and 95% confidence interval (solid error bars) are shown. (C) Plot of the volume estimation error as a function of OCT-measured organoid volume for all organoids pre-treatment. Volume estimation error is the percent error ($|(OCT - 2D\;EST.)/OCT|*100$) between the OCT-measured organoid volume and 2D-estimated organoid volume (as shown in (A)). (Spearman’s ρ = 0.37, p < 0.01) (D) The mean intensity projection (top-down view) of each representative PCO and a B-scan chosen from the center (denoted with a red dashed line) of each organoid. These PCOs were chosen to be representative of low (0.42), medium (0.65), and high sphericity (0.93), see Eq. (2). Scale bar: 50 µm. (E) Plot of the volume estimation error for the three representative PCOs (shown in d) at pre-treatment and multiple time points post-treatment. The high sphericity organoid had the lowest error (9% at 48 hours) and low sphericity organoid had the highest error (38% at 48 hours post-treatment), with the medium sphericity organoid in between (25% at 48 hours post-treatment). (F) Plot of the volume estimation error as a function of organoid sphericity for all PCOs at 48 hours post-treatment (Spearman’s ρ = -0.85, p < 0.01). Table 1 shows the number of PCOs for each patient.

Download Full Size | PDF

Figure 5(C) quantifies the volume estimation error (percent error) between the OCT-measured volume and 2D-estimated volume for every PCO pre-treatment (n = 124 PCOs for P1, and n = 98 PCOs for P2). A linear model was fit to the volume estimation error as a function of measured PCO volume. There was a moderate positive correlation (Spearman’s ρ = 0.37, p < 0.01) between the volume estimation error and the measured PCO volume, which shows that estimating PCO volume from 2D diameter measurements becomes increasingly inaccurate with larger PCOs.

Next, the relationship between PCO morphology (sphericity, see Eq. (2)) and the volume estimation error (percent error) as a function of time after treatment was examined. Figure 5(D),(E) shows three representative PCOs and their volume estimation error as a function of time. In addition, Fig. 5(F) shows the volume estimation error at 48 hours as a function of sphericity. Three representative PCOs were chosen to be representative of low, medium, and high sphericity (Fig. 5(D)). Figure 5(D) shows the mean intensity projection (top-down view) of the entire PCO and a representative B-scan chosen from the center of the PCO (denoted with a red dashed line) for each PCO. Figure 5(E) shows the volume estimation error for each of these representative PCOs at pre-treatment and at multiple time points post-treatment. As expected, the high sphericity PCO (sphericity = 0.93) has the lowest volume estimation error (9% at 48 hours post-treatment), whereas the PCO with the lowest sphericity (sphericity = 0.42) has the highest error (38% at 48 hours post-treatment). The medium sphericity PCO (sphericity = 0.65) had a volume estimation error (25% at 48 hours post-treatment) that fell between the other two PCOs. To see if this observation held for all PCOs used in the study, Fig. 5(F) shows the volume estimation error at 48 hours post-treatment as a function of organoid sphericity for all PCOs from both P1 and P2. There is an inverse relationship between volume estimation error and organoid sphericity (Spearman’s ρ = -0.85, p < 0.01), with higher organoid sphericity resulting in lower volume estimation errors.

3.4 Patient 1 organoids exhibit growth rate heterogeneity

Additional analysis was performed to explore heterogeneity in PCO growth rate. Patient 1 had PCOs with solid or hollow morphologies (Fig. 5(D), middle and bottom panels, respectively). Patient 1 had 103 PCOs with a solid morphology and 21 with a hollow morphology. Patient 2 only had PCOs with solid morphologies (Fig. 5(D), top panel). Figure 6(A) shows single-component Gaussian distribution fits to the pre-treatment normalized PCO volumes for each of the Patient 1 treatment groups at 48 hours. Stratifying the Patient 1 PCOs into hollow and solid sub-types qualitatively showed that the hollow PCOs in each treatment group had higher growth rates relative to the treatment group mean (Fig. 6(B)). Further quantitative analysis was performed by calculating the Z-score for each PCO relative to all other PCOs in that treatment group (Fig. 6(C)). PCOs with positive Z-scores had growth higher than the mean of all other PCOs in that treatment group, while PCOs with negative Z-scores had growth lower than the mean of all other PCOs in that treatment group. This analysis quantitatively shows that hollow PCOs had higher pre-treatment normalized organoid volume at 48 hours compared to solid PCOs (Fig. 6(C), p < 0.01).

 

Fig. 6. Volumetric Growth Rate Heterogeneity. Heterogeneity analysis of pre-treatment normalized organoid volume at 48 hours for Patient 1 (P1) PDCOs. (A) Single-component Gaussian distributions were separately fit to the pre-treatment normalized PCO volumes for each of the P1 treatment groups at 48 hours. (B) Dot plots of the pre-treatment normalized organoid volume at 48 hours for PCOs in each of the P1 treatment groups. PCOs were stratified into hollow and solid morphology sub-types based on visual inspection of the images. Hollow organoids for each group are plotted using open circle markers and solid organoids are plotted using solid square markers. CN: cyanide, Cis.: cisplatin, Pac.: paclitaxel. (C) The Z-score [($X - {\mu _{Treatment}})/{\sigma _{Treatment}}]\; $for all pre-treatment normalized P1 PCOs at 48 hours were calculated and stratified into hollow or solid morphological sub-types. The Z-score enables comparison between treatment groups. Statistical significance was assessed using the non-parametric Mann-Whitney U test (** for p < 0.01). Distributions of the Z-scores were plotted from the dot plot data (gramm Package, MATLAB). Table 1 shows the number of PCOs for Patient 1.

Download Full Size | PDF

4. Discussion

Clinicians base the selection of a treatment for an individual cancer patient on biomarkers (e.g., mutational profile, medical imaging, histology, and stage) with limited ability to predict response from heterogeneous populations of tumor cells. Drug screening with PCOs is a promising approach to optimize treatment decisions for individual cancer patients. However, non-destructive techniques to quantify drug response on an individual PCO level are needed to capture intra-patient heterogeneity. Here, we provide the first demonstration of 3D single-organoid tracking using OCT to assess volumetric growth and drug response in PCOs. Single-organoid tracking enables organoid-level measurements of drug response over time that improves sensitivity over bulk measurements across an entire well. OCT is inherently non-destructive, rapid, and label-free, which is attractive for capturing heterogeneity in drug response across numerous PCOs from the same patient.

PCOs typically vary in size and morphology due to the 3D culturing technique, so single-organoid tracking is needed to account for this baseline heterogeneity [42,43]. OCT has been used to image drug response in organoids by tracking volumetric growth, measuring the amount of cell death, or quantifying changes in cellular-level dynamics [10–12,14]. None of these previous OCT studies tracked each organoid over time, but instead measured a single time point or calculated averages over all organoids in a well. Here, we demonstrate a framework to link the same PCO across multiple OCT imaging time points. These methods can be used to measure drug response in PCOs that accounts for baseline PCO variability. In addition, statistical techniques such as mixed-effect models can use these time-course data to account for organoid-level heterogeneity and provide improved sensitivity to drug response [40,44]. Our results show that single-organoid tracking of PCOs reduces variability and increases statistically significant differences between the treatment groups.

To demonstrate the benefit of OCT for drug response measurements in PCOs, we compared OCT direct volume measurements to 2D volume estimation that mimics standard brightfield measurements. Interestingly, while there was error for organoids from P1 and P2, the distribution of actual and estimated PCO volumes were only significantly different for PCOs from P2. Further examination of volume estimation error as a function of sphericity showed that error decreased as sphericity increased. PCOs from P2 had lower sphericities, which may explain the difference in the volume distributions between OCT-measured volume and 2D-esimtated volume. Our results show that assessing the volumetric growth rate of PCOs offers improved accuracy compared to traditional 2D estimates similar to brightfield diameter measurements.

Lastly, P1 PCOs exhibited two morphological sub-types, solid and hollow, which were identified by visual inspection of the OCT images. The hollow morphology is consistent with cystic morphologies observed in previous studies of colorectal PCOs, and is thought to reflect the genetic and phenotypic in vivo heterogeneity of the patient’s tumor [30]. Hollow PCOs grew more over 48 hours compared to solid PCOs, which may be a result of differential gene expression between the two morphologies or the lack of an apoptotic core [29,30]. Note that differences in the density of hollow compared to solid organoids were not taken into consideration in the volume estimation, which assumes that all organoids are dense and non-cystic. Future studies could set a threshold on the OCT signal to only quantify the cellular regions of cystic organoids in the volume calculations. Our results show that 3D single-organoid tracking and OCT can be used to identify growth rate heterogeneity, and differences in growth rate between hollow and solid PCOs.

This study illustrates a framework for robust drug response measurements of PCOs using OCT and 3D single-organoid tracking that accounts for pre-treatment heterogeneity in PCO volume. OCT can provide rapid volumetric images of organoids, which could be used to characterize diverse drug response across numerous PCOs derived from the same patient. The framework developed in this manuscript can be extended to evaluate other aspects of the PCO system, including basic scientific questions about the growth dynamics and heterogeneity of PCOs.

In addition, 3D single-organoid tracking could be further developed to enable multimodal optical imaging of the same PCOs over time. We envision that multiple optical modalities, such as autofluorescence imaging [26,29,45–50] and functional extensions of OCT [11–13,51–55], could be used with 3D single-organoid tracking to monitor each PCO over time. Autofluorescence imaging of the metabolic coenzymes NAD(P)H and FAD could provide label-free metabolic information of each organoid and has already been demonstrated to be a powerful tool for assessing PCOs [26,29,45–50]. Fluorescence lifetime of NAD(P)H and FAD reports on the binding activity of these enzymes, and can provide additional information on the specific metabolic program of cells [47,50,56–59]. Images taken using other modalities at higher magnification [60–62] could provide cellular and sub-cellular information for each PCO, which could then be integrated with OCT to provide a comprehensive readout of organoid health and drug response. Images at multiple scales could also provide insight into the relationship between the OCT signal in PCOs and its corresponding cellular and sub-cellular organization. Lastly, because OCT is non-destructive, methods such as genomic sequencing, transcriptomics, proteomics, and metabolomics could be used in tandem to provide complementary information on drug response [2,8,63,64].

We believe that OCT coupled with 3D single-organoid tracking can be applied to improve clinical treatment selection, basic science, and drug development.