Deep-learning analysis of micropattern-based organoids enables high-throughput drug screening of Huntington's disease models.

Organoids are carrying the promise of modeling complex disease phenotypes and serving as a powerful basis for unbiased drug screens, potentially offering a more efficient drug-discovery route. However, unsolved technical bottlenecks of reproducibility and scalability have prevented the use of current organoids for high-throughput screening. Here, we present a method that overcomes these limitations by using deep-learning-driven analysis for phenotypic drug screens based on highly standardized micropattern-based neural organoids. This allows us to distinguish between disease and wild-type phenotypes in complex tissues with extremely high accuracy as well as quantify two predictors of drug success: efficacy and adverse effects. We applied our approach to Huntington's disease (HD) and discovered that bromodomain inhibitors revert complex phenotypes induced by the HD mutation. This work demonstrates the power of combining machine learning with phenotypic drug screening and its successful application to reveal a potentially new druggable target for HD.

Introduction

The neurological space is in desperate need for new drugs, but clinical successes have been extremely low. For rare neurological disorders, only 5% of all drugs entering clinical trials are estimated to ultimately reach Food and Drug Administration approval in the United States (Wong et al., 2019). The failure to efficiently provide first-in-class medicines has been linked to the target-based methodology that has fueled most drug-discovery programs for decades (Swinney and Anthony, 2011). This approach starts from a defined molecular target that is hypothesized to be linked to the etiology of the disease. A radically different methodology, the phenotypic approach, proposes to use cell-based phenotypic readouts and does not rely on prior knowledge of the disease pathogenesis. Thus, it allows for an unbiased identification of new mechanisms of action and has great potential for delivering new first-in-class drugs, especially for poorly understood neurological disorders (Moffat et al., 2017; Swinney, 2013; Zheng et al., 2013).

With the quick development of organoid-based technologies, we are now able to generate powerful pre-clinical models of human diseases for phenotypic high-throughput screening (HTS). While still lacking key physiological features such as an immune system, blood-brain barrier, or aging hallmarks, these models incorporate the complex interplay of signaling and morphogenesis in a multi-fate, multi-tissue environment in a human background and thus reproduce the molecular and cellular defects induced by mutant Huntingtin (HTT) linked to Huntington’s disease (HD) more accurately. They therefore generate better models of human diseases, with the potential of discovering new drugs and new mechanisms of action for diseases with large unmet needs (Akkerman and Defize, 2017; Dutta et al., 2017; Lancaster and Huch, 2019; Ranga et al., 2014; Sirenko et al., 2019). However, while neural organoids have been used to mimic aspects of a number of diseases, typical three-dimensional (3D) organoids display low standardization and scalability, preventing their integration into HTS campaigns. Overcoming these limitations is an area of intense research. One convenient solution is to rely on neural spheroids (Sirenko et al., 2019), which are more standardized but do not recapitulate the radial organization of the cortex. Another solution, technically demanding, is to introduce automatization and tissue-clearing methodologies (Renner et al., 2020). It is thus imperative to keep devising new methods that allow for a large number of organoids to be generated simultaneously with high reproducibility as well as to develop new computational tools for the detection of potentially subtle disease phenotypes that have to be extracted despite the inherent variability of multi-cellular biological processes.

Here, we report the design of a generic phenotypic screening methodology at the organoid level and the discovery of new modulators of phenotypes of HD, a fatal dominant autosomal neurodegenerative disease that is caused by an increase in the number of CAG repeats, which expands a polyglutamine (polyQ) tract in the HTT protein (Ross and Tabrizi, 2011). While being a degenerative disease, HD has recently been shown to alter human neurodevelopment in human fetal samples (Barnat et al., 2020), suggesting that early human organoid models that reproduce aberrant signaling and morphogenesis in HD constitute a promising approach for discovering new mechanisms that are also relevant at later stages of the disease. Here, we leveraged the recently devised “neuruloids” that use micropattern-based differentiation to create organoids mimicking the ectodermal compartment during human neurulation (Haremaki et al., 2019), with the possibility for easy upscaling while retaining excellent reproducibility. We then performed a drug-discovery screen aimed at reversing a complex phenotype that we have previously reported for HD in the neuruloids, using a previously characterized isogenic series of human embryonic stem cell (hESC) lines with graded increases in CAG lengths (Ruzo et al., 2018). Finally, we developed a deep-learning computational pipeline to analyze the screening results and to quantify, for each compound, its efficacy at reversing the disease phenotype back to normal as well as its adverse effects. Using this combination of tools, we discovered that specific bromodomain inhibitors can efficiently revert HD phenotypes to wild type (WT) and can also alleviate neuronal susceptibility to apoptosis in human HD neurons in vitro, highlighting a potential new druggable target for HD that should be further evaluated.

Results

Development and validation of a micropattern-based organoid screening platform

We and others have recently reported the generation of organoids that reproduce the entire ectodermal compartment of the human embryo during its fourth week of development (Britton et al., 2019; Haremaki et al., 2019). These self-organized structures, called neuruloids, reproducibly and accurately generate all four ectodermal lineages (neural, neural crest, placode, and epidermis) in the correct layered 3D organization as observed in vivo, based on a minimal protocol. Moreover, using an isogenic set of cell lines carrying different mutations associated with HD, we have demonstrated that HD neuruloids exhibit a complex, stereotypical phenotype. We hypothesized that these findings could be scaled up for usage as a drug-screening application which would reveal drugs that reverse the disease state back to the normal state (Figure 1A). Our motivation is to use these early developmental organoids to discover molecules that can rescue the deleterious impact of the HD mutation in the context of a complex neural tissue. Successful hit molecules will then be tested in alternative HD-specific secondary assays, with the aim of discovering molecules rescuing a breadth of the pleiotropic deleterious functions, as yet mostly unknown, that are associated with mutant HTT. Such molecules will rescue deleterious functions of mutant HTT and might eventually lead to the development of drugs for preventing disease onset.

Figure 1

High-throughput screening strategy and validation

(A) Schematic illustration of the high-throughput screening (HTS) strategy. Organoids carrying the mutation for Huntington’s disease (HD) are treated with drug compounds and then analyzed to establish whether they revert the phenotype to the wild type (WT).

(B) A deep neural network can be used as an efficient classifier to distinguish between the phenotypes and to determine that degree of reversal.

(C) Example wells of a 96-well plate for WT and HD organoids (diameter 700 μm) stained for the early neural progenitor marker PAX6 and the actin marker phalloidin.

(D) Example WT and HD images and averaged images over all control organoids (diameter 700 μm).

(E) Classification accuracy of the neural network trained on control WT and HD images. The accuracy is determined by validating on a set of images that were not used during training (individual images: n = 1,464, nWT = 700, nHD = 764; averaged per well: n = 58, nWT = 29, nHD = 29).

(F) Comparison with other discrimination methods in terms of distance between WT and disease and the resulting Z′ factors show superior accuracy of the neural network approach compared with image clustering using UMAP or feature segmentation using random forest classification (each data point represents average per well, all n = 29).

(G and H) Validation of the neural network approach for identifying phenotypic reversal using an HTT knockout and overexpression (OE) (n = 54, nWT = 10, nKO = 11, nKO+OE = 33). The neural network correctly identifies the rescued HTT−/− + HTT-OE as WT.

Precise separation between negative and positive controls is a requirement for HTS platforms in order to minimize the rate of false positives and false negatives. This is commonly measured by the Z′ factor, a quantity smaller than 1 that considers the separation between controls as well as the relative spread of the data (Zhang et al., 1999). The use of organoids in HTS is limited by inherent noise in multi-cellular processes, but first and foremost by the current high inter-organoid structural variability. We hypothesized that the combination of the neuruloid, a highly reproducible micropattern-based organoid, and the image classification power of deep neural networks (Figure 1B) would allow a strict classification of phenotypes despite the noise resulting from underlying biological processes related to multi-cellular system morphogenesis, and would thus allow us to meet the stringent statistical requirements associated with HTS. In our case, this refers to a strict separation between WT and HD neuruloids, and thus a high Z′ score.

To scale up organoid production, neuruloids were grown on 700-μm-diameter, disk-shaped micropatterns at the bottom of 96-well plates, resulting in approximately 35 neuruloids per well, 25 of which were not cut by the edge of the well and were thus suitable for downstream analysis. While the original published protocol used 500 μm diameter, we chose commercially available 700-μm-diameter micropatterns to perform the primary screen before validating hits in the originally reported 500-μm-diameter setting. For the larger 700-μm-diameter patterns, the neuruloids showed a slightly less compact neuroepithelial core compared with our previous results with 500-μm patterns; however, this did not affect reproducibility, the HD phenotype, or our ability to differentiate between WT and HD neuruloids. Immunofluorescence for the neural marker PAX6 and the actin marker phalloidin served as a readout for both fate patterning and structural organization (Figure 1C). Consistent with our previous report, we observed a clear phenotype of the HD neuruloids compared with WT with high reproducibility, in particular an expansion and change of shape of the PAX6+ neural core, highlighted in individual neuruloids as well as in the average of channel intensity between all neuruloids of a similar well (Figures 1D and S1). We next asked whether we could computationally separate the two phenotypes with a low enough false detection rate for usage in a drug screen. Leveraging the power of deep convolutional neural networks (CNNs) for image classification tasks, in particular their success at extracting and combining an optimal set of image features, we trained a CNN for recognizing immunofluorescent images of WT and HD neuruloids (Figures 1B and S1; STAR Methods). This approach led to a near-perfect separation of WT and HD phenotypes at the level of individual neuruloids as well as for the average score per well (Figure 1E). To quantitatively evaluate the performance of the CNN-based discrimination approach in the context of HTS, we calculated the corresponding Z′ factor and compared it with other machine-learning techniques (Figure 1F). A general consensus is that a high-quality screen has 0.5 < Z′ < 1, while Z′ < 0 is unacceptable (Zhang et al., 1999). Strikingly, our CNN approach provided an excellent distance between positive and negative control, resulting in a Z′ factor close to the optimal value of 1, thus representing a state-of-the-art screening platform. Conversely, alternative classification approaches (Figure 1F), e.g., uniform manifold approximation and projection (UMAP; Becht et al., 2019) clustering of the raw images or feature segmentation using Ilastik (Berg et al., 2019), a random forest classifier, resulted in a Z′ factor below 0, rendering the approach unusable for screening.

Ultimately, this tool should not only recognize different phenotypes but also quantify biological rescue. We therefore took advantage of the previously reported strong phenotype observed in neuruloids with a knockout (KO) of the Huntingtin gene, HTT−/− (Haremaki et al., 2019), to design a proof-of-concept rescue experiment where WT Huntingtin protein (HTT) was overexpressed in the HTT−/− background. As previously reported, the HTT−/− neuruloids showed an enlarged PAX6 center compared with WT (Figure 1G). Moreover, the PAX6+ domain of HTT−/− neuruloids was not organized around a closed central cavity, as highlighted by N-cadherin and phalloidin stains. Constitutive expression of WT HTT in the HTT−/− background led to a clear rescue of these defects (Figure 1G). We first trained our classifier to distinguish between WT and HTT−/− neuruloids. As expected, when new images were presented to the network, WT neuruloids had a probability of belonging to the WT class (prob. WT) close to 1, while HTT−/− neuruloids had prob. WT close to 0 (Figure 1H). This demonstrates successful training of the network. When we used the images of the HTT−/− neuruloids with constitutive expression of WT HTT, the probability of belonging to the WT class reverted back to 1, demonstrating a rescue of the HTT−/− phenotype toward the WT (Figure 1H). In conclusion, we have established a tool that not only can distinguish WT and mutant phenotypes with the highest accuracy but also can be used to quantify phenotypic rescue.

Determining adverse effects of compounds

We next asked whether we could determine the potential adverse effects of compounds on the neuruloids (Figure 2A). These adverse effects would globally encompass potential compound toxicity, but also off-target effects which would affect neuruloid disease phenotypes that are not in the direction of the WT “rescued” domain. Obtaining an efficient measure of adverse effects simultaneously with the degree of rescue would be highly advantageous in ranking therapeutic candidates. However, there are many different ways in which compounds could adversely affect a neuruloid, not all of which are known a priori. Therefore, it is not possible to train a classifier for this task. In addition, although neural network classifiers perform well for in-class data, they are known to be inaccurate for data belonging to classes to which it has not been exposed in the training phase. In order not to bias our analysis to pre-defined classes, we thus chose to measure adverse effects using unsupervised learning based on convolutional autoencoders (Figure 2B). These neural networks extract the most important features of the data by compressing images unbiasedly to a vector of defined length, the latent vector (Figure 2B). The latent vector can then be used as a representation of a phenotypic space, in which distances from control neuruloids can be used to determine the adverse effects of a compound (Figure 2C): the distance between WT and HD neuruloids is interpreted as the disease direction along which the therapeutic rescue will happen, while the orthogonal distance is interpreted as a measure of adverse effects. To test the accuracy of this approach, we performed a pilot experiment applying compounds with known toxic effects to the HD organoids at different concentrations and quantitatively compared this value with 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT), a conventional colorimetric assay for cytotoxicity. As expected, we found control WT and HD neuruloids to be well separated into two different clusters in the phenotypic space (Figure 2D, left). After computationally removing the disease direction defined by the center of the WT and HD clusters, we were able to measure the adverse effects of compounds as the distance from the central cluster in the space defined by the remaining dimensions (Figure 2D, center and right). Since this distance is a distance to a distribution of control neuruloids, we used the Mahalanobis distance (Mahalanobis, 1936), which takes into account that the reference control distribution may not be uniformly distributed in the latent space, and that therefore distances should be rescaled by the variance of the control distribution in each dimension. We found that the dimensionality of the latent space did not have a significant effect and used a fixed latent space of dimension of 64 (Figures S2A–S2D). To compare the results from the autoencoder quantitatively with the MTT assay, we defined a distance of more than three times of the average distance of control neuruloids as toxic. Distances were then scaled to lie between 0 and 1 and, for each well, compared with the light absorption of the colored MTT solution. A comparison between this method and the MTT assay for toxicity showed very good agreement (Figure 2E; example images shown in Figure 2F). We therefore conclude that our method can precisely measure phenotypic rescue as well as the adverse effects of a given compound.

Figure 2

Strategy for determining adverse effect and phenotypic space

(A) Schematic illustration of the method for quantifying the adverse effect of a compound.

(B) Overview of the autoencoder architecture used for encoding organoid data into a latent vector.

(C) Illustration of the phenotypic space spanned by the latent vectors of the autoencoder.

(D) adverse effect estimation procedure and phenotypic space of a screen with 1,065 compounds (each data point is average of one well; n = 1,257, nWT = 96, nHD = 96, nHD+compound = 1,065).

(E) Comparison of adverse effect scores for different drugs and varying concentration obtained by the method described here and the conventional MTT assay shows high accuracy of our method (high scores correspond to low adverse effect, error bars indicate standard deviation over triplicates).

(F) Example phalloidin stains for the drugs and concentrations shown in (E). Panels missing are associated with high toxic compound concentrations that resulted in cell death and no structure to be imaged. Micropattern diameter, 700 μm.

Small-molecule screen identifies compounds that revert neuruloid phenotype of Huntington’s disease

Our methodology allowed us to perform a forward screen with the aim of discovering compounds that can revert the HD neuruloid phenotype. To obtain unbiased coverage of a large number of different mechanisms of action, we used a structurally diverse set of 1,065 bioactive compounds at 10 μM. The results of the screen are presented in the rescue-adverse effect space in Figure 3A. The phenotypic space, which was used to estimate the adverse effect of compounds, is shown in Figure 3B. The phenotypic space for untreated control wells of all 14 compound plates demonstrating excellent reproducibility across plates is shown in Figure 3C, and rescue and adverse effect resolved per screening plate is presented in Figures 3D and 3E. Quality control for the screen revealed extremely low variability in the classification of the control plates (Figure S2E) and untreated control wells in each assay plate (Figure S2F).We defined hits as more than 95% rescue and less than three times the mean adverse effect distance of untreated control wells. Our data show that five hit compounds could revert the phenotype back to WT without presenting adverse effects to the tissue (see STAR Methods for hit description). Example images for the action of each of these compounds are shown in Figure 3F.

Figure 3

Screen for compounds that rescue the HD phenotype

(A) Phenotypic rescue and adverse effect scores for screen on the bioactive compounds library. Each data point represents the average score per well. Hit compounds are identified as high rescue (>0.95) and low adverse effect (<3) (n = 1,481, nWT = 96, nHD = 96, nHD+compound = 1,065, nadditional HD control = 224).

(B) Phenotypic space used for estimating the adverse effect for compounds.

(C) Control wells from 96-well compound plates show no substructure in the latent space, confirming high reproducibility across plates (each data point is average of one well; n = 224, nper plate = 16, nplates = 14).

(D and E) Phenotypic rescue (D) and adverse effect (E) for the 14 individual 96-well plates.

(F) Example images and details for hit compounds.

(G) Dose response for rescue and adverse effect for the hit compound bromosporine.

(H) Example images for the quantification of the dose response in (G).

All micropattern diameters, 700 μm.

Interestingly, the hit list contains bromosporine (hit #1 in Figure 3F), a broad-spectrum inhibitor for acetyl-lysine reader bromodomains (BRDs). Since transcriptional and epigenetic alterations have been extensively documented in HD (Bassi et al., 2017; Glajch and Sadri-Vakili, 2015; Sharma and Taliyan, 2015; Zuccato et al., 2010) but bromodomains have not previously been linked to HD pathologies, we chose to further investigate BRD inhibitors. We determined the half-maximal effective concentration (EC50) of bromosporine to be 118.8 nM and the half-maximal adverse effect score on the HD neuruloids as 760.7 nM (Figures 3G and 3H). There is therefore almost an order of magnitude in concentration separating efficacy and adverse effect. Moreover, when the adverse effect was quantified on the WT neuruloids, we found the half-maximal value to be 1.1 μM. Therefore, the effect of bromosporine at submicromolar concentrations is specific to the CAG expansion, as it shows efficacy on HD neuruloids while not affecting WT neuruloids.

Because bromosporine is a broad-spectrum inhibitor, with reported half-maximal inhibitory concentrations (IC50) of 0.41 μM, 0.29 μM, 0.122 μM, and 0.017 μM for BRD2, BRD4, BRD9, and CECR2, respectively, we then asked whether we could identify which target within this panel was responsible for the rescue. We used a chemogenic approach, taking advantage of the relatively large number of bromodomain inhibitors that were recently developed. Our panel consisted of 19 small-molecule inhibitors of BRDs, which we assayed at 10 μM for rescue of the HD neuruloids. The results of this BRD-inhibitor screen resulted in the identification of BI-7273, a specific BRD7/9 inhibitor, as the only inhibitor reversing the HD phenotype with low adverse effect at 10 μM (Figure 4A). Since BRD9 is the common target between our primary hit, bromosporine, and BI-7273, we then tested whether other BRD9-specific inhibitors could rescue HD neuruloids with low adverse effect for a range of concentrations. In total, we found that five BRD9 inhibitors—bromosporine, BI-7273, BI-9564, dBRD9, and I-BRD9—all showed submicromolar EC50 associated with low adverse effect (Figure 4B). This raises the interesting hypothesis that BRD9 is the relevant target in the reported rescuing effects.

Figure 4

Screen of bromodomain inhibitors and additional validation assays

(A) Target, rescue, and adverse effect for specific bromodomain inhibitors.

(B) Dose response for rescue and adverse effect of several bromodomain inhibitors (associated z-AUC values are shown in Figure S4D).

(C) Validation of efficacy of bromosporine and BI-7273 in 500-μm-diameter micropatterns.

(D) Rescue efficiency determined using individual channels.

(E) Rescue efficiency determined using a combination of channels (n = 350).

(F) Phenotypic space obtained from the latent vectors of an autoencoder confirms phenotypic rescue by compounds as quantified in (E) (n = 435).

(G and H) Validation in cortical neurons shows that BRD inhibitors rescue the effect of BDNF removal in HD neurons, with two concentrations tested per compound, 1 μM (blue) and 5 μM (red). Data points are individual fields of view collected in multiple wells; n = 185. p values were calculated using the two-tailed Mann-Whitney U test (Benjamini-Hochberg corrected); n.s. (not significant), p > 0.05; ∗∗∗p < 0.0005. Scale bar, 50 μm.

Bromodomain inhibitors rescue HD phenotypes in multiple assays

Finally, we sought to fully characterize the rescuing properties of bromodomain inhibitors in the HD neuruloids, as well as in secondary assays. We chose to focus on the broad inhibitor bromosporine as well as the more specific BRD9 inhibitor BI-7273. We first performed rescue experiments using our previously published neuruloid protocol on 500-μm-diameter micropatterns. We chose four markers that were shown to reveal the main phenotypic signature of HD neuruloids: collagen IV, PAX6, the neural crest marker SOX10, and phalloidin. We also included an additional HD cell line carrying a 72-length CAG repeat. Immunofluorescence images were analyzed using a classifier and autoencoder as described above. We found that both drugs rescued the 56CAG phenotype back to WT, confirming our screening results. Moreover, both drugs had a similar rescue effect on the 72CAG lines. Analyzing the individual contributions of the channels to the rescue, we found that while all individual channels by themselves show some degree of rescue (except for collagen IV), they also display greater variability in the rescue score. This is not unexpected, since we have previously found that the HD phenotype consists of subtle combinations of effects, such as the spatial expression of markers relative to each other. Consistently, training a neural network on WT and 56CAG lines using PAX6, phalloidin, and DAPI gives a better separation between WT and HD and shows that the drugs rescue both HD lines (Figures 4C–4E). The addition of collagen IV to these channels leads to a similar result (Figures S3C–S3E). Similarly, using the phenotypic space of the neuruloids obtained by an autoencoder, we find that this also indicates a complete rescue of both 56CAG and 72CAG phenotypes (Figure 4F). Rescue is also indicated by more conventional analysis of individual features such as the radius of the central PAX6 area (Figure S3F).

Currently, the most prominent strategy for HD therapeutics is based on reducing total or mutant HTT levels (Leavitt et al., 2020; Li et al., 2019; Marxreiter et al., 2020; Shannon, 2020; Tabrizi et al., 2019). To investigate the mechanism of action by which BRD inhibitors rescue the HD neuruloids, we asked whether they reduced the levels of HTT protein. Interestingly, we found that only bromosporine reduced HTT levels, but at higher concentrations than the EC50 measured on HD neuruloids (Figure S4). Therefore, the mechanism of action of iBRDs is not based on the reduction of HTT levels.

Finally, we asked whether iBRDs could show a positive effect in assays that more directly mimic neurodegeneration in HD. We thus turned to the susceptibility to apoptosis that can be measured in HD neurons upon brain-derived neurotrophic factor (BDNF) removal (Mattis et al., 2015). We differentiated 20CAG and 72CAG hESC lines toward cortical neurons for 45 days (Shi et al., 2012). BDNF was removed from WT and HD neuron culture media for the following 7 days. As previously reported, neurons grown from HD cells showed upregulation of apoptosis markers as measured by the TUNEL assay, in contrast to 20CAG neurons (Mattis et al., 2015). Strikingly, treatment with the BRD inhibitors bromosporine, I-BRD9, and BI9564 showed significant and in some cases complete reduction of apoptosis levels (Figures 4G and 4H), with no effect on MAP2+ projection length (Figure S4C). Interestingly, BI-7273 had no activity in this assay. Importantly, these results translate our findings from a developmental context to a paradigm of neurodegeneration and strengthen our findings that iBRDs represent a promising class of molecules that rescue a range of HD phenotypes in vitro, and justify further work aimed at exploring the potential of BRDs as new targets for HD.

Discussion

This study lays the groundwork for combining deep neural network and bioengineered human microtissues to carry out drug screens on neural organoids. This is achieved by the combination of a highly reproducible, scalable organoid platform allowing the easy generation of large image databanks required for leveraging the power of data analysis schemes based on deep neural networks. The method has a broad applicability and can be used with any kind of micropattern-based self-organized structure, potentially including 3D organoids that can be grown from a micropattern seed. We also foresee widespread applicability of this method with future developments that focus on standardizing organoid cultures.

We have found that disease phenotypes can be subtle, and we therefore require specifically engineered reproducible organoids to be able to identify such phenotypes reliably. This distinction would be more challenging in the first generation of brain organoids that mimic the complexity of multi-fate tissues, but with limited reproducibility. However, even with the limited organoid variability associated with micropatterning, phenotypic analysis using deep neural networks was required to meet the stringent statistical requirements of HTS. These powerful tools can extract all existing features that constitute a complex phenotype while simultaneously accounting for the natural variability of biological data, leading to significantly improved classification accuracy compared with methods that use pre-defined features. This method also abolishes the need to engineer specific features that may not be applicable in all situations. For example, approaches such as Cell Painting (Bray et al., 2016) contain a range of features that rely on the segmentation of cells, which we found to be impossible to perform reliably in an organoid context because of the dense and complex 3D arrangement of the cells. On the other hand, a current limitation of our approach is the difficulty of extracting the combination of features that a neural network learns. However, improving the interpretability of networks is an active area of research (Olah et al., 2018), and we have shown that simple approaches such as identifying the contribution of individual imaging channels can already give insight into which features are most useful for classification. In addition, we have introduced an unbiased phenotypic space using autoencoders. This approach offers great potential for obtaining deeper insight into the effect of diseases and molecules on phenotypes, since the directions in the latent space can be interpreted using the decoder part of the autoencoder. We have already successfully used these neural networks for determining the adverse effect of compounds, and a more detailed analysis of the phenotypic space is a promising avenue for future research.

The pleiotropic function of HTT is still far from being elucidated, and how mutant HTT leads to the disease pathogenesis is mostly unknown. Unfortunately, this is a common situation for most rare neurological diseases. The lack of therapeutic routes for HD is further exacerbated by an over-reliance on murine models, which for many neurological disorders translate poorly to the human. While the ideal pre-clinical discovery tool for HD would be the creation of a mature striatum showing neurodegenerative hallmarks integrated within an HTS screen, this is still far beyond our reach. Conversely, our method allows precise quantification, unbiased discovery of new mechanism of actions, and easy scalability in the context of an early developmental transition. While being far from modeling patients at the time of clinical onset, this has many advantages. (1) The body of evidence pointing toward a developmental component of HD is increasing, and with it the hypothesis that the disease has its roots in development; (2) this allows convenient modeling of the pre-disposition state of HD implied by the ubiquitous presence of mutant HTT from fertilization onward. If we could normalize mutant HTT functions in the pre-manifest phase, we could potentially push back clinical onset. (3) Our approach allows integrating many of the unknown functions of HTT and mutant HTT within a self-organized structure that integrates many biological phenomena and results in a phenotype easily learnable by deep neural networks. This has led to the agnostic discovery of a new class of molecules, iBRDs, which not only revert the effects of mutant HTT in developmental organoids but also reduce neurotoxicity induced by the HTT CAG expansion in vitro. This highlights how an unbiased phenotypic approach can deliver molecules efficient at rescuing many of the core toxic functions of mutant HTT without relying on a priori hypotheses. Moreover, a crucial advantage of our approach over existing pharmacogenomic methods that focus on transcriptomic endpoints such as the Connectivity Map (Lamb et al., 2006) or the Library of Integrated Network-based Cellular Signatures (Stathias et al., 2020) is that organoids reproduce the complex interactions of multiple tissues and the associated interplays of multiple signaling pathways and morphogenetic processes and, hence, represent a more integrated biological readout.

There are currently no approved cures for HD, and clinical strategies exclusively aim at reducing HTT levels. This lack of alternative mechanisms of action highlights our poor understanding of HD, as a significant body of research also points toward HD as being at least partially due to a loss of function, or a more complex interplay between mutant and WT HTT. Overall, this lack of therapeutic handles stresses a dramatic need for new discoveries. Here, using a completely unbiased approach, we reveal that BRD inhibitors should be further explored and validated as therapeutic target for HD.

BRDs are the principal readers of ε-N-acetyl-lysine marks. They therefore have an important role in the targeting of chromatin-modifying enzymes to specific genomic sites (Muller et al., 2011). Progressive transcriptomic dysregulation is an established hallmark of HD, with epigenetic deregulation as determined by DNA methylation or histone modifications becoming a prevailing feature. In this context, multiple histone deacetylase inhibitors have shown beneficial effects in HD mouse models (Jia et al., 2015; Muller et al., 2011; Siebzehnrubl et al., 2018; Thomas et al., 2008). Our discovery extends these findings by showing the potential of epigenetic readers for HD besides mark erasers such as histone deacetylase. Moreover, the BRD family is highly druggable, and iBRDs that have been extensively developed for cancer-related applications and are thought to carry important potential for neurological disorders have been poorly investigated. In this context, the BET BRD JQ1 has been recently shown to provide no therapeutic benefit in mouse models of HD (Kedaigle et al., 2020). However, the main target of JQ1 is BRD4, and our results point toward other family members such as BRD9, with the limitation that no genetic validation of these targets was performed in our study that could disentangle potential small-molecule off-target effects and would conclusively identify the precise family member responsible for the rescue. We therefore envision specific iBRDs to drive exciting new therapeutic developments for HD.

In the future, the neuruloids used here will be replaced by more mature and differentiated neuronal structures with high reproducibility. An exciting prospect is to also incorporate additional readouts, such as transcriptomic data, into the phenotypic space. This has the potential to lead to the generation of novel hypotheses beyond classical pharmacology and has the potential to greatly expand our insight into the mechanisms and possible treatments of a wide range of diseases.

Limitations of the study

We were not able to reliably segment nuclei in the neuruloids because their density is much higher than in typical 3D organoid structures. We have therefore limited ourselves to benchmarking our results against tools that do not necessitate nuclear segmentation. For validation of the screening results, there are further approaches that we have not undertaken here, e.g., western blots or qPCR for key cell types. At the level of our biological findings, BRD9 has not been validated by other means than small-molecule-based inhibition; a genetic knockdown approach would strengthen the validation of this target in our assays. Finally, given that no animal work has been performed here, the validity of BRD9 as an HD target remains to be observed in vivo.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Fred Etoc (fred@rumiscientific.com).

Materials availability

The HTT overexpression line generated during this study is available from the lead contact upon request.

Experimental model and subject details

Cell lines and clones

All experiments were performed using genetically engineered clones derived from the female RUES2 (NIH #0013) parental cell line, which was created in the Brivanlou lab and is listed in the NIH Human Embryonic Stem Cell Registry. The RUES2-HD allelic series used in this study was characterized and published previously (Ruzo et al., 2018).

Cell culture

All hESC lines were grown in HUESM medium that was conditioned with mouse embryonic fibroblasts and supplemented with 20 ng/mL bFGF (MEF-CM). HUESM is the growth media that we used and differs from the standard H1 medium described by Thomson’s group in the use of DMEM rather than DMEM/F12. Cells were tested for mycoplasma at 2-month intervals using the MycoSEQ Mycoplasma Detection System (ThermoFischer). Cells were grown on tissue 10 cm2 culture dishes coated with 1.5 mL Geltrex (Life Technologies) solution diluted 1/100 from stock in DMEM/F12 for cell maintenance and propagation, and a 1/50 dilution for monolayered corrtical neural differentiation.

Micropatterned cell culture

Micropatterned glass cover slips (CYTOOCHIPS™ Arena A, Arena 500 A, Arena EMB A) were first coated with 10 μg/mL of recombinant Laminin-521 (BioLamina, LN521-03) diluted in PBS+/+ (+/+ meaning with Calcium and Magnesium, Gibco) for 3 h at 37°C. Micropatterns were placed face-up on to Parafilm in a 10 cm dish, then 800 μL laminin solution was added to the micropattern. After 3 h at 37°C, the coated micropattern was transferred to a 35 mm dish with 5 mL of PBS+/+. Laminin was then removed with six serial dilutions in PBS+/+ (dilution 1:4) before two complete washes in PBS+/+. The coated micropattern was kept in PBS+/+ at 37 °C. Cells were seeded as follows: cells growing in MEF-CM on culture dishes were washed once with PBS-/- (Gibco), then treated with Accutase (Stem Cell Technologies) for 5 min. Cells were then triturated with a pipet to ensure single cell suspension and Accutase was diluted out with 4× HUESM medium supplemented with 20 ng/mL bFGF and Rock inhibitor Y27632 (10 μM, Abcam ab120129). Cells were then further diluted with the same medium and 5 × 105 (or as indicated) cells in 3.0 mL of medium were placed over the micropattern in a 35 mm tissue culture dish, then incubated at 37°C. After 3 h, the micropattern in a dish was washed once with PBS+/+ and then replaced with HUESM with 10 μM SB431542 and 0.2 μM LDN 193189. At Day 3, the medium was replaced with HUESM with 10 μM SB431542 and BMP4 (10 ng/mL). At Day 5, the medium was replaced with the same fresh medium and incubated until day 7.

Method details

Immunofluorescence

Micropattern coverslips were fixed with 4% paraformaldehyde (Electron Microscopy Sciences 15713) in warm medium for 30 min, rinsed three times with PBS-/-, and then blocked and permeabilized with 3% normal donkey serum (Jackson Immunoresearch 017-000-121) with 0.5% Triton X-100 (Sigma 93443) in PBS-/- for 30 min. Micropatterns were incubated with primary antibodies for 1.5 h, washed three times in PBS-/- for 5 min each, incubated with secondary antibodies conjugated with Alexa 488, Alexa 555, Alexa 594 or Alexa 647 (1/1000 dilution, Molecular Probes) and 10 ng/mL of 4′,6-diamidino-2-phenylindole (DAPI, Thermo Fisher Scientific D1306) for 30 min and then washed two times with PBS-/-. For double staining with antibodies from the same species, Alexa 488 Fab fragments (Jackson Immunoresearch, 715-547-003) and Fab fragment IgG (Jackson Immunoresearch, 715-007-003) were used. Coverslips were mounted on slides using ProLong Gold antifade mounting medium (Molecular Probes P36934).

Imaging

All confocal images were acquired on a Zeiss Inverted LSM 780 laser scanning confocal microscope with a 10X or 20X water-immersion objective. 4 equally spaced Z planes spanning the height of each neuruloid were acquired and their maximum projection used as input for the neural networks.

High-throughput screening

We adapted our protocol for neuruloids formation on micropatterned chips to 96 well plates. Briefly, one day before cell plating, micropatterned 96 well plates with 700 μm disks (CYTOO) were coated with laminin for 3 h using a wellmate dispenser and a solution of 100ul per well. Plates were then washed 3 times with 200 μL PBS+/+ per well, wrapped in parafilm and stored overnight at 4 degrees. On the day of cell seeding, a single cell suspension at a density of 0.18 M cells/mL was prepared with the same specification as in our protocol for chips. The final volume of the resuspension was calculated to be 20 mL per plate and adding an extra 20 mL to account for the dead volume of the cell dispenser. For seeding, the wellmate was used for aspirating the PBS+/+ from the plates and dispensing 200 μL/well of cell suspension. Just after seeding, the plates were placed on a flat bench and laterally shaken on the surface in a cross-like shape to homogenize the cell suspension in each well after which the plates were left untouched on the bench for 1 h for the cell suspension to homogeneously sediment and attach on the glass substrate. After 1 h, the plates containing cells were returned to the incubator for 3hrs. Following cell seeding, the plates were processed according to our neuruloids differentiation protocol in an EL406 robot allowing to process 12 plate at a time for both media changes at days 3 and 5 of differentiation, but also immunostaining. During media changes, compounds were applied at a 10 μM final concentration using an automated Biomek FXp (Beckman Coulter) equipped with dual (96/384-well and Span8) pipetting arms and pin tool (VP Scientific) for compound pinning form the compound library plate to the assay plates, or manually using a manual floating pinning tool from V&P Scientific. After staining, plates were sealed with black stickers and kept in the fridge before imaging with an InCell Analyzer 2000 high-content imager and 4X lens allowing to image full wells with four tiles. A single Z plane was acquired, using the auto-focus on the PAX6 channel. After acquiring all images, individual wells were stitched and individual colony extracted using a custom Matlab code described previously (Etoc et al., 2016). For each experiment realized in 96 well plates, 2 control plates of untreated 56CAG and 2 control plates of untreated 20CAG were used in addition to the assay plates that were contacted with compounds. Moreover, the wells on the side of each treated plates were left unperturbed as control wells. Hit molecules are: hit#1: Bromosporine, hit#2: BI-D1879, hit#3: MK-2461, hit#4: Ponatinib, hit#5:

Neuronal differentiation and BDNF assay

Isogenic lines were subjected to a default neural induction protocol adapted from a protocol published previously (Shi et al., 2012). In short, cultures were seeded in confluent adherent culture, at 1×105 cells/cm2 and fed every day with N2B27 serum-free medium with TGFβ inhibition (10 μM SB-431542, Sigma and 0.2 μM LDN-193189, StemGent) applied for the first 10 days and 5 ng/mL FGF8 (R&D Systems) from day 12–22 to maximize the frequency of CTIP2 (BCL11B)-positive cells, dissociated and seeded on adherent substrate (polyornithine/laminin, Ibidi) at day 14, continued in culture with added BDNF (R&D Systems, 10 ng/mL), IGF1 (R&D Systems, 10 ng/mL), cAMP (1 μM, R&D Systems) and ascorbic acid (200 μM, Sigma), until dissociation and reseeding on a 96 well plate with glass bottom at day 40. At this stage, for BDNF removal assays, samples were kept in medium with or without all growth factors removed for 7 days together with or without small molecules, before fixation and analysis at day 47. The TUNEL assay was then performed with a kit from Abcam (ab66110) and counterstained with DAPI and MAP2. Multiple images from 4 separate wells per condition were taken at 20× magnification (exact numbers listed in the legend of Figure 4H). TUNEL fluorescent intensity was measured at the location of nuclei as indicated by the DAPI stain. DAPI-positive nuclei were identified using Otsu thresholding implemented in scikit-image, and TUNEL intensity was summed at the location of DAPI-positive pixels and normalized by the DAPI intensity. TUNEL intensity was then normalized to the mean of the untreated controls.

Quantification and statistical analysis

Quantification of HTT and mutant HTT levels

We performed the neuruloid experiments in 96 well plates treated with compounds or untreated in control wells. At the end of the experiment medium was removed, samples were washed with PBS once, before removing any liquid left form the cells, freezing the plates in dry ice and shipping to Charles River Laboratory for mutant and WT HTT levels quantification using the MSD method (Macdonald et al., 2014).

Neural network analysis

For both data cleaning and classification, we used a well-established image classification pipeline based on pretrained residual neural networks (He et al., 2016) as implemented in the Fastai library (Howard and Gugger, 2020), which is based on the PyTorch machine learning framework (Paszke et al., 2019). For all applications, we found that already the smallest of the standard residual network architectures (resnet18) gave excellent classification results. We used 3-channel images and normalized them according to the same normalization that was used in the original residual network training performed on the Imagenet dataset (mean and standard deviation per channel (0.485, 0.456, 0.406) and (0.229, 0.224, 0.225), respectively). The last layer of these pre-trained networks was removed and replaced with the default untrained combination of the Fastai library (concatenated adaptive average pool and adaptive maximum pool, batch norm, drop-out and additional dense layers that reduce to the final two output nodes). The training data was split randomly into a training and validation set with a ratio of 70% to 30%. For the phenotypic screen, control data consisted of a total of 4879 images (2333/2546 WT/HD), which were split into 3415 (1633/1782 WT/HD) training and 1464 validation images (700/764 WT/HD). The trained models were then applied to a total of 32173 images from the compound plates. Models were first trained with all layers frozen except for the untrained last layers, using one-cycle learning (Smith, 2018), a maximum learning rate of 0.003 and the Adam adaptive learning algorithm (Kingma and Ba, 2015). Subsequently, all layers were unfrozen and trained with differential learning rates between 10^-6 and 10^-4. Training was performed until the training loss did not decrease further. Overfitting was avoided by using dropout layers and data augmentation, including random flips, rotations, zooms, contrast changes and affine transformations. The training accuracy did not exceed the validation accuracy, indicating that our models were not overfitting.

For data cleaning, we chose to first train a neural network on images of nuclear DAPI stains that would recognize well-formed organoids and would filter out images with no cells or with detached or merged organoids. For this, we prepared a training set of 3698 annotated images and trained a neural network as described above. The results, shown in Figure S1, allowed us to efficiently pre-process all similar datasets. For data cleaning of the 500 μM neuruloid data presented in Figures 4C–4F, we directly used the latent space of the autoencoder to detect abnormally formed colonies. These show up as a separate cluster in the latent space and can easily be removed from the analysis. We also found that, because of subtly varying imaging conditions for the different samples in Figures 4C–4F, matching the intensity histograms of images to the histogram of a randomly chosen reference image improved the classification accuracy. This was done using the match_histograms function implemented in scikit-image.

To obtain the phenotypic space, we implemented a custom autoencoder. Since we could not use pretrained networks as for the classifier and had to train the network from scratch, we chose a small architecture consisting of two convolutional layers (with 32 and 64 filters respectively, and associated ELU activation, max pooling and batch normalization layers). The results after those operations were mapped to the latent space, after which the data passed through a similar inverse architecture (two transposed convolutional layers with 64 and 32 filters, respectively) and a final sigmoid layer. The autoencoder was trained using the mean squared error and Adam algorithm with a learning rate of 2e-5. We found that the dimensionality of the latent space did not have a significant effect (Figure S2).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
PAX6 Rabbit Polyclonal	BioLegend	901301
PAX6 Mouse monoclonal	BD Biosciences	561462
SOX10 Goat Polyclonal	R&D	AF2864
COL4 Rabbit Polyclonal	Abcam	ab6586
MAP2 Chicken Polyclonal	Abcam	ab5392
Chemicals, peptides, and recombinant proteins
Laminin-521	BioLamina	LN521-03
Rock inhibitor Y27632	Abcam	ab120129
Accutase	Stem Cell Technol.	07920
Triton X-100	Sigma	93443
SB-431542	Sigma	s4317
LDN-193189	StemGent	04–0074
BDNF	R&D	248-BDB
IGF1	R&D	291-G1
Experimental models: Cell lines
RUES2 stem cell line	Rockefeller University	NIH#0013
RUES2-HD allelic series	Ruzo et al., (2018)	N/A
RUES2-HTT KO with WT overexpression	This paper	N/A
Software and algorithms
Organoid screen analysis software	This paper	https://github.com/jjmetzger/organoid_screen_analysishttps://doi.org/10.5281/zenodo.7009263
Scikit Image		https://scikit-image.org/