3D Morphological Characterization of Human Cortical Organoids Using a Customized Image Analysis Workflow

1. Introduction

The advent of research using human induced-pluripotent stem cells (hiPSCs) has seen personalised models of neurodevelopmental disorders made possible for the first time [1,2].

HiPSCs multiply continuously in monolayer cultures while retaining the ability to differentiate into nearly any human cell type [1,2]. They originate as mature somatic cells, commonly fibroblasts obtained from skin samples, and are converted into pluripotent stem cells using reprogramming factors [1,2,3]. Several different approaches are available for differentiation, including dual-SMAD inhibition or the overexpression of transcription factors such as Neurogenin-2 [4].

Derived from hiPSCs, three-dimensional (3D) brain or cerebral organoids have emerged as a ground-breaking technology with the potential to combine traditional laboratory findings with personalised clinical treatments [5,6]. These self-organising 3D structures closely mirror the neuronal organisation and functional activity of the developing human brain [2,5,6]. Depending on the protocol applied, they may represent cells from all brain regions, like cerebral organoids or the patterning of specific regions, such as the cortex, yielding so called corticoids. The ability to produce cortical-specific organoids is crucial for investigating disease mechanisms of severe developmental disorders like developmental and epileptic encephalopathies (DEE) [7,8], Autism spectrum [9,10] and Angelman syndrome [11]. In addition to pathophysiology of neurodevelopmental disorders, these models can be used for therapeutic screening in a tissue-relevant context.

Despite the promise of organoid technology, small changes made by individual laboratories when implementing differentiation protocols contribute to structural heterogeneity in organoid models. Previous efforts to standardize organoid differentiation have aimed to reduce variability [12], but any adopted protocol still requires minor adjustments based on specific cell lines and procedures. Notably, the size, shape and viability of cortical organoids can vary significantly both within and between batches. To address this heterogeneity, we sought to create an image analysis workflow specifically designed to examine large datasets and benchmark the selected protocol. This imaging workflow would serve as an important exercise in standardization.

Evidence of cellular diversity in organoid models has primarily been drawn from measures of cell function, gene expression and protein abundance; however, these give little insight into architectural patterns of cell growth [13,14]. Recently many groups have focused on high-throughput imaging methods to assess organoid morphology [15,16,17,18]. However, at the cellular level, high-resolution imaging techniques are still needed to measure the integrity of maturing organoids, which is challenging when applied to cortical organoids [13,14]. Additionally, the need for large sample datasets to obtain statistically valid comparisons creates a bottleneck in the imaging process. High-content imaging of organoid models offers a promising solution to accelerate tissue analysis, especially as these models become more structurally complex as they mature.

Brain organoid maturation facilitates synaptic and network development, necessitating extended culturing periods, which ultimately results in larger tissue sizes [19,20]. Investigating the cellular architecture of these larger tissues often involves sectioning tissue and employing immuno-histochemistry (IHC), a widely used procedure for organoid analysis [14]. Other methods, such as tissue clearing, 2-photon imaging, and light-sheet imaging, require specialized expertise and imaging systems, making them less commonly used [14,21,22]. Our objective therefore, was to develop a widely applicable image analysis method that could benefit various research groups.

One major challenge associated with large avascular organoids is their susceptibility to developing a necrotic core. A comprehensive structural analysis method must include a measurement approach that filters out dying cells and the necrotic core region [23]. The absence of a vascular system limits the diffusion of oxygen and nutrients from the surface of organoids, leading to the formation of a dead cell core as the organoids grow [24,25]. Despite the emergence of methods to improve nutrient diffusion, physical constraints limit nutrient flow to approximately 400 µm and in long-term 3D organoid cultures, hypoxic core susceptibility remains an issue [26,27,28]. Changes in the volume of dead cell cores or the percentage of dead cells may also reflect features of the disease model [29,30,31]. Image-based analyses have traditionally relied on manual delineation of non-viable tissue boundaries, a process that is subjective, time consuming and impractical for large sample sizes.

Another critical measurement we sought to address was obtaining cell counts for proteins expressed outside the nucleus or in peri-nucleic regions. While robust co-localization imaging methods with cell nuclei are available, antibody labelling without spatial overlap – labelling non-nuclear proteins – presents a significant analytical challenge. Assigning these non-nuclear structures to a specific cell is non-trivial, as illustrated by antibodies marking filamentous structures. For example, microtubule-associated protein (MAP-2) is expressed in neuronal perikarya and dendrites, which in cerebral organoids have complex boundaries and lack stereotypical and layer-specific projection patterns [32,33,34]. Alternatively, mean fluorescence intensity is often used to measure overall antibody expression; however, this approach does not account for tissue size nor reveals the number of cells expressing specific markers [21,33,35,36]. This can lead to misinterpretations, as low cell density and high fluorescence signal might resemble low density and weak fluorescence signal when averaged over a large area. By enabling more precise cellular-level statistics, our approach to estimating cell numbers enhances research utilizing organoid models of human disease.

Overall, our objectives were to enhance whole organoid image analysis by providing a way to identify and separate areas of cell death in non-viable regions and quantify cell number and density of non-nucleic antibody markers in viable tissue areas. The resultant analysis workflow incorporated machine-learning-based cell nuclei segmentation implemented via Cellpose v2.2 (Howard Hughes Medical Institute, Ashburn, VA, USA) providing masked images that were further processed in MATLAB via a custom Graphical User Interface (GUI). The analysis workflow was developed to facilitate processing large datasets of tiled and z-series images via progression through a semi-automated step-sequence.

2. Materials and Methods

2.1. Ethics

Ethics for the project was approved by the Austin Health Human Research Ethics Committee (HREC/16/Austin/472, Melbourne, Australia) and conducted in adherence to the established guidelines set forth by the Austin Hospital, University of Melbourne, and the National Health and Medical Research Council (NHMRC) of Australia.

2.2. Human iPSC Model

Organoids were derived from two control cell lines, a hiPSC line derived from a neonate (Cat #WC026i-5807-3, WiCell, WI, USA) and an embryonic stem cell (ESC) line (Cat #H9 or WA-09, WiCell, WI, USA). Three batches of cortical organoids were generated from each cell line using the protocol established by Velasco et al. [12]. Each batch comprised 12 organoids at two developmental stages: six were collected on day 120 (4-months) post-induction and six on day 180 (6-months).

2.3. Immunohistochemistry

At 4 and 6-months post-induction, a minimum of two technical replicates were immersed in 4% paraformaldehyde and incubated for 2 hours at 4°C, then washed three times in 1X phosphate buffered saline (PBS), and placed in tissue-Tek O.C.T. Compound (Cat # 25608-930); then snap-frozen and stored at -80°C. Serial cryo-sections were cut at 20µm thickness (Cat # CM1950; Leica Biosystems, Buffalo Grove, IL, USA), collected on SuperfrostPlus glass slides (Cat # I6172PLUS, Thermo Fisher, Waltham, MA, USA) and stored at -80°C.

Immunohistochemistry was performed according to Yakoub and Sadek [37]. All primary and secondary antibodies are listed in Table 1. Sections were defrosted for 15 minutes at room temperature (RT), then blocked using 10% normal goat serum (NGS; Cat # 16210064; Invitrogen, Carlsbad, CA, USA) and 0.3% Triton X-100 in PBS, for a minimum of 2 hours, at RT. Primary antibodies were diluted in an antibody vehicle solution (AVS) of 1% NGS and 0.1% Triton X-100 in PBS) and sections incubated overnight at 4°C, in the dark. Table 2 details the number of sections included in our analysis. Slides were washed in a solution of 0.05% Tween20 and PBS and incubated with the secondary antibodies (diluted in AVS) for one hour, at RT (see Table 1). A second wash using a solution of 0.05% Tween20 and PBS preceded a 5-minute incubation with 6-diamidino-2-phenylindole (DAPI). Two final PBS washes were undertaken then sections were mounted using Invitrogen ProLong Gold Antifade mounting medium (Cat # P36934, Thermo Fisher, Waltham, MA, USA) and covered with Menzel-Glaser glass coverslips (Cat # 11911998, # 1.5H, Thermo Fisher, Waltham, MA, USA). Slides were subsequently stored at -80°C until imaged.

To characterize the morphology of the organoids with respect to protein markers indicative of mature expression patterning, three antibodies were used. In mammalian cells S-100 (β-subunit) (S100β) protein is predominantly expressed by mature astrocytes and distributed throughout the cytoplasm [38,39,40,41]. Microtubule-Associated Protein 2 (MAP2), a cytoskeletal protein mainly expressed in the soma and dendrites of mature neurons [32]. Gamma-aminobutyric Acid (GABA)ergic neurons, an inhibitory neurotransmitter chiefly located in the cytoplasm and dendrites of neurons [42,43]. A total of 88 sections were included in the final analysis, as outlined in Table 2. Additional organoid sections (n=6) were included to verify non-viable cellular and core features using cleaved-caspase-3 antibody and TUNEL assay kit (Table 1).

2.4. Imaging

Images of each organoid slice were acquired using a 20X/0.8NA objective on a Leica THUNDER Imager DMI8 fluorescence microscope (Leica Microsystems, Germany). Image z-stacks with a z-step size of 0.571µm and a pixel size of 0.323µm xy, were acquired to ensure adequate sampling in line with recommendations for post-processing using Leica Application Suite X (LASX, version 3.8.1.26810, Leica Microsystems, Germany). A tile scan function was used to rapidly image the entire organoid slice, with a 10% tile overlap. All sections were labelled with DAPI for nuclei segmentation (excitation wavelength used for imaging 391 nm) and either S100β (ex 479 nm), or MAP2 or GABA (ex 554 nm). The images were subject to small volume computational clearing (SVCC) using the LASX software [44], before being exported. Settings included feature scale - 6000 (nm) and 30% background scale. In addition, confocal images were also acquired on a subset of organoid slices using an LSM900 scanning confocal microscope 20X/0.8NA (Carl Zeiss Inc., Germany). Confocal imaging time was considerably longer at approximately 6 hours (parameters optimized for Nyquist sampling), versus 10 minutes on Leica THUNDER (parameters optimized for LASX processing as noted above). The Leica THUNDER system therefore presented a significant advantage to achieve high throughput and was used for analysis workflow development due to capturing sufficient image quality to segment and analyze antibodies for cell counts. All images were acquired within the Florey Neuroscience Microscopy Facility (The Florey Institute of Neuroscience and Mental Health, Melbourne, Australia).

2.5. Image Data Analysis

A custom set of MATLAB (Version R2022b, Mathworks, Natick, MA, USA) functions were developed for the purpose of image analysis. These incorporated a range of steps, including tile stitching, nuclei and peri-nucleic antibody segmentation, cell-state characterization, the identification of the non-viable cell core boundary and applicable calculations. A description of the image analysis workflow is incorporated in the results section outlined below. All analysis scripts are available on request. The analysis workflow was executed on a Dell PowerEdge R820 Server, equipped with four Intel Xeon E5-4620 processors (2.20GHz, 16M Cache, 7.2GT/s QPI, Turbo, 8 Core) and supported by a RAM capacity of 512GB.

2.6. Statistical Framework

All statistical analyses and figures relating to the visualization of data were completed using GraphPad PRISM software (Version 9.5.1, GraphPad Software Inc., La Jolla, CA, USA). Numeric values are reported to three significant figures. and data is presented as mean ± Standard Error of the Mean (SEM). To compare groups, multiple t-tests with Welch’s correction were conducted to account for unequal variances. To compensate for the potential inflation of Type I errors, the False Discovery Rate (FDR) – using the Benjamini-Hochberg method (Q = 5%) – was also determined. Comparisons with an FDR-adjusted p-value (q) less than 0.05 were considered discoveries.

3. Results

Cortical organoids obtained from two stem cell control lines were collected at two time points, after 4 and 6-months of differentiation. A workflow for their analysis is presented in Figure 1 and included cortical organoid sectioning (Figure 1A-C), immunohistochemistry and image acquisition (Figure 1D-G). To quantify the cellular composition of organoids, we developed a customized image analysis pipeline.

3.1. 3D Analysis of Organoid Slices

As previously outlined, fluorescently labelled cortical organoid slice images were captured using Leica Application Suite (LASX). To preserve the inherent 3D tissue structure, image processing of z-stacks, which comprised 2D tiled images taken across the x-y plane, were reconstructed as a 3D spatial model of the organoid sections. Automated analysis of the entire organoid section was performed to capture the intrinsic structural heterogeneity and provide a more representative sampling approach to quantifying organoid morphological features such as cell density, viability, and antibody expression across developmental stages.

3.2. Identification of Tissue Boundaries

Automated and accurate detection of the outer boundary of each organoid ensured a clear demarcation of the tissue (Figure 2A-B). The identification of tissue artefacts such as tears, folds, non-viable core and their exclusion from further analysis was implemented within the image processing workflow (Figure 2C-D). By automating the detection of sectioning and other artefacts, this process helped filter out non-biological features, ensuring more reliable measurements. Figure 2D is a schematic of the viable region after refinement of the section boundaries and illustrates only the viable cells as described in the next section.

3.3. Nuclei Segmentation

The analysis workflow enabled accurate segmentation of nuclei (Figure 3) and provided the organoid cell counts. Nuclei segmentation was based on DAPI-labelling. Algorithm selection was a crucial step, as tightly packed cellular structures within cortical organoids can be highly confounding for cell counts. Several tools for 3D object segmentation were explored, including Imaris (Oxford Instruments, Switzerland), Cellprofiler (Broad Institute, Cambridge, MA, USA), Cellpose v2.2 and Image J v1.54f. [45]; example segmentation results shown in Figure 3A-D. Each package was evaluated for speed, reliability, and ability to adapt to a heterogenous sample. Compared to other packages tested, the machine learning framework utilized by Cellpose v2.2 demonstrated proficiency when faced with complex and varied samples and was embedded into the analysis workflow. 3D segmentation was completed using the nuclear model with a minimum mask size of 15 pixels and an initial object diameter of 30.8 pixels. Anisotopy was defined as the ratio of z-step size (µm) to pixel size (µm) for each image. The resulting segmented DAPI-labelled images, outputted as 16-bit Tag Image File Format (.tif) files, provided a basis for subsequent analyses within the workflow.

3.4. Cell-State Classification

Nuclei were then classified into four distinct groups: cell fragments, non-viable, viable, and clustered or dividing cells (Figure 3E). This classification was initially based on nuclei volume, calculated from the segmented images outputted by Cellpose v2.2. The volumes were fitted to a Gaussian Mixture Model (GMM) with three components (Figure 3E). The GMM distinguished between very small volumes ($\overline{x}$ = 8.57µm³), representing cell fragments or improperly segmented nuclei; moderate volumes, representing single cells; and large volumes ($\overline{x}$ = 548µm³), corresponding to dividing or clustered cells. Fragments were excluded from further analysis (Figure 3F).

DAPI stained nuclei of apoptotic or dead cells typically exhibit lower volumes and higher fluorescent intensities due to condensed DNA and cytoplasm compared to their viable counterparts [46,47]. These features were leveraged for the application of a second GMM, used to classify the group of moderate volumes as either viable ($\overline{x}$ = 219µm³) or non-viable ($\overline{x}$ = 149µm³; apoptotic/dead) cells. The estimated mean volume of nuclei within the non-viable and viable regions were provided as initial parameters for the GMM, which incorporated both volume and intensity data. For the purposes of clarity, dying and dead cells are described as non-viable for the remainder of the report. The sequential application of GMMs resulted in each nucleus being allocated to one of the four final groups (Figure 3F).

The results obtained for the nuclei feature detection and cell classification were confirmed using a separate set of organoid sections containing a non-viable cell core. To identify the apoptotic pathway, we included cleaved-Caspase-3 antibody and Terminal deoxynucleotidyl transferase dUTP nick-end labelling (TUNEL) assay (n=3 each). Classification of non-viable cells, at high density within the organoid core, was confirmed by both cleaved-Caspase-3 antibody (Supp. Figure 1A-C), and TUNEL assay (Supp. Figure 1D-F), which showed positive staining concentrated within the organoid core. To verify, at a cellular level, non-viable versus viable cells, the mean fluorescence intensity and diameter of nuclei for either DAPI/TUNEL positive (non-viable) and DAPI⁺ve/TUNEL^-ve (viable) was calculated using Image J (n=20 each). Nuclei intensity plots were obtained using a line-segment 1 µm-width (Supp Figure 1H, I). Classification of non-viable cells using DAPI-nuclei staining confirmed that small diameter (< 6mm), high signal intensity nuclei were TUNEL positive; while lower signal intensity DAPI⁺ve cell nuclei (>6mm) were TUNEL negative, consistent with viable cell classification used in the analysis workflow.

3.5. Identification of Non-Viable and Viable Regions

To obtain organoid measurements including for example, the core size and cell density, the analysis workflow was developed to automatically delineate viable and non-viable areas.

3.5.1. Non-Viable Regions

Automatically demarcating the non-viable cell core is challenging due to the intermingling of dead and viable cells throughout the organoid. This phenomenon is underpinned by several factors, including the gradient of nutrient diffusion, natural rates of cell death, incurred damage, and variations in hypoxia resistance [48,49]. To demarcate the presence of a core and approximate core boundary, a novel protocol was applied. Binary masks were generated of the non-viable and viable nuclei in each segmented image tile, and Maximum Intensity Projections (MIPs) were created for each of the two masks. Local cell densities were then calculated for both the non-viable and viable MIPs using a sliding window convolution. The necrotic density map was then subtracted from the viable density map, and regions with a high ratio of non-viable cells were identified using a threshold-based method, marking the necrotic core. Finally, morphological smoothing and hole-filling operations were applied to refine the necrotic core boundary.

3.5.2. Viable Region

The viable regions in all samples were calculated from the total area of the section, delineated by the outer boundary, minus the region containing the non-viable cell core (example Figure 2C).

3.5.3. Non-Viable Cells Within the Viable Region

Nuclei centroids positioned within the viable region of the samples were extracted for analysis. The relative frequency of cell death within the viable region was subsequently calculated as the percentage of non-viable cells relative to the total number of cells within the viable region.

3.6. Estimation of Antibody Positive Cell Count

This section outlines the portion of the analysis workflow used to quantify the number of cells expressing S100β, MAP2, and GABA antibodies. As a brief overview, first the antibody-positive voxels were identified by thresholding the antibody-labelled images. Next, a peri-nuclear shell was obtained by dilating the DAPI-mask corresponding to each viable cell nuclei. Finally, the percentage of antibody-positive voxels within the shell was measured. A schematic to distinguish the spatial classification of nuclei and peri-nucleic markers is given in Figure 4A-D for this part of the analysis workflow and described in detail below.

3.6.1. Antibody Positive Voxels.

The measurement of antibody expression involved the application of a cutoff to the antibody image. An adaptive threshold was selected according to Otsu’s three-class method, a classification method based on intensity, with the middle class assigned to the background. This is an established approach that was selected for its robustness when applied to antibody-stained images with multi-modal characteristics [50,51]. The result was a binary image of all voxels identified as positive for the antibody, from which the percentage volume occupied by the antibody was calculated.

3.6.2. Delineation of Peri-Nucleic Region

The masked images of the segmented nuclei were dilated using a spherical structuring element three voxels (~1µm) in diameter (Figure 4A-B). This resulted in a mask of the three-dimensional “shell” around the nuclei, or the peri-nucleic region (Figure 4C).

3.6.3. Identification of Antibody-Positive Cells

In order to distinguish antibody, “positive” cells from antibody “negative” and to account for cells whose nuclei were positioned in close proximity to the antibody labelling expressed by other cells; a cutoff method was applied. The cutoff value was guided by two experts in the domain, who visually classified cells from a single 3D tile for each antibody as either positive or negative (Figure 4D). In Figure 4E are example images of the antibody labelling highlighting cells which were manually identified as positive. Positive cells were defined as those with a 3D fluorescence signal surrounding a DAPI-labelled nucleus. The mean percentage of antibody signal within the peri-nucleic shell of visually identified positive cells was calculated and used to set the minimum cut-off value for classifying antibody positive cells in the image processing workflow. The thresholds were established as 31.5% for GABA, 23.2% for MAP2, and 10.4% for S100β, with any cell exceeding the respective cut-off considered positive. The percentage of antibody-positive cells within the viable regions were then calculated.

3.7. Volume Estimation

To estimate the volume of the slice we combined the x-y pixel map with the z-coordinates of the segmented nuclei. This was performed as a heatmap spatially averaging this information across sub-regions, 200x200 pixels in size. Overlaying the pre-defined inner boundaries of the non-viable core and outer boundary of the organoid slice their volumes could be estimated.

3.8. Temporal Analysis of Organoid Morphology and Maturation

The analysis workflow was applied to measure the developmental trajectory of control cortical organoids between 4 and 6-months post-induction. Measurements included the volume of the non-viable core, cell densities, and three cell populations, comprised of astrocytes (S100β) and neurons, which included all neurons (MAP2) and GABAergic interneuron sub-type. In total, 88 sections from 31 organoids were analysed, as outlined in Table 2.

3.9. A Proportional Decrease in Non-Viable Cells Was Observed from 4 to 6-Months Post-Induction

The structural health of the organoids was assessed by the relative size of the non-viable cell core. A non-viable cell core was identified in 40 (86.9%) of those sampled at 4-months post-induction and 32 (61.5%) of those sampled at 6-months post-induction. There was no change in the total volume of the organoid sections examined over time (Figure 5A). However, analysis of the volume of the non-viable core as a percentage of the total volume was significantly reduced with a mean of 26.4% versus 18.8% at 4 month and 6 months post-induction respectively (Figure 5B); Difference = 7.6%, SE of difference = 3.1%, t(62.9) = 2.5, p = 0.0168, q = 0.0176.

The cellular health of the organoids was evaluated by calculating the percentage of non-viable cells relative to the total cell count within the viable region of the samples. No discovery was made when comparing mean ratios at 4 months ($\overline{x}$ = 22.0%) and 6 months ($\overline{x}$ = 21.8%) post-induction (Figure 5C); Difference = 0.19%, SE of difference = 3.1, t(66.26) = 0.096, p = 0.924, q = 0.485.

3.10. Cell Density in Cortical Organoids Was Higher at 4-Months Versus 6-Months Post-Induction.

Cell density was utilized to measure cell proliferation and evaluated in the viable region of cortical organoids. The viable cell density in this region was higher at 4 months ($\overline{x}$ = 31.3 cells/100µm³) compared to 6 months ($\overline{x}$ = 22.8 cells/100µm³) post-induction (Figure 5D); Difference = 8.4, SE of difference = 0.794, t(64.4) = 10.6, p < 0.000001, q = < 0.000001.

3.11. Quantifying Expression of Astrocyte and Neuronal Markers.

3.11.1. S100β Expression Was Unchanged from 4 to 6-Months Post-Induction.

The maturation of individual astrocytes was assessed using the ratio of S100β-positive cells to all viable cells. At 4 months the mean percentage of S100β-positive cells was 11.1%, compared to 13.6% at 6 months post-induction (Figure 6A). No discovery was made with this comparison; Difference = -2.5%, SE of difference = 1.7, t(25.8) = 1.4, p = 0.162, q = 0.346. The mean density of S100β-positive cells was also similar at 4 ($\overline{x}$ = 2.9 cells/100µm³) and 6 months ($\overline{x}$ = 2.5 cells/100µm³) post-induction (Figure 6B); Difference = 0.376, SE of difference = 0.404, t(31.7) = 0.93, p = 0.36, q = 0.378.

3.11.2. No Change in MAP2 Expression from 4 to 6-Months Post-Induction.

Temporal assessment of neuron maturation was made using the percentage of MAP2-positive cells in organoids at 4 and 6 months. As a proportion of all viable cells, MAP2-positive cells showed a decrease from 8.7% to 7.4% from 4 to 6 months post-induction respectively, which did not result in a discovery (Figure 6C); Difference = 1.3, SE of difference = 2.3, t(14.8) = 0.577, p = 0.573, q = 0.502. Similarly, no discovery was made in relation to the differences in density of MAP2-positive cells from 4 ($\overline{x}$ = 1.95 cells/100µm³) to 6 months ($\overline{x}$ = 1.3 cells/100µm³) post-induction (Figure 6D); Difference = 0.602, SE of difference = 0.445, t(13.8) = 1.35, p = 0.198, q = 0.346.

3.11.3. The Proportion of GABA-Positive Cells Increased, But No Change in GABA-Positive Cell Density 6-Months Post-Induction.

The differentiation of GABAergic neurons was assessed by the percentage of GABA-positive cells at 4 ($\overline{x}$ = 8.6%) and 6 months ($\overline{x}$ = 14.5%) post-induction. A discovery revealed a significant increase in the proportion of GABA-positive cells over time (Figure 6E); Difference = -5.9%, SE of difference = 1.2, t(36.8) = 4.9, p = 0.00002, q = 0.0001. An increase in the density of GABA-positive cells was observed between 4 ($\overline{x}$ = 2.07 cells/100µm³) and 6 months ($\overline{x}$ = 2.4 cells/100µm³) post-induction, though no discovery was made regarding this change (Figure 6F); Difference = -0.286, SE of difference = 0.258, t(26.0) = 1.1 p = 0.278, q = 0.365.

4. Discussion

We have established an image analysis workflow to serve as a standardized method for the assessment of organoid viability and cytoarchitecture. The workflow was tested on a dataset of more than 1.5 million cells and demonstrated remarkable efficacy in identifying specified cell types and states, as well as providing temporal and spatial information when applied to cortical organoid development.

To enhance our analysis beyond single-plane imaging, we captured the entire organoid slice in 3D, revealing details of the cellular organization that traditional maximum intensity projection methods omit. Additionally, our automated analysis of the complete section enabled us to explore the regional diversity within the organoids across each group. This comprehensive approach provided deeper insights into cellular features such as density, viability, and antibody expression at various developmental stages.

To implement high-throughput imaging we utilized a widefield fluorescence microscope system. The advantage of this imaging approach provided 3D acquisition of entire sections in a fraction of the time needed for confocal microscopy when using similar acquisition parameters. The widefield system adequately captured the nuclear marker and the peri-nuclear environment, which were the primary focus of our analyses. These findings suggest that widefield microscopy is a suitable alternative when a balance between resolution and higher acquisition speed is required.

The Leica Thunder system has three computational clearing algorithms available. Here we implemented Small Volume Computational Clearing (SVCC). This algorithm applies adaptive deconvolution before removing unwanted background [44]. Image datasets captured using this method provided sufficient cellular detail to facilitate automated cell counts. The resultant analysis workflow was designed to process widefield 3D image tiles; however, it could also be applied to any fluorescence image datasets for example, optically-sectioned whole tissue samples using light-sheet microscopy [21,22].

At present, research using organoid models lacks a well-defined benchmark for 'normal' growth and development. We have established an “in-house” method to benchmark the elected organoid protocol tailored to the specific needs of cytoarchitectural analysis, an exercise critical to identify, interpret and compare deviations between control and disease models.

The analysis workflow’s initial task was to quantify non-viable and viable cells using the nuclear label. Non-viable or necrotic cells have previously been identified in organoid image-based studies by their higher fluorescence intensity and smaller diameter nuclei in comparison to functional cells [52,53,54]. While previous studies focused on a single attribute, we adopted a more comprehensive approach by implementing a two-tiered Gaussian Mixture Model that analyzes both volume and intensity. This approach allowed for more complex stratification of candidate cells. We show that nuclei classified using the integrated approach as non-viable cells was confirmed using cleaved-caspase 3 and TUNEL staining, consistent with previous reports [12,36,47,55]. The approach also allowed automatic delineation of the non-viable core boundary from the viable region, which was a significant advancement.

To evaluate organoid viability, the nucleic viability classifications were applied to obtain quantitative “in-house” assessment of cortical organoids including the non-viable core volume and percentage non-viable cells in the viable region. We found a significant decrease in the volume of the non-viable core between 4 and 6-months post-induction. In contrast to this finding, previous research has proposed the “dead cell” core increases over time, caused by the limited diffusion of oxygen and nutrients from the organoid surface in the avascular model [14,48,56]. To interrogate the current finding, we quantified two additional factors likely to affect core size including the densities of viable and dead cells within the viable region. Our results indicated a reduced viable cell density within the viable region, from 4 to 6-months. Second, we found the viable region showed a decreased number of non-viable or dead cells from 4 to 6-months, which is consistent with reported cell death rates in cerebral organoids plateauing as they near synaptic maturity [57]. Taken together these results indicate diffusion rates across the viable region may be higher at 6 months. Additional transcriptional analysis would help to determine whether changes in gene expression involved in aerobic respiration and apoptotic genes are consistent with current findings [58]. Overall, the results suggest that evaluating the non-viable core and viable region’s growth, at a cellular level can help to understand how these factors might influence the health and development of cortical organoids in control and disease models.

To quantify the number of protein labelled cells, the analysis workflow was designed to capture each cell’s peri-nuclear protein expression level. Micro-structure in mature cortical organoids is complex in terms of both tightly packed heterogeneous cells and sub-cellular protein expression patterns, which make automating and obtaining cell counts for non-nucleic markers challenging. Therefore, our approach was to use each viable cell’s DAPI-mask and create a dilated 3D volume or peri-nuclear shell.

With respect to quantifying organoid maturation, our assessment between 4 and 6-months post-induction showed no significant differences in S100β expression, indicating that astrocytic maturation remains stable throughout this period. This finding was somewhat unexpected, given the literature describing a continual increase in proliferation from the beginning of astrogenesis in cerebral models [12,20,59,60]. Two possible explanations for these findings arise: First, it may be that the rate of astrocytic maturation is variable over time and that this variation has not previously been captured. Second, although widely applied as a measure of astrocytic maturation, S100β is not expressed in all mature astrocytes [41,61]. Therefore, quantitative changes in maturation based on S100β may be limited, and additional markers like glial fibrillary acidic protein (GFAP) could provide further insight.

The percentage of all cells characterised as mature MAP2 positive neurons were similar from 4 to 6-months post-induction; while the percentage of GABAergic neurons was significantly increased. These findings suggest a potential shift in the neuronal population towards GABAergic neurons, consistent with previous studies that revealed a shift occurs at approximately 3-6 months post-induction [9,12]. These observations emphasize how the image analysis workflow can be used to evaluate cell number in models of neurodevelopmental disorders, particularly where GABAergic interneurons are thought to have a significant impact [9,62,63].

Future application of the image analysis workflow could include machine learning development. By capturing spatial and single-cell-level information from sample images, including protein signatures, the workflow can also be expanded to incorporate the cells’ molecular profile. This spatial mapping enables detailed profiling of each cell across large replicate datasets. Such advancements align with proposed developments in machine learning, where immunohistochemistry analysis of patient-derived organoids could serve as a valuable tool for personalized drug screening platforms [64,65].

In conclusion, we established a method to benchmark the preferred organoid protocol “in-house”. Functionalities were built into a single, customizable script and incorporated changes to improve their theoretical robustness, tailored to the specific needs of cytoarchitectural cortical organoid analysis. Implementing the current workflow would help laboratories to establish their own benchmark; a critical exercise for identifying, interpreting and comparing deviations across disease and control models.