We established an HCS pipeline that includes label-free temporal imaging of PDOs ( Figure 1A ) coupled with quantitative image analysis using a linear classifier ( Figure 1B ) and data compilation and visualization ( Figure 1C ). To execute this workflow, we utilized PDOs from our biobank generated from CRC patient samples, which includes primary and liver metastatic tumors.

Organoid set up ( Figure 1A ): Briefly, to generate PDOs that recapitulate the morphology of the tissue of origin ( Supplemental Figure 1 ), tissue samples were obtained post-surgery, processed, seeded in extracellular matrix, and expanded for future use (see Methods section for details). To set up the screening assay, organoids were first digested to a single cell suspension before being seeded into a 96 well plate. Then, using an HCS platform, the samples were imaged at multiple timepoints in brightfield to minimize phototoxicity and photobleaching.

Supervised machine learning algorithm used to classify organoids as live/dead based on phenotypic features ( Figure 1B ): Image analysis was subsequently performed on the maximum intensity projections of multiple z-scan images using a machine learning algorithm that enables users to build a linear classifier by identifying regions of interest (ROIs) that are part of distinct groups. Using a trained feature-based textural machine learning algorithm we divided image regions into two classes: organoid ROIs and background (the ML algorithm was trained on the segmented and annotated images, see the Methods section for details). Morphological and textural features were measured for each identified object within the organoid class ( Supplemental Table 1 ). STAR morphology features encompass Symmetry properties, Threshold compactness, Axial and Radial properties. Spot-Edge-Ridge (SER) textural features are based on Gaussian derivative images measuring pixel intensity patterns within each ROI. Distributions of all 25 STAR and SER features measured across 6 different PDOs, in media or treated with staurosporine (positive control - apoptosis inducer), are depicted in Supplemental Figure 2 (see PDO counts and tissue site in Supplemental Table 2 ). At day 3, morphological features are patient-specific and consistent across replicate wells ( Figure 2A ) and unsupervised clustering of the data identified clusters that matched the patient or origin rather than number of days in culture ( Supplemental Figure 3 ). Using this ML approach, we can detect morphometric similarities and differences across PDOs.

Statistical analysis and development of data visualization tool ( Figure 1C ): For all 25 textural and morphological features, the measurements for each detected PDO were first summarized (mean value) across all detected objects within a well. Next, the mean and standard deviation were computed from technical replicate wells for each treatment and time-point. A Shiny-based web tool, the Organoizer, was developed to process the data and produce plots for organoid-survival and monitor changes in features over the course of the treatment period (27).

To create a training data set, we manually classified 179 objects, 80 live (untreated media control) and 99 dead (5 μM staurosporine treated positive control), across 41 images of organoids derived from six different patients at various time points. When applied to new experimental data each detected PDO was assigned to either the “dead” or “live” class by the linear classifier based on 9 significant morphology and texture features chosen by the algorithm to delineate between live and dead PDO categories: SER Valley, SER Edge, SER Ridge, Profile 2/2, Threshold Compactness 60%, Axial Small Length, Area, Ration Width to Length, and Threshold Compactness 50% ( Figure 2B ). Representative images illustrate selected PDO textural and morphological features (i.e., Profile 1/2, and 2/2, SER Edge Ridge and Valley), found to be distinct between live and dead PDOs ( Figure 2C ). The signal to noise ratio of the classifier algorithm is expressed as “goodness” based on the distance of the data points from the classifier line, which is visualized using a scatter plot ( Figure 2D ).

As a baseline, we compared our classifier to the visual assessments of trained cell biology experts. Trained cell biologists are adept at visually assessing the health of their cell cultures using brightfield microscopy, however manual classification not only limits throughput but also introduces inter-observer variability (23). To generate a ground truth by visual assessment, we asked 9 scientists to blindly classify images of individual organoids (18 PDOs) as either “live” or “dead” ( Supplemental Figure 4 ). Statistical evaluation of inter-rater reliability indicated only moderate agreement between visual classifications (Fleiss’ κ = 0.451, z = 11.5, p = 0.00; perfect agreement κ=1), highlighting the inherent variability in subjective manual classification. Using the data set containing the 18 manually classified organoids, we also performed live/dead classifications based on vital dye (VD) intensity from DRAQ7 staining. For our purposes we defined a dead PDO as one that contained at least one DRAQ7⁺ area ≥ the area of a nucleus. For each organoid we compared the expert consensus classification against DRAQ7 staining results and our linear classifier ( Figure 3A ). We found 78% (14/18) concordance between the linear classifier and expert majority, 61% (11/18) between the linear classifier and DRAQ7, and 61% (11/18) between the expert majority and DRAQ7. In instances where all experts agreed, concordance between the linear classifier and the expert majority increased to 100% (8/8), however expert majority concordance with DRAQ7 only reached 62% (5/8). The strong agreement between the expert classifications and the linear classifier reinforces machine learning as a valuable approach for 3D organoid phenotyping. An important note, the 18 PDO images classified across methods (i.e., experts/ML/VD) were not included in the training set to ensure that there is no “leakage” of information between the training and the testing set.

To further evaluate the performance of the linear classifier using time series data, we compared our algorithm classifications with those made using DRAQ7 for PDO-12620 ( Figure 3B ). We normalized the proportion of live/total PDOs at each timepoint to the proportion of live PDOs at day 0 for untreated and staurosporine treated conditions determined by ML or DRAQ7 staining. For the staurosporine treated group, the ML classifier detected a reduction in the number of live organoids on day 1, whereas DRAQ7 shows a comparable reduction past day 3 ( Figure 3C ). In the untreated group, many organoids are classified as live by both the ML and the VD, with 80% concordance between the two methods. However, the two methods started to diverge upon treatment, with concordance in staurosporine treated organoids dropping to 60% ( Figure 3D ). Additionally, our ML approach allows us to follow individual organoids over time to determine their vital status ( Figure 3E ). To further clarify the discrepancy between the ML-based and DRAQ7-based classifications we tracked individual PDOs (PDO-13154, N~900) over 7 days (days 0,1,3 and 7) and determined their vital status at each time point by each method ( Supplemental Figure 5 ). While the ML classification of dead PDOs increased with time, the population of dead PDOs determined by DRAQ7 decreased. The discordance is most likely due to clearance of cellular debris over time where a PDO previously defined as dead is now DRAQ7 negative. Given staurosporine is an inducer of apoptosis, this result suggests that our ML method may identify dead or dying PDOs more accurately including those that have lost their ability to retain DRAQ7.

Tumors evolve over time in response to various stimuli, such as organ-specific microenvironments and drug perturbations. Our approach allows us to characterize the dynamic drug responses of PDOs from both primary and metastatic CRC tumors over time. To accomplish this, we treated PDOs with standard chemotherapy agents: irinotecan, a topoisomerase I inhibitor, and oxaliplatin, an alkylating agent. To interrogate drug specific phenotypic responses, we used heatmaps to examine morphological and textural features within the dead class of PDOs over the course of drug treatment. ( Figure 4A ). Across all PDOs, the feature pattern in the staurosporine-treated group is distinct from the chemotherapy groups. For PDO-12415, the symmetry features in the media control stand out with generally higher feature values compared to the drug treated groups. PDO-12527 and PDO-12911 showed an increase in area and radial mean over time for all treatments, however the symmetry properties of all the PDOs did not show distinct variation across treatment or time. Importantly, the increase in area and radial mean could be attributed to the loss of structure and spreading of dead organoids in response to drug rather than proliferation.

Using our Shiny-based visualization tool we generated box plots of extracted features of interest over time ( Figure 4B and Supplemental Figure 6 ). We chose to highlight Regional Threshold Compactness 60% due to the unique patterns observed across patients; however, our visualization tool is capable of displaying all features captured. In addition, we plotted dose response curves for each drug treatment ( Figure 4C ). While PDO-12415 (bottom left panel) showed a limited response to irinotecan at 20 μM, a much stronger effect was measured with 40 μM irinotecan starting on day 1. Oxaliplatin at 20 μM also elicited a response with the normalized proportion of live/total dropping 60% by day 7. Analysis of PDO-12527 (bottom middle) revealed a similar response to both oxaliplatin and irinotecan at 20 μM, while irinotecan at 40 μM more effectively killed the PDOs. It appeared that PDO-12911 (bottom right) did not respond to oxaliplatin at 20 μM, as the proportion of live/total PDOs was comparable to the negative control across all 7 days. A slight difference in response timing was seen between the two doses of irinotecan, where the 40 μM dose showed a stronger response at day 3 compared to 20 μM, however by day 7 both showed a response below 40%. Despite temporal variations in response, all three PDOs showed a 60-70% reduction in the proportion of live/total PDOs by day 7. Taken together our ML approach identified PDOs that responded to chemotherapy early in the dosing regimen, highlighting the ability to capture patient-specific drug responses.

Organoid growth medium consists of base medium (ADMEM/F12 with 10% FBS, 1% penicillin/streptomycin, 1% Glutamax, and 1% HEPES) supplemented with 1X N2 (Sigma Aldrich, 17502048), 1X B-27 (Sigma Aldrich,17504044), 1mM N-Acetylcysteine (Sigma Aldrich, A7250) 50 ng/ml EGF (Life Technologies, PGH 0313), 100 ng/ml Noggin (Tonbo, 21-7075-U500), 10 mM nicotinamide (Sigma, N0636), 500 nM A 83-01 (Calbiochem, 616454-2MG), 10 μM SB202190 (Sigma 47067), and 0.01μM PGE2 (Sigma Aldrich, P5640).

Tissue digestion solution consists of 1.5 mg/ml collagenase (Millipore, 234155), 20 μg/ml hyaluronidase, (MP Biomedicals 100740) and 10 μM Ly27632 (Sigma Y0503).

Tumor tissue was received from consented patients following Institutional Review Board (IRB) approval at the Norris Comprehensive Cancer Center of USC, Los Angeles CA. Tissue was washed with PBS, minced and digested for 30 minutes at 37°C. Digest suspension was filtered using a 100 μm strainer to remove large residual pieces of tissue, then centrifuged at 189 x g for five minutes. Pellet was washed in DMEM/F12 media (ThermoFisher, 11320033) supplemented with 10% FBS three times and single cells were re-suspended in BME (Culturex® Reduced Growth Factor Basement Membrane Matrix Type 2, Trevigen, 3533-005-02). Cell/BME mixture was plated in 24 well plates with 60 μl per well and incubated upside down at 37°C until solidified (10-20 minutes). Then 500 μl of organoid growth media was layered on top and media was changed as needed. To passage organoids, BME was dissociated with 500 μl/well TrypLE (ThermoFisher Scientific, 12605028) for 5 minutes at 37°C. Organoid suspension was pooled and centrifuged at 450 x g, the pellet was re-suspended in BME and re-plated in a 24 well plate. PDOs used for experiments were ≤ 20 passages in culture.

PDOs were harvested from BME using Gentle Cell Dissociation Reagent (Stemcell technologies, 07174), pooling all wells and incubated on ice for 45 minutes then centrifuged for 5 minutes at 189 x g. Supernatant was removed and the pellet was re-suspended in 50% TrypLE with 10 µM Y-27632 (Stemcell Technologies, 72302), incubated at 37°C for 10-15 min with occasional agitation. Alternatively, PDOs were harvested using 500 μl/well TrypLE, incubated at 37°C for 30 minutes. For both methods, PDOs were centrifuged for 5 minutes at 189 x g. Supernatant was removed and the pellet was re-suspended in 1 ml of organoid base medium then filtered through a 40 μm strainer to remove aggregates. Flow through was centrifuged at 189 x g for 5 minutes and the pellet was re-suspended in BME. A 96 well μ-Plate (Ibidi, 89646) was coated with 5 μl of BME/well and incubated at 37°C until BME solidified. The μ-plate was then seeded with 5 μl BME/cell mixture at a concentration of 1000 cells/μl and topped with 70 μl organoid growth medium.

Organoid ROIs were counted, and ROI-level morphological metrics were averaged on a per-well basis at each timepoint and for each class (“dead” vs. “live”). The mean and standard deviation were then computed from replicate wells with the same treatment conditions. Response curves were computed as either the proportion of “live’ ROIs over “dead” ROIs (ratio) or as “live” ROIs over all ROIs (proportion). Optionally, the proportion or ratio of “live” organoids can be normalized to the proportion (or ratio, respectively) on the first day of measurement (usually day 0). Boxplots for each feature across timepoints were also generated. The code is available at https://github.com/eitm-org/organoid_drug_response.

Heatmaps were generated by averaging feature values per well, then taking the average value across wells to get one average value for each unique group. Rows and columns were grouped using hierarchical clustering and rows were scaled using the heatmap package in the R statistical computing language. All analyses were performed using the R statistical language (v. 4.1.0) using the following packages: cowplot (v.1.1.1) (35), eulerr (v. 6.1.1) (36), ggridges (v. 0.5.3) (37), ggthemes (v. 4.2.4), here (v. 1.0.1), irr (v. 0.84.1) (38), knitr (v. 1.33), networkD3 (v. 0.4), pheatmap (v. 1.0.12), plater (v. 1.0.3), readxl (v. 1.3.1), reshape2 (v. 1.4.4), scales (v. 1.1.1), tidyverse (v. 1.3.1) and viridis (v. 0.6.1) (27, 39–49).