FATE-MAP predicts teratogenicity and human gastrulation failure modes by integrating deep learning and mechanistic modeling

Domain: chemistry
Published: 2026-02-19
Editorial outcome: minor_revision

Abstract
Gastrulation, a critical developmental stage involving germ layer specification and axes formation, is a major point of failure in human development, contributing to pregnancy loss and congenital malformations. However, due to ethical constraints and anatomical differences in animal models, the failure modes underlying human gastrulation remain poorly understood. To elucidate these failure modes, we introduce FATE-MAP (Failure Analysis and Trajectory Evaluation via Mechanistic-AI Prediction), an integrated platform that combines high-throughput perturbations of human 2D gastruloids with quantitative phenotypic mapping, predictive deep learning, and mechanistic morphogen modeling. Analyzing over 2,000 drug-treated human 2D gastruloids, we mapped a phenotypic morphospace that separates canonical patterning, in which primitive-streak fates are correctly specified and radially organized, from failure modes, defined as departures from this organization and marked by a loss of a required fate and/or radial symmetry. To predict and interpret patterning outcomes, FATE-MAP combines a transformer linking chemical structure to phenotype with PDE simulations of morphogen transport and cell fate specification, and projects both outputs onto the experimentally defined morphospace. Applying this framework, we flagged two clinical molecules as potential teratogens and identified two parameters, cell density and SOX2 stability, that form orthogonal morphospace axes along which canonically patterned gastruloids systematically vary. FATE-MAP thus provides a roadmap for decoding human developmental trajectories and accelerating safe therapeutic discovery.

Introduction
Human gastrulation, which involves two symmetry-breaking events and the specification of the germ layers, has been hypothesized to act as a go, no-go decision for the embryo, serving as an early test of the fidelity of developmental processes that are essential for proper formation later.1,2 As many as 30% of human pregnancies are predicted to fail at this critical stage,3,4 yet because of anatomical differences between human and other vertebrate models of gastrulation,5 as well as ethical constraints on experimentation with human embryos,6 we lack a comprehensive view of the key factors that govern human gastrulation and how perturbations to these factors lead to failure.

Inborn errors of gastrulation can lead to either spontaneous abortion or congenital malformations.7,8 Notably, congenital malformations cause more deaths than all pediatric cancers9 and have a higher mortality rate than that of stroke, Alzheimer’s disease, and diabetes combined.10,11 Despite this societal burden, little progress has been made in prevention or prediction, with 50–80% of congenital malformations having no known cause.12,13 Teratogenicity testing has traditionally relied on exposing model organisms, such as rats, mice, and zebrafish, to potential teratogens and assessing gross morphological changes. Critically, such models do not capture many human-specific teratogens as mice and rats are notoriously resistant to thalidomide-induced birth defects,14,15 highlighting the need for quantitative methods to predict even minor harm to the developing embryo.

To address these challenges, we developed FATE-MAP (Failure Analysis and Trajectory Evaluation via Mechanistic-AI Prediction), an integrated experimental and computational pipeline that combines high-throughput drug screening with phenotypic mapping, mechanistic modeling, and deep learning models to identify the failure modes of human gastrulation, uncover their mechanistic origins, and improve our ability to forecast human-specific developmental risks. Central to this approach is the construction of a morphospace,16,17 a low-dimensional, quantitative map of morphological phenotypes, that captures both canonical and disrupted trajectories of human gastrulation. Here, we use morphospace specifically to describe the diversity of 2D patterning outcomes in micropatterned stem cell colonies, not to imply full embryonic morphology or organization.

Generating such a map requires an experimental model system that faithfully recapitulates early human development while enabling high-throughput perturbation and analysis. Recent advances in stem cell–based embryo models provide this opportunity, 18 offering experimentally tractable systems to probe human developmental trajectories under defined conditions.

While 3D models such as blastoids and synthetic embryo models (SEMs) can recapitulate aspects of early development, they often suffer from variability and low formation efficiencies.19-25 On the other hand, 2D or 2.5D models that leverage lithographic micropatterned surfaces achieve near uniform formation of 103-104 constructs per experiment,26-28 offering highly scalable and reproducible systems. These models accurately reproduce the disc-like geometry of the epiblast, albeit at the cost of altering the topology of surrounding extraembryonic tissues such as the amnion. Among these, the 2D gastruloid, a micropatterned disc of embryonic stem cells, serves as a reliable model of the human primitive streak, differentiating into radially patterned germ layers and recapitulating all major cell types of the gastrulating human embryo.29 Leveraging this platform, we generated an experimental morphospace from ~2,000 gastruloids exposed to a diverse library of small molecules targeting developmental pathways. Using nonlinear dimensionality reduction, we embedded the resulting phenotypic profiles into a low-dimensional map, where both canonical and disrupted trajectories formed structured regions of morphospace, revealing discrete failure modes marked by symmetry loss or germ layer absence.

Although nonlinear dimensionality reduction effectively organizes phenotypic diversity and reveals failure modes, it lacks interpretability: the approach offers no direct path for predicting outcomes from molecular inputs, and the developmental parameters that determine a colony’s position in morphospace remain obscured. FATE-MAP addresses this limitation through a two-pronged modeling strategy. First, we developed a transformer-based neural network that predicts morphological phenotypes directly from chemical structure. This model enables in silico forecasting of developmental risk for untested compounds and provides a scalable, prospective approach to identify human-relevant teratogens based on molecular features alone. However, this chemical structure-to-phenotype model operates as a black box; it maps inputs to outputs without explicitly modeling the developmental processes that generate phenotypes. To provide mechanistic insight, we developed complementary partial differential equation (PDE)-based models of morphogen diffusion and cell fate transitions that simulate patterning outcomes from first principles.

To unify these predictive and mechanistic approaches, we trained a neural network to embed simulated phenotypes from our mechanistic model and deep learning model into a shared, experimentally defined morphospace. This integrated framework allows direct comparison between model predictions and observed outcomes, enabling both validation of predictions and mechanistic interpretation of developmental failure modes. Together, these complementary strategies demonstrate the power of combining chemical structure-based modeling with biologically grounded simulations to enhance predictive accuracy and reveal the underlying rules of developmental fate patterning.

Results
A high-throughput platform for 2D gastruloid phenotyping
We sought to leverage the 2D gastruloid model in a high-throughput screen to map the morphospace of human gastrulation based on prior demonstrations that the model mimics the cell fate decisions of human gastrulation and produces well-defined spatial regions of differentiation that correspond to cell fate specification of the primitive streak (Fig. 1a-b). We chose to perturb BMP4-initiated gastruloids with compounds from a library of 210 drugs annotated for their activity against stem cell signaling pathways (Fig. 1c). While it is possible that a given drug may have multiple protein targets, our approach is target-agnostic and only requires a large diversity of perturbations to capture the variation in outcomes. To quantify phenotypes, we chose to use immunofluorescence staining of germ layer markers: GATA3 for amniotic ectoderm, Brachyury (BRA) for mesoderm, and SOX2 for the undifferentiated epiblast. We imaged ~10 colonies per drug condition and used a custom image segmentation algorithm to identify the levels of each cell fate marker in every nucleus of each colony (Supplementary Fig. 1a), collecting both cell fate and location for ~2 million cells across 2,025 colonies.

We leveraged the fact that colonies are generally radially-symmetric to build a data structure for downstream pattern analysis. Inspired by approaches that vectorize animal behavior,30 we reasoned that we could meaningfully compress colony morphology through averaging cell fates over a set of 50 concentric annular bins, each ~5 µm in width, ranked by their position, from the edge to the center of each colony (Fig. 1d). This yielded a 150-dimensional vector for each colony containing the average azimuthal signal for GATA3, BRA, and SOX2, as a function of colony position (Fig. 1e). Although this embedding does not explicitly preserve geometric distances between adjacent bins, it preserves the exact radial ordering, enabling spatially coherent comparison of radial expression patterns across colonies (Supplementary Fig. 1b-c). As controls, we included several BMP4-only and untreated (no-BMP4) colonies from multiple independent wells, as well as a Wnt-activating (CHIR-98014) positive control which, as expected, differentiated the entire colony into mesoderm (Supplementary Fig. 1d). These controls demonstrated that patterning phenotypes were reproducible across multiple wells and days of running the assay (Supplementary Fig. 1e). Across the broader compound screen, plate-to-plate variability was minimal compared to drug-induced effects (Supplementary Fig. 1f), confirming that the 2D gastruloid model analyzed in this way meets the criteria for constructing a morphospace map.

Human gastruloid morphospace
Upon visual inspection of the drug-treated gastruloids, we noticed a variety of deviations from the canonical patterning observed in the controls. We sought to develop an unsupervised algorithm that would (1) project each colony down to two dimensions, and (2) identify groups of perturbations that result in similar phenotypes so that they could be evaluated as potential failure modes. Inspired by previous methods,31 we coupled t-SNE to an unsupervised clustering procedure which applies watershed segmentation to the two-dimensional embedding after convolution with a Gaussian kernel (Supplementary Fig. 1g; Methods). This method does not impose a fixed number of clusters, as other clustering algorithms do (k-means, gaussian mixture models, spectral clustering, etc.), and thus produces clusters that reflect the local variations in phenotype without supervision (Fig. 1f, top). We identified 12 clusters (C1-C12) which we call morphological phenotypes.

To visualize the phenotype of each cluster, we generated composite images by taking the mean pixel intensity for each cell fate marker from all colonies in each cluster (Fig. 1f, bottom). This allowed us to build a coarse-grained understanding of colony patterning across the entirety of morphospace. Clusters could be divided into two major categories. The canonical patterning region (C2–C10) contains colonies that retain all three fate markers in concentric radial order (SOX2 central, BRA intermediate, GATA3 peripheral) and differ primarily in the relative sizes of these domains. BMP4-only controls map exclusively within this region (C2, C3, C4, C5, C8, C9; Supplementary Fig. 2a). In contrast, clusters C1, C11, and C12 lie outside the canonical region and exhibit either the loss of a required cell-fate marker (C1: loss of SOX2; C11: loss of BRA) and/or loss of radial symmetry (C12). Since either of the aberrations captured in these clusters would result in a failure of gastrulation, we classify these clusters as failure modes.

Because failure modes can arise from either differentiation failure or compound-induced cytotoxicity, we quantified total nuclear counts for every colony and identified a small subset of drugs that reduce viable cell number below 50% of the mTeSR-only control (Supplementary Fig. 3). This analysis separates genuine patterning failures from cytotoxicity-driven collapse within the failure-mode regions. Representative drug treatments from each cluster are shown in Fig. 1g. For more examples, see Supplementary Fig. 2b, and for the entire dataset, visit https://max-wilson.mcdb.ucsb.edu/research/gastruloid-morphospace.

Predicting morphological phenotype from chemical structure
Having defined a morphospace of human gastrulation and identified discrete failure modes, we next sought to develop a predictive model capable of mapping a chemical structure directly to morphological phenotype and morphospace location. To do so, we designed a two-step modeling framework in which chemical structure is first converted into a predicted phenotypic vector, and then embedded into the low-dimensional morphospace (Fig. 2a). Each step is performed by a separately trained neural network.

To build the chemical structure-to-phenotype model (Fig. 2a, left), we began by embedding each compound’s SMILES representation using ChemBERTa, a pretrained transformer model trained on ~10 million PubChem compounds.32 This yielded a 768-dimensional molecular feature vector for each drug, 𝑥⃑ , which we used as input to a neural network, 𝜙(𝑥⃑), that was trained to predict the 150-dimensional phenotype vector of the gastruloid (Methods). Ten percent of the compounds were randomly held out as a blinded test set to assess generalization.

To map phenotype vectors into morphospace, we trained a second neural network to approximate the relationship between the 150-dimensional phenotype profile, 𝑦⃑ , and its corresponding t-SNE coordinates (t-SNE , t-SNE ). While t-SNE lacks a formal inverse function, neural networks can be used to learn the embedding (i.e., parametric t-SNE). We trained a separate neural network to minimize the distance between predicted and true t-SNE coordinates based on experimental data (Fig. 2a, right; Methods). The resulting network, 𝜓(𝑦⃑), accurately reconstructed the original experimental morphospace and was used for all subsequent projections (Supplementary Fig. 2c).

The combined model accurately predicts 150-dimensional phenotype vectors and corresponding morphospace positions for both the training and blinded compounds, capturing canonical and failure-mode clusters. Across 188 training compounds, the mean RMSE was 0.17, with a mean cosine similarity of 0.90. Importantly, on average 74% of features per compound fell within one standard deviation of the empirical replicate distribution, indicating that most predicted feature values lie inside the observed biological variability. On the blinded set (n = 21), performance remained strong (mean RMSE = 0.24; mean cosine similarity = 0.85), with 54% of features per compound within one standard deviation of replicate variation. Aggregating all 209 compounds, the overall means were RMSE = 0.17, cosine similarity = 0.90, and 72% within one standard deviation. Representative predictions for all 12 clusters are shown in Fig. 2b, including both phenotype profiles and projected morphospace coordinates.

Model performance was maintained despite considerable chemical dissimilarity between the blind and training compounds (Supplementary Fig. 4a–b), indicating robust generalization to structurally diverse molecules. To explore this further, we clustered all experimentally screened compounds based on structural similarity and compared these groupings to morphology-defined phenotypic clusters (Supplementary Fig. 5). Despite the presence of chemically coherent scaffold-based clusters, we observed minimal spatial overlap between structure- and phenotype-defined groupings. This suggests that even compounds with similar scaffolds can produce distinct developmental effects, reinforcing the need for neural network–based approaches to uncover structure–function relationships that scaffold similarity alone may obscure.

Because the most consequential classification task is distinguishing failure-mode phenotypes from canonical patterning outcomes, we evaluated performance at two levels: (i) strict 12-cluster classification, which requires assigning each compound to the correct morphospace cluster, and (ii) binary failure-mode discrimination, treating the task as failure-mode vs. canonical patterning classification. Across the 50 compounds whose experimental centroids fall in failure-mode clusters, the model predicted all 50 to reside within the failure-mode region, yielding a binary failure-mode sensitivity of 1.00. Within this region, strict cluster identity was correct for 39 of 50 compounds (78%), with most misassignments reflecting shifts toward C11, the densest failure-mode attractor in morphospace. Overall binary classification accuracy was 0.71 (0.72 in the training set and 0.67 in the blind set). These results indicate that the model is highly reliable for detecting whether a compound induces a failure-mode phenotype, even if the precise sub-cluster label within the failure-mode region is occasionally uncertain. Of note, drugs were assigned to clusters not by visual inspection of their position in morphospace, which can be ambiguous, but using the same watershed-derived morphospace boundaries defined during the original clustering, ensuring consistent, data-driven classification across experimental and predicted phenotypes (see Supplementary Note 1 for details on the morphospace classification boundary and thresholding approach).

FATE-MAP predicts teratogenic risk from chemical structure
To evaluate the translational potential of FATE-MAP for predicting teratogenicity, we first examined drugs from our initial screen with known pregnancy-related risks. Several established teratogens,33-35 including isotretinoin, TTNPB, and bexarotene, were experimentally determined to fall within cluster 11 and shared a common morphological phenotype characterized by loss of mesodermal identity (BRA expression) (Fig. 3a). These colonies frequently displayed lighter central regions that reflect vertical protrusions caused by increased cell density under micropattern confinement (Supplementary Fig. 6). The observation that known teratogens induce failure-mode phenotypes suggests that morphospace location can be leveraged to identify additional high-risk compounds directly from chemical structure.

To benchmark this structure-to-phenotype framework against established teratogenicity assays, we assembled a curated panel of 33 FDA-approved drugs with well-characterized teratogenicity outcomes (positive or negative), supported by human clinical data (as reported in FDA labeling) and/or in vivo embryo–fetal development (EFD) studies in rats or rabbits, and for which results from both the TeraTox iPSC multilineage differentiation assay and the mouse embryonic stem cell test (mEST) have been reported.36 For each compound, FATE-MAP first generated a 150-dimensional predicted phenotype vector from its chemical structure and then projected this vector into the experimentally defined morphospace using the neural network approximation of the t-SNE embedding (Fig. 3b). Teratogenicity was inferred directly from morphospace location using the same cluster boundaries applied to experimental data: compounds whose predicted embeddings fell within the failure-mode region (clusters C1, C11, and C12) were classified as teratogenic, while those mapping to canonical patterning clusters (C2–C10) were classified as non-teratogenic (Fig. 3c). Notably, only 4 of the 33 benchmarking compounds were included in our experimental dataset. Across the full 33-drug panel, FATE-MAP achieved an overall accuracy of 72.7%, outperforming the mEST assay (51.5%) and matching the TeraTox assay (72.7%), while exhibiting higher balanced accuracy (74.8% for FATE-MAP vs. 66.3% for TeraTox) (Fig. 3d, Supplementary Data 1). Because FATE-MAP assigns phenotypes to clusters using fixed morphospace boundaries derived from experimental data, classification does not require re-testing each compound, enabling scalable, fully in silico teratogenicity prediction.

To further assess predictive performance relative to existing in silico developmental toxicity models, we evaluated a set of 20 compounds held out from training across three models: FATE-MAP, the CAESAR developmental toxicity predictor, and the Procter & Gamble (P&G) teratogenicity model. For CAESAR and P&G, we used the reported model outputs provided in their published evaluations,36 ensuring a direct comparison without retraining or re-tuning either model. Across this compound set, FATE-MAP achieved an accuracy of 75.0%, matching the CAESAR model (75.0%) and substantially outperforming the P&G model (35.0%), while also demonstrating the highest balanced accuracy (85.7% for FATE-MAP vs. 81.2% for CAESAR) (Supplementary Fig. 7a–b; Supplementary Data 2). Together, these results show that FATE-MAP captures developmental toxicity signals that generalize beyond the training distribution and performs competitively with both experimental and computational state-of-the-art approaches.

To evaluate FATE-MAP’s ability to prospectively identify teratogens, we selected two clinical-stage compounds, ulixertinib (an ERK1/2 inhibitor) and nirogacestat (a γ-secretase inhibitor), from our initial screen. These compounds were both experimentally confirmed to reside in cluster 11 (Supplementary Fig. 8a-b) and had no established developmental risk data in humans. Importantly, they were also excluded from model training, yet were still accurately predicted to fall within cluster 11 based solely on chemical structure (Fig. 3e). Notably, the predicted failure-mode phenotypes for both compounds are consistent with the known essential roles of ERK and Notch signaling in embryogenesis, where disruption of either pathway leads to severe developmental defects and embryo–fetal loss in animal models.37,38 While direct validation in human embryos is not possible, we applied a zebrafish-based developmental toxicity assay as an independent in vivo model to test these predictions, alongside isotretinoin as a positive control. Zebrafish are widely used in pharmaceutical developmental toxicology pipelines and offer higher-throughput, cost-effective testing compared to traditional mammalian models. Prior studies have shown that zebrafish assays achieve comparable or superior predictive performance to mouse-based testing (74.19% accuracy and 87.50% sensitivity for zebrafish versus 67.74% accuracy and 75.00% sensitivity for mice).39 At sub-micromolar to micromolar doses, all three compounds induced characteristic teratogenic phenotypes, including craniofacial defects, body axis shortening, and yolk sac edema (Fig. 3f, Supplementary Fig. 8c-h). Although ulixertinib and nirogacestat produced the expected posterior truncation in zebrafish, axial tissues such as somites and notochord remained present, consistent with partial, rather than complete, disruption of primitive streak–derived mesoderm and reflecting species-specific differences in pathway sensitivity. Overall, these results confirm that FATE-MAP can prospectively identify previously unrecognized teratogens, offering a scalable in silico alternative to in vivo developmental toxicity testing

Cell density shapes mesoderm positioning and fate
We next turned our attention to determine the causes of variance in the canonical patterning region. From previous observations on the importance of 2D gastruloid density in determining patterns,28,40-42 we hypothesized that mesoderm band width could be regulated by cell density. Because cell seeding produces inherent variation in cell density across colonies, BMP4-treated controls provide a natural spectrum for analyzing density-dependent effects. Indeed, the distribution of these control colonies in morphospace is explained by their density (Supplementary Fig. 9a), with a moderate correlation between density and mesoderm peak position (R2 = 0.69). Ranking colonies by density further revealed a clear expansion of mesoderm toward the colony center as density decreased (Supplementary Fig. 9b). This trend also held across cluster averages: lower-density clusters showed more central mesoderm expression, with a stronger correlation (R2 = 0.96) when plotted by cluster mean (Supplementary Fig. 9c). Ranking clusters by edge-to-center BRA profiles reinforced this relationship (Supplementary Fig. 9d), indicating that cell density explains a major source of phenotypic variance.

Given the impact of cell density on patterning, we hypothesized that it might be strong enough to override the effects of drug treatment. If correct, we would expect colonies to be rescued back into the canonical patterning region by increasing density or driven into failure modes by decreasing density. Supporting this, we observed (-)-Blebbistatin-treated colonies, typically found in failure mode clusters C1 and C12, regained canonical phenotypes at higher densities (Supplementary Fig. 9e, left). Conversely, DEL-22379–treated colonies, normally found at the edge of the canonical patterning region, were pushed into failure mode C1 at lower densities (Supplementary Fig. 9e, right). These results suggest that density is not only correlated with but also causal in shaping gastruloid patterning outcomes, and may serve as a target for rescuing disrupted developmental trajectories.

Modeling density-dependent morphogen signaling dynamics
Various studies have developed models of the effects of individual morphogens on cell fates.43-45 While combinatorial morphogen approaches have recently been applied to model cell fate decisions in non-patterned stem cells46 and 2D gastruloids,47 we lack a comprehensive model that incorporates combinatorial signaling dynamics and how those dynamics are decoded into cell fates. Such a model would enable an unbiased analysis of patterning outcomes in the gastruloid based on their fate patterning. Building off of previous PDE models of morphogen signaling,47 we utilized a reaction-diffusion description of BMP, Wnt, and Nodal dynamics (Fig. 4a, Methods). To model the fate choices, we incorporated equations that describe how morphogens are decoded into each of the three cell fate markers, GATA3, BRA, and SOX2 (Fig. 4b, Supplementary Fig. 10a, Methods). We note that existing models of cell fate in the 2D gastruloid incorporate entire signaling histories but only focus on a single morphogen,43 or incorporate multiple morphogens but neglect dynamics and make heuristic predictions based on instantaneous morphogen levels.47

Based on our observations of the importance of cell density, we sought to include its effects into our mechanistic model. Because hESCs form a tight monolayer, but also cannot grow beyond the bounds of the micropatterned ECM onto which they are plated, multiple morpho-mechanical features of cells in the gastruloid change as a function of density (e.g., aspect ratio, cortical tension, receptor availability, volume, etc.). We hypothesized that Wnt signal secretion, with its reported sensitivity to changes in local mechanics,48 is regulated by cell density. Mathematically, we captured this by introducing a function, h(d), that linearly decreases all sources of Wnt ligand creation as density, d, increases. See Methods for details of model implementation and justification of parameter choices.

To test if we accurately captured the effects of density, we simulated morphogen dynamics and cell fate determination of three colonies that range in densities, from 750 cells (~10th percentile) to 1150 cells (~90th percentile). Kymographs of simulated signaling revealed a faster Wnt wave that further permeates the colony interior as density decreases (Fig. 4c). As a result, Nodal, which is created as a function of Wnt, also follows this trend. Ultimately, this leads to a deeper penetration of BRA into the center of the colony with a corresponding decrease in SOX2 because BRA is a function of both Wnt and Nodal (Fig. 4d, Supplementary Movie 1). Experimental colonies with similar densities matched simulation results across a wide range of densities, even at the extremes of density and complete inhibition of a morphogen signal (Supplementary Fig. 10b-c). These findings suggest that our comprehensive model accurately captures the interactions between morphogen signaling, cell density, and cell fate markers.

A central implication of our model is that Wnt signaling is inversely proportional to cell density. To test this prediction, we sought to experimentally measure endogenous Wnt signaling in 2D gastruloids as a function of cell density. We leveraged a CRISPR-tagged β-Catenin (β-cat) hESC line49 (Methods) to observe Wnt signaling in real time. We varied initial seeding density and took live-cell, time-lapse images of 8 colonies per condition over the 48 hr differentiation (Fig. 4e). To extract the β-cat signaling activity we wrote a custom image processing algorithm that extracts the active, non-membrane β-cat from the inactive, membrane-bound pool, as has been done previously47 (Methods; Supplementary Fig. 10d). Both timelapse videos from individual colonies (Supplementary Fig. 10e-g, Supplementary Movie 2) and mean kymographs, subset by colony density, clearly demonstrate the impact of cell density on the Wnt wave (Fig. 4f). While the model captures density-dependent Wnt wave propagation, absolute signaling amplitudes differ from experiment, reflecting subtle differences in effective signaling timescales between simulation and experiment, as well as technical limitations in fully isolating active β-catenin. Finally, we fixed and stained the same colonies for which β-cat was imaged and confirmed that their BRA staining also decreased with increasing density (Fig. 4g). Overall, these experiments validate our choice of model architecture, confirm the effects of cell density on Wnt signaling, and suggest that our cell fate rules accurately capture the decoding process that occurs in the gastruloid.

Parameterizing gastruloid morphospace
We previously trained a neural network 𝜓(𝑦⃑) to approximate the original t-SNE embedding from the 150-dimensional phenotype vector. We now leverage this function to assign interpretability to morphospace by projecting simulated gastruloids with known developmental parameters into the same space. This allows us to understand how perturbations move colonies through morphospace, enabling mechanistic interpretation of drug action. Importantly, this analysis is independent of chemical structure and does not rely on the structure-to-phenotype model used elsewhere in the study.

To begin, we simulated gastruloids across a grid of two key parameters: cell density and SOX2 stability (Fig. 5a). SOX2 stability is a model parameter that reflects the resistance of SOX2 expression to mesodermal signals and captures how GATA and BRA modulate epiblast identity. Although not directly measurable, this parameter provides a mechanistically interpretable readout that broadly reflects drug mechanism of action and influences a colony’s position within morphospace. We projected all simulated colonies into morphospace using 𝜓(𝑦⃑) and found that these projected phenotypes formed a smooth grid spanning the canonical patterning region (Fig. 5b–e). The resulting structure reveals two interpretable axes aligned with these parameters, while not precluding additional contributing factors. To assess whether other coefficients could generate comparable trajectories, we systematically swept all model parameters over biologically plausible ranges; only density and SOX2 stability consistently produced extended, coherent axes in morphospace (Supplementary Figs. 11-12). Importantly, this organization persists across different t-SNE initializations, confirming that the parameter axes are not artifacts of a particular embedding (Supplementary Fig. 13).

We next asked whether this structure could be used to predict experimental outcomes directly from developmental parameters. To do so, we fit polynomials to the simulated grid to define density and SOX2 stability axes, then interpolated the SOX2 stability values for each drug-treated gastruloid from the location of its experimental centroid. We paired this interpolated value with the experimentally measured mean density of each drug treatment and used them as inputs to our morphogen model to simulate colony patterning from first principles. The predicted colonies were then projected into morphospace using 𝜓(𝑦⃑) and compared to the original experimental centroids (Fig. 5f). Notably, the simulated colonies closely matched both the morphospace locations and spatial phenotypes of the drug-treated gastruloids, validating the accuracy of the underlying model.

This approach demonstrates that the majority of phenotypic variance in the canonical patterning region can be described by just two parameters: cell density and SOX2 stability. By projecting simulated outcomes onto the same manifold as experimental data, we not only assign interpretable axes to morphospace, but also validate the predictive power of our mechanistic model, providing a framework to uncover how specific perturbations reshape developmental trajectories.

Linking black-box and mechanistic models in development
While our AI-based structure-to-phenotype model predicts 2D gastruloid phenotypes with high accuracy (Fig. 6a), it functions as a black box and does not directly reveal which developmental processes are perturbed. To provide mechanistic insight, we developed a complementary model grounded in morphogen diffusion and fate transitions (Fig. 6b). By projecting simulated phenotypes into the same learned manifold as experimental data, this model assigns interpretable structure to morphospace and enables causal reasoning about how perturbations modulate underlying parameters. In practice, the mechanistic model codifies the rules that govern the structured, canonical region of morphospace and, when supplied with prior mechanistic hypotheses (e.g., morphogen inhibition; Supplementary Fig. 10c), can also reproduce specific failure modes. The AI model, in turn, predicts where a compound will land in morphospace directly from chemical structure. Together, the two are complementary: the AI localizes compound effects, and the mechanistic model interprets whether they reflect canonical displacements or topology-altering failures.

Quantitatively, the black-box model (chemical structure → black box → phenotype) and the mechanistic model (initial conditions → PDE’s of morphogen spread and cell fate transition → phenotype) generated predictions with similar accuracies (Black-box model: mean RMSE = 0.17, mean cosine similarity = 0.90, 72% of predicted features within 1 standard deviation of replicate variation; Mechanistic model: mean RMSE = 0.18, mean cosine similarity = 0.93, 50% of predicted features within 1 standard deviation of replicate variation; Supplementary Data 3). However, the mechanistic model allowed us to predict which molecules were affecting Wnt signaling and cell density versus which molecules were affecting SOX2 stability.

Even deviations between mechanistic model predictions and experimental results revealed valuable insights. The morphogen spread model assumes fixed rules for BMP, Wnt, and Nodal signaling. However, certain compounds appear to override these assumptions. For example, isoxazole 9 (ISX-9)–treated colonies had relatively high cell counts (mean = 1,020 cells), which would predict moderate BRA expression in simulation; experimentally, however, they exhibited strong BRA expression (Supplementary Fig. 14a). This suggests an additional mechanism at play; indeed, isoxazole 9 has been shown to activate Wnt/β-catenin signaling by covalently binding Axin1 and stabilizing β-catenin,50 thereby enhancing mesodermal gene expression. Conversely, longdaysin-treated colonies had lower cell counts (mean = 892 cells), which would predict high BRA expression, but showed reduced BRA experimentally (Supplementary Fig. 14b), consistent with its reported Wnt-inhibitory action via disruption of LRP6 phosphorylation and reduced nuclear β-catenin accumulation.51 These discrepancies reveal when a perturbation alters the fundamental rules of morphogen signaling, enabling the model to identify biologically meaningful outliers rather than merely recapitulating the input.

These findings highlight an additional utility of the morphogen spread model as a tool for inferring mechanisms of action and identifying pathway-level effects of small molecules, offering a mechanistic complement to traditional phenotype-based screening in drug discovery.

Methods
Ethics Statement
All experiments involving human pluripotent stem cells were performed in accordance with applicable regulations and the principles outlined in the 2025 ISSCR Guidelines for Stem Cell Research and Clinical Translation.

Cell Culture and Differentiation
Cell lines
All experiments were conducted using H9 (WA09) human embryonic stem cells obtained from WiCell (Cat No. WB66446). For the initial drug screen, H9 cells were utilized to assess the impact of various drug treatments on patterning outcomes. Additionally, a CRISPR-Cas9 tagged β-Catenin (β-Cat) line was developed using H9 cells specifically for the quantification of Wnt signaling dynamics, enabling the study of active, non-membrane-bound β-Catenin in live cells.

Maintenance and culture conditions
For routine maintenance, H9 cells were cultured on Corning Matrigel hESC-Qualified Matrix, LDEV-free (Corning, Cat #L: 354277) coated dishes, and grown in mTeSR Plus feeder-free maintenance medium (STEMCELL Technologies, Cat #: 100-0276). Matrigel-coated dishes were prepared by coating them overnight at 4°C, followed by warming at room temperature for 1 hour before cell seeding.

Micropatterned gastruloid differentiation
For micropattern cell culture, custom ordered 500 µm diameter micro-patterned 96 well glass bottom dishes (CYTOO Inc, Cat #: A500P650) were first coated with CellAdhere Laminin 521 (STEMCELL Technologies, Cat #: 200-0117) at 10 ug/mL for 2 hours at room temperature. Established protocols for 2D gastruloid patterning were used. 28 Briefly, the wells are serially diluted with ice-cold calcium and magnesium free PBS (Thermo Fisher Scientific, Cat #: J67802.K2). Cells already resuspended in growth medium with 5 μM ROCK inhibitor (RI), Y-27632 (STEMCELL Technologies, Cat #: 72307) are then immediately plated upon removal of PBS. For cell seeding onto the micropattern, cells growing on Matrigel are lifted using ACCUTASE (STEMCELL Technologies, Cat #: 07920). Cells are centrifuged and 5x105 cells per well are resuspended in RI containing media and added to the micropatterned plate. After 4 hours the medium is replaced with mTeSR plus without ROCK inhibitor and incubated overnight.

Drug screening protocol
For the drug screen, we utilized the Stem Cell Signaling Compound Library (Selleck Chemicals, Cat #: L2100), containing 210 distinct drug conditions. The unannotated compound IB-2744 was obtained separately from Specs (Cat. #AK-968/15360521). Each compound was provided at an initial concentration of 10 mM, either dissolved in DMSO or water, based on solubility requirements. The compounds were stored and handled according to the manufacturer's recommendations to maintain stability and activity. To prepare the working solutions, the compounds were aliquoted into 96-well plates and diluted to a final concentration of 10 μM in mTeSR Plus medium supplemented with BMP4 (STEMCELL Technologies, Cat #: 78211.1) at 50 ng/mL and without RI. Following the initial preparation of the micropatterned colonies, as described in the Micropatterning section, the cells were incubated overnight in mTeSR Plus medium. The next day, the media was carefully aspirated, and the BMP4 + drug mixture was added to the micropatterned wells, exposing the cells to both the differentiation cue (BMP4) and the specific drug conditions simultaneously. The plates were then incubated at 37°C with 5% CO₂ for 48 hours, allowing the cells to respond to the combined BMP4 and drug treatment.

Wells are fixed with 4% paraformaldehyde and rinsed twice with PBS then permeabilized with -80 °C ethanol for two minutes. Wells are then incubated with primary antibodies (Dilution ratio: 300:1): anti-mouse GATA3 (Thermo Fisher Scientific, Cat #: MA1-028, RRID: AB_2536713), anti-goat Brachyury (R&D Systems, Cat #: AF2085, RRID: AB_2200235), and anti-rabbit SOX2 (Cell Signaling Technology, Cat #: 3579, RRID: AB_2195767), at 4 °C overnight, washed three times in PBS for 5 minutes each, and incubated with secondary antibodies (Dilution ratio: 1000:1) Alexa Fluor™ 568 Donkey anti-Mouse (Thermo Fisher Scientific, Cat #: A10037, RRID: AB_11180865), Alexa Fluor™ 647 Donkey anti-Goat (Thermo Fisher Scientific, Cat #: A-21447, RRID: AB_2535864), and Alexa Fluor™ 488 Donkey anti-Rabbit (Thermo Fisher Scientific, Cat #: A-21206, RRID: AB_2535792) and DAPI nuclear counterstain (Thermo Fisher Scientific, Cat #: D3571) for 1 hour at room temperature before being washed 3 times with PBS.

Genetic Engineering and Reporter Line Generation
CRISPR-Cas9 fluorescent tagging
For CRISPR-Cas9 editing in H9 human embryonic stem cells, we utilized Lipofectamine Stem Transfection Reagent (Thermo Fisher Scientific, Cat #: STEM00015), to transfect a plasmid expressing a gRNA targeting the N-terminus of the CTNNB1 gene, a homology template with ~180bp of both upstream and downstream homology to the targeted insertion site flanking the tdmRuby2 construct (Supplementary Data 5), and the Cas9 enzyme (IDT) according to manufacturer’s directions. Edited cells were passaged into a single cell suspension and sorted, into ROCKi (5 µM in mTeSR), for tdmRuby2 positive clones using a Sony SH800 FACS system. Clones we isolated and expanded prior to verification of successful edition through PCR followed by Sanger sequencing.

Imaging and Image Analysis
Imaging and quantification of gastruloid patterning
Fixed cell imaging of drug-treated gastruloids was carried out using a Nikon W2 SoRa spinning-disk confocal microscope equipped with incubation chamber maintaining cells at 37 °C and 5% CO2. We imaged in the four channels corresponding to DAPI, Alexa488, Alexa555 and Alexa647 conjugated antibodies. Images were exported as .tiff files and analyzed using the custom MATLAB and Python software, enabling quantification of the spatial distribution and intensity of cell fate markers across various drug treatments.

We developed custom MATLAB software to extract both raw and radial intensity profiles from imaged colonies. The pipeline begins by isolating the DAPI channel, converting it to a binary mask, and removing background noise using morphological filtering. A standard nuclei detection algorithm involving fill, dilation, and connected component analysis is then applied to generate a mask and identify nuclear positions. This is followed by secondary segmentation using enhanced binarization and a regional maxima filter to refine nuclei detection. Intensity values for each germ layer marker are extracted by applying the nuclei mask to the corresponding channels. Colony edges are determined by computing the centroid of all nuclei and applying a boundary detection algorithm to identify each cell’s proximity to the colony edge. The colony is divided into 50 concentric radial bins, and average marker intensities are computed per bin.

Live imaging and quantification of density-dependent mesoderm specification
Live-cell imaging of 2D gastruloids seeded at varying densities (10–100%) was performed using the CellVoyager CQ1 Benchtop High-Content Analysis System (Yokogawa Electric Corporation). Cells were seeded at 5 × 10⁵ per well in ROCK inhibitor-containing media, diluted to the desired concentration. Eight colony positions were selected and saved per well to enable consistent tracking across time. Imaging was conducted under incubation conditions (37°C, 5% CO₂), with z-stacks acquired every 15 minutes using brightfield and 561 nm laser illumination with autofocus enabled. Images were saved as .tiff files and analyzed using custom Python software to assess Wnt signaling dynamics.

Following live imaging, colonies were fixed, immunostained for BRA, and re-imaged using the previously saved positions to directly compare live signaling dynamics with endpoint fate outcomes. Fixed-cell imaging was conducted on the same CQ1 system, acquiring z-stacks under identical incubation and autofocus settings. Images were analyzed using custom Python and MATLAB software to quantify germ layer marker expression across seeding densities.

Wnt signaling quantification and kymograph generation
To quantify active Wnt signaling, we developed a custom Python-based image processing pipeline to isolate non-membrane (active) β-catenin from the membrane-bound pool. After background subtraction on .tiff images, a Gaussian blur reduced noise, and adaptive thresholding generated a binary membrane mask. This mask was dilated and subtracted from the original image to isolate cytoplasmic and nuclear β-catenin. Morphological operations further refined the mask, and mean fluorescence intensity was calculated to quantify Wnt signaling activity. To assess temporal dynamics, the same pipeline was applied to time-lapse images, with each colony divided into 50 radial bins. Mean intensities were computed per bin across frames to generate kymographs capturing spatial and temporal patterns of β-catenin signaling.

Zebrafish Teratogenicity Assay
Zebrafish maintenance and compound treatment
All experiments were conducted using wild-type AB zebrafish (Danio rerio) maintained at 28 ± 1 °C with a 14-hour light and 10-hour dark cycle. All zebrafish husbandry, breeding, and experimental procedures were conducted by a commercial research organization (ZeClinics, Barcelona, Spain) in accordance with institutional guidelines and applicable national and European regulations governing the use of animals in research. Fertilized zebrafish embryos were collected in E3 1X medium (60X stock is prepared by dissolving 34.8 g NaCl, 1.6 g KCl, 5.8 g CaCl ·2H O, 9.78 g MgCl ·6H O in 2L H O and adjusting the pH to 7.2 using NaOH) in Petri dishes, and abnormal or unfertilized embryos were discarded at 6 hours post-fertilization (hpf). Healthy embryos were grown up to 96 hpf and exposed to ulixertinib, nirogacestat, and isotretinoin at five concentrations (0.01 µM, 0.1 µM, 1 µM, 10 µM, 100 µM) starting at 6 hpf. 4-Diethylaminobenzaldehyde (DEAB) (MilliporeSigma, Cat #: D86256-100G) was used as a positive control at concentrations of 0.1 µM, 0.3 µM, 1 µM, 3 µM, and 10 µM.

Phenotype quantification and analysis
Mortality was assessed at 24 h and 96 h, and larvae were imaged at 96 hpf using the automated VAST system (Union Biometrica). Dead, unhatched, and incorrectly detected larvae were excluded from analysis, resulting in a variable number of larvae analyzed per condition. Teratogenic effects were quantified by assessing both qualitative and quantitative phenotypes. Qualitative phenotypes included body curvature abnormalities, yolk edema, necrosis, and craniofacial defects, while quantitative phenotypes such as lateral area, eye diameter, pigmentation, and body length were measured using FIJI software. Continuous data were transformed to binary outcomes based on interquartile range thresholds.

LC50 and EC50 values, representing the concentration at which 50% of mortality is observed and the concentration at which 50% larvae population show a teratogenic phenotype, respectively, were calculated using dose-response models via Proast software, allowing for the determination of the teratogenic index (TI) as the ratio between LC50 and EC50 for the most sensitive phenotype. The Lowest Observed Effect Concentration (LOEC) was identified for each compound, with effects within 20% of biological range considered non-significant.

Data Analysis and Computational Modeling
Dimensionality reduction and clustering
Bin-wise features were normalized per plate using robust scaling (25th–75th percentiles) to mitigate batch effects and outliers. The normalized single-colony profiles were embedded into two dimensions with t-SNE. To define morphological phenotypes, we applied a Gaussian kernel density estimate (KDE) over the t-SNE space of all experimental colonies, identified local density maxima, and then performed watershed segmentation to partition the morphospace into discrete, contiguous regions (Supplementary Fig. 1g). This unsupervised procedure does not require pre-specifying the number of clusters and yields cluster boundaries that reflect the empirical structure of the embedded data. While alternative clustering methods (e.g., Spectral or Leiden) produce similar numbers of clusters, the watershed-defined boundaries offer a key practical advantage in that they can be directly reused to assign new predicted phenotypes without re-running clustering, ensuring consistent and unbiased classification across experimental and in silico data.

Simulation of reaction-diffusion models
The reaction-diffusion models were simulated using a custom Python script to capture the dynamic interactions of signaling pathways in the system. The simulation involved creating a 2D spatial grid to represent the environment and updating the concentrations of key signaling molecules over time. Standard numerical techniques were applied to compute diffusion, reaction, and degradation processes. Custom functions handled boundary conditions and interactions between signaling molecules.

Simulation of cell fate specification
Cell fate simulations were conducted using custom Python scripts designed to integrate the dynamics of morphogen signaling with rules for determining cell fate outcomes. The process involved calculating fate maps based on our cell fate equations to determine the expression levels of GATA3, BRA, and SOX2 proteins. Visualization of fate maps was performed by generating heatmaps of normalized expression levels using custom color maps to represent each fate (GATA, BRA, and SOX). Separate images were created for each morphogen, and a composite image was generated to illustrate the combined spatial distribution of all fates within the colony. Additionally, radial intensity profiles were computed by binning the colony into concentric rings and averaging the protein expression levels within each ring.

Neural Network-Based Modeling
Chemical structure-to-phenotype prediction (ChemBERTa)
To predict gastruloid phenotypes directly from chemical structure, we used a two-stage model based on molecular embeddings generated from ChemBERTa. ChemBERTa was downloaded from the HuggingFace repository seyonec/PubChem10M_SMILES_BPE_450k, which was pretrained on over 10 million SMILES strings using masked language modeling and byte pair encoding (BPE). SMILES strings were tokenized and embedded with the HuggingFace Transformers library, yielding a fixed-length molecular representation for each compound. The dataset was split 90/10 into training and blind test sets, with each compound represented by a ChemBERTa embedding paired to its corresponding radial phenotype vector.

In the first stage, we trained a neural network classifier to map ChemBERTa embeddings to logits over morphological clusters. The classifier consisted of two fully connected layers (512 and 256 units) with ReLU activations and dropout, followed by a linear output layer. Training used the Adam optimizer (learning rate = 0.001) with a combined cross-entropy and cost-sensitive loss designed to penalize misclassification of failure-mode clusters. Optimization was performed using 5-fold cross-validation with an early-stopping criterion based on validation accuracy (patience = 30 epochs). Final cluster assignments were obtained from fold-averaged logits. Each compound was then assigned the mean 150-dimensional phenotype vector of its predicted cluster.

In the second stage, we modeled structure-dependent deviations from these cluster means. For each cluster, a ridge regression model was trained to predict the residual between the cluster mean and the observed 150-dimensional phenotype vector. The final output phenotype vector for each compound was obtained by adding the predicted residual to the cluster mean, and accuracy was evaluated by comparing predicted and experimental radial profiles and morphospace locations.

Approximation of t-SNE embedding
To approximate the t-SNE embedding of high-dimensional phenotypic data, we trained a neural network, 𝜓(𝑦⃑), to predict low-dimensional projections from bin-wise fate marker distributions. The network comprised three fully connected layers with 128, 64, and 2 neurons, respectively, and ReLU activations. Input data were split 80/20 into training and validation sets.

Training was performed for up to 1,000 epochs using the Adam optimizer (learning rate 0.001) with mean squared error loss. A learning rate scheduler was used to reduce the step size when validation loss plateaued, and early stopping (patience = 50 epochs) prevented overfitting. Model accuracy was assessed by comparing predicted versus true t-SNE embeddings using Jensen-Shannon Divergence (JSD) and visual overlay of density plots.

Mathematical Modeling
Reaction-diffusion model for morphogen spread
Inspired by recent studies examining BMP-Wnt-Nodal dynamics in 2D gastruloids,47 we developed a simplified reaction-diffusion model tailored to our experimental observations. The model specifically addresses the influence of cell density on mesoderm patterning, as we observed that lower cell density shifts the mesoderm peak inward and causes the mesoderm region to expand outward. This density-dependent modulation was incorporated into our model, accounting for the dynamic interactions between BMP, Wnt, and Nodal, and aligning with observed spatial and temporal morphogen patterns.

BMP was modeled with the following equation:
𝑏(𝑥,𝑦,𝑡) = 1 , 𝑡 ≤ 𝑡
𝑙𝑜𝑐
𝜕𝑏
{ = 𝐷 ∆𝑏 , 𝑡 > 𝑡 , |(𝑥,𝑦)| ≤ 𝑟 (1)
𝑏 𝑙𝑜𝑐 𝑙𝑜𝑐
𝜕𝑡
𝑏(𝑥,𝑦,𝑡) = 1 , 𝑡 > 𝑡 , |(𝑥,𝑦)| ≤ 𝑟
𝑙𝑜𝑐 𝑙𝑜𝑐
The BMP signaling dynamics are modeled as a diffusive species within a localized region, r , which varies based on cell density. Initially, BMP concentration is set to 1 across the entire domain, representing the uniform activation phase. After t , BMP diffusion is restricted within r , simulating the feedback inhibition observed experimentally, where BMP signaling diminishes centrally but remains active at the periphery.40 In this context, D , represents the diffusion coefficient of BMP, r , denotes the radius defining the localized BMP diffusion modulated by cell density, and t indicates the time point when BMP signaling becomes localized. The boundary for BMP signaling, r , is defined by 𝑟
𝑟𝑎𝑑
∗(1−𝑒−𝑎∗𝑑), where and r
rad
is the radius of the stem cell colony. Here, a is a constant that scales the effect of cell density, and d represents the cell density; this formulation reflects our observation that the outermost cells spread more at lower densities, expanding the peripheral region in which BMP signaling remains active. The spread of Wnt was modeled by:
𝜕𝑢
= 𝐷 ∆𝑢+ℎ(𝑑)(𝑓(𝑢)+𝑘 𝑏+𝑐 )−𝑘 𝑢 (2)
𝑢 E 𝑏 𝑢 𝑢
𝜕𝑡
Where:
𝑑
ℎ(𝑑) = 1−𝑐 ( −1) (3)
𝑚
𝑐𝑑
𝑓(𝑢) =
𝑠𝑢𝑢2
(4)
1+𝐾𝑢𝑢4
The Wnt signaling dynamics incorporate diffusion and autocatalytic activation, modulated by the density-dependent scaling function h(d). The term f(u) captures the autocatalytic production of Wnt, with inhibition effects included through a saturation term in the denominator to limit excessive Wnt production at high concentrations. Here, D , is the diffusion coefficient of Wnt, h(d) is the density-dependent scaling function, s denotes the rate of Wnt autocatalytic production, K is the saturation constant controlling Wnt production, k represents the contribution of BMP to Wnt production, c is a constant source term for Wnt, and k is the degradation rate on Wnt. The scaling function h(d) modulates Wnt reaction terms based on cell density, where c is a scaling factor, d is the cell density, and c is a constant representing the mean cell density from the screen. The spread of Nodal was modeled by:
𝜕𝑣
= 𝐷 ∆𝑣 +𝑔(𝑢,𝑣)−𝑘 𝑣 +𝑐 (5)
𝑣 𝑣 𝑣
𝜕𝑡
Where:
𝑔(𝑢,𝑣) =
𝑠𝑣𝑣
𝑢 (6)
1+𝐾𝑣𝑣2
The Nodal signaling dynamics incorporate diffusion and production influenced by Wnt activity, as captured by the term g(u,v), which includes a saturation effect to reduce Nodal production at high concentrations. Unlike previous models, we chose not to model Nodal as autocatalytic due to evidence indicating that Nodal spreads primarily through a relay mechanism rather than an autocatalytic loop.80 In this context, D , denotes the diffusion coefficient of Nodal, s , is the rate of Nodal production influenced by Wnt, K is the saturation constant reducing Nodal production at high concentrations, c represents the constant source term for Nodal, and k is the degradation rate of Nodal. All parameter values used in the model can be found in Supplementary Data 6.

Morphogen-based cell fate determination model
We developed cell-fate determination rules based on the integration of BMP, Wnt, and Nodal signaling dynamics. Each fate map was generated by simulating the morphogen concentrations over time, followed by applying specific computational rules to map the spatial location of fate-markers, GATA3, BRA, and SOX2, throughout the colony. Specifically, GATA3 expression was modeled by:
∫𝑏(𝑥,𝑦,𝑡)𝑑𝑡
𝑮𝑨𝑻𝑨𝟑 = (7)
1+𝑒−𝑘1(√𝑥2+𝑦2−𝑟𝑙𝑜𝑐)
The expression of GATA3 is modeled as a response to the integral of BMP signaling over time, reflecting the cumulative impact of BMP exposure on cells.43 The denominator accounts for the preferential availability of transforming growth factor β receptors at the edge of the colony,40 with the steepness parameter modulating the influence of distance from the colony edge. Here, √𝑥2+𝑦2 represents the radial position within the colony, r defines the localized region where BMP signaling is most concentrated, and k adjusts the drop-off rate of GATA3 expression as distance increases from the colony edge increased. BRA expression was modeled by:
1 𝜕𝑣(𝑥,𝑦,𝑡)
𝑩𝑹𝑨 = 𝑑𝑡 (8)
1+𝑒−𝑘2(𝑢(𝑥,𝑦,𝑡)−𝑊𝑛𝑡𝑡ℎ𝑟𝑒𝑠ℎ 𝜕𝑡
BRA expression is governed by the combined influence of Wnt and the rate of change of Nodal signaling, capturing the dynamics where Nodal alone does not induce mesoderm formation.81 Here, the Wnt concentration must reach a certain threshold, Wnt in order for the cells to respond to changes in the rate of Nodal concentration. Specifically, unlike BMP, which is concentration-dependent, our experimental observations indicate that many Nodal targets depend on the rate of concentration change rather than its absolute levels. Here, k represents the steepness of the BRA activation response to Wnt, Wnt denotes the threshold level of Wnt required for BRA expression, and is the rate of change of Nodal concentration over time. SOX2 expression was modeled by:
𝑺𝑶𝑿𝟐 = 𝐶 (𝑐 +𝑐 𝑮𝑨𝑻𝑨𝟑+𝑐 𝑩𝑹𝑨) (9)
𝑠 1 2 3
SOX2 expression is modeled as inversely proportional to the levels of GATA3 and BRA, aligning with our observations that SOX2 expression typically decreases as cells begin to differentiate. The stability constant, C , modulates the persistence of SOX2 under varying conditions, reflecting our experimental findings that certain drug treatments allow cells to maintain high SOX2 expression despite concurrent BRA expression. In this equation, C , C , and C are constants that adjust the relative influence of GATA3 and BRA on SOX2 suppression. All parameter values used in the model can be found in Supplementary Data 7.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this work, the authors undertake a formidable effort to screen and quantify the effects of 210 compounds on the development and differentiation of 2D gastruloids from hPSCs, a system well established, primarily by the Warmflash lab. They quantify three fluorescent markers in 10 colonies per treatment, ultimately quantifying images from over 2,000 gastruloids at the single-cell level, a remarkable scale for this kind of analysis. The fluorescence quantification is converted to a vector describing the radial intensity profile for each marker.\n\nThe idea of generating a mean behaviour image for each cluster is nice, seems to work well in these radially symmetric colonies, and provides an intuitive understanding of each response type. The authors cluster the gastruloid vectors into 12 clusters, assigning 9 or them to normal behaviour and 3 clusters to distinct failure modes (loss of Bra/Sox2 or loss of radial symmetry). Note that converting a profile with a geometric meaning (expression along the radial axis) into a vector (simply by binning) looses all geometric information. For example, under standard metrics the vector 1000 will be equidistant from 01000 and from 00001, while in its geometric interpretation it is clearly closer to the first one. All the same, the clusters obtained seem to be meaningful. The authors may want to discuss this point, and perhaps point to limitations of this embedding. Of course, when we move to network learning in the next step, the network model can re-learn the geometric connections between output nodes.\n\nIn the next step, the authors encode the 210 chemical compounds as vectors, and train a network model to predict the gastruloid vector (output) from this encoded input. They also train an additional (smaller) network to embed the predicted vector in tSNE “morphospace embedding”. I assume the prediction errors were done on the original 150-dim vector space, correct? Why is it important to embed into the 2-dim tSNE space? It is only used for visualization purposes. As the name says, it is only a 2-dim embedding of a 150-dim vector, which itself is a representation of a 2-d circular pattern. If the purpose is mostly to evaluate the error (as depicted in Fig. 2b) for the 150-dim space, the embedding process only adds some more distortion to it. It is more informative to show some quantification of the prediction errors (as depicted in the red/green/blue plots), how they behave in training vs. test set, etc.\n\nIn the sample images in Fig. 3a we see some bald regions within the colony (which we do not see in the examples in Fig. 2. What are these? Why do we only see them here?\n\nThe authors then run 947 FDA approved drugs through their prediction pipeline, finding 100 to fall within the no-Brachyury cluster. The association statistics corresponds with pre-known risk levels (as teratogens) of those drugs. A few comments on this approach/result:\n- Looking at Fig. 3c, it is not clear how association of the new compounds to cluster 11 was determined. Some of them (including several category X compounds) seem just as close (or closer) to clusters 8 and 9. Or is it because the association was computed in the 150-d space, and the 2-d projection is just misleading? The same is true for Ulixtertinib prediction.\n- In this specific example, the association to cluster 11 describes failure in the form of no Brachyury expression. Did we have to go through the sophisticated experimental model (2D gastruloid) and heavy computational framework? Wouldn’t we get a similar (or better) prediction simply screening Bra-expression levels (i.e. a 1d readout instead of 150d readout) in e.g. simple diff protocol applied to a 2D culture? It is worth at least discussing this point.\n\nThe authors then devise a mechanistic model, that models the spatio-temporal dynamics of Wnt, BMP and Nodal activity, and tie it with the measureable marker expressions (that serve as a proxy for cell fate) through a simple logic-based model. They demonstrate how such a model can predict the spatial patterns observed, as well as their dependence on cell density. They experimentally validate one implication from the model simulations, which is the inverse dependence of Wnt signaling on cell density.\n\nThey then parametrize morphospace by projecting gastruloids predicted by the mechanistic model while varying two parameters: cell density and SOX2 stability. It is not clear (or well motivated) why SOX2 stability was chosen – this is not a parameter that can be directly controlled experimentally, and it is not clear to what degree it naturally varies between the gastruloids. The authors claim that Sox2 stability (Cs) “emerged as a key regulator of fate boundary positions in our model.“ This claim is a bit misleading – the model was devised in a way that gives Cs a high weight, determining the resulting patterns. For example, the stabilities of Gata3 or Brachyury do not appear is parameters in the model at all. This should be better explained and motivated.\n\nThe idea of parametrizing the abstract UMAP embedding by varying interpretable parameters is nice, as it adds to the interpretability of the morphospace.\n\nThe authors present the AI and mechanistic models as if they were two alternatives to choose from. This is a bit misldeading, as the mechanistic model does not incorporate the potential effects of teratogens. If I understand the rationale of the authors, in integrating (and comparing) the two, they assume that the effect of the teratogens can be summed up through their effect on just 2 parameters (cell density and SOX2 stability). What is the rationale for that? Are the effects of the compounds on density simply due to toxicity? (this could be cross-checked against simple growth curves in 2D cell culture for these compounds). And if they affect Sox2 stability (as rationalized in the Discussion), don’t they have similar effects on GATA3 and Brachyury? Or on Wnt, Nodal and BMP proteins?\n\nIn summary, the effort to generate a prediction system for the effects of compounds on human embryonic development is commendable. The AI-generated predictions seem to generate useful predictions. The integration with the mechanistic model, while adding mathematical gravity to the manuscript, is not very clear currently.", "reviewer_2": "Summary:\nThis paper develops a framework for high throughput drug screening and phenotype prediction for 2D gastruloids. Improving on previous work describing teratogenicity screening in vitro, the authors develop an integrated approach of deep learning and mechanistic modeling to predict and interpret phenotypes. The paper contains several nice contributions to the field, including 1) the machine learning model mapping phenotypes into a (somewhat) interpretable 2D space, 2) a large dataset made public of the effect of different drug treatments on 2D gastruloids, 3) an extension of the model from Chhabra et al 2019 to include cell density and explicit modeling of cell fate. However, a major weakness is a lack of rigor and detail, including poorly motivated choices in approach. For a paper about what is essentially a statistical testing and prediction framework, rigor seems particularly important.\n\nMajor comments:\n- The clarity and accuracy of the abstract could be improved. First, it claims to elucidate “failure modes” but neither the abstract nor the rest of the paper clearly defines what a failure “mode” is. If what is meant by this is that the phenotypes can be captured by cell density and SOX2 stability, perhaps that should be made more explicitly. It also says the paper “revealed” cell density and SOX2 stability as determinants of patterning, but both were already known, the author do in fact cite previous work showing the importance of density (e.g. Etoc et al) and SOX2 stability (Camacho-Aguilar et al, Teague et al). Perhaps corroborate is a better word and the discussion would ideally more explicitly incorporate this past work.\n- Rigor should be improved:\no The paper states “10 colonies per drug condition were imaged”. However, it is not clear whether these all came from the same experiment and if so if they came from the same well. Inter-colony variation within a well is very small, between wells on the same plate small, between experiments or different plates within an experiment can be significant. To understand if phenotypes of two drugs are significantly different, or whether the prediction is correct within the experimental uncertainty, it is essential to quantify the inter-experiment phenotypic variation under the same drug perturbations and to make clear whether all drug treatment data was collected in the same experiment and if not to show if there are differences in the control condition between experiments and determine if there are any batch effects.\no The untreated (BMP only) control conditions is not shown in most figures but should be for comparison. For example, in which cluster does if fall in Fig 1g and where is it in 1f, and in Fig 6c – how does Raloxifene HCl compare to untreated? Can it be said to be significantly different? The same for all other figures. If I understand it correctly, Supplemental Fig. 5a shows the BMP only control clusters and their natural variation. This covers a large part of the map and the majority of the ”Canonical Cell Patterning” clusters in Fig.1g, suggesting that those clusters should perhaps be a single cluster and that most drug treatments within that cluster cannot be distinguished from control.\no More generally, it is not clear the clusters are meaningful. The authors use a custom method based on watershed segmentation. It is not clear why they need this, and the method is not described in sufficient detail to be reproducible. To motivate the approach, it is stated that an advantage is that the number of clusters is not predetermined but there are many methods for which that is true, such as Leiden/Louvain. Furthermore, one can optimize cluster number for methods with fixed cluster number using cluster analysis. Regardless, cluster analysis should be performed here to determine if the obtained clusters are statistically meaningful, e.g. calculating Silhouette coefficients, etc (https://en.wikipedia.org/wiki/Cluster_analysis). This reviewer would not be surprised if the number of statistically and biologically meaningful clusters is a lot smaller than the 12 that are discussed.\no Error bars on radial profiles. The analysis throughout the paper only considers the mean radial profile for each colony, but there are error bars on these radial profiles that may also change by perturbation, e.g. the variance in expression for some distance from the edge goes up when the patterns become less symmetric, and generally one could imagine propagating these error bars to the morphospace to get a sense of how significantly different the profiles represented by dots there are.\no In testing the prediction framework, the separation into training and testing data is essential. How was the 10% of compounds held out for test selected? Randomly? It would be good to see N-fold cross-validation to see how much performance varies as different compounds are held out for validation.\no The authors use RMSE as a measure of goodness of fit, but these numbers have no absolute meaning. Cosine similarity is also hard to interpret as good or bad. More meaningful would be to see if prediction and data agree within the experimental error.\no In figure 3, the authors compare their prediction to a “ground truth” drug classification. They should show at least a confusion matrix with accuracy, % false positive, false negative, perhaps also check other measures like precision. Systematically showing %enrichment of drug classes in clusters would also be good. Looking at Figure 3b it looks well mixed.\n- For figure 3, the authors focus on ulixertinib and nirogacestat for verification. Why are these not shown among the data points in Figure 2? Looking at the Supp. Fig. 4ab, the colony images show dark blue lines streaking along and across the colonies. These do not look like nuclear stains for the expected markers. What are they? If they are not cells, how do they affect analysis? are they excluded? If they are cells, why does their pattern look so strange?\n- As stated above, the modeling in Fig. 4 is one of my favorite parts of this paper, but it does raise some questions that if addressed would further strengthen the manuscript. The authors model Wnt production as density-dependent and verify the predicted Wnt signaling dynamics. This nicely captures the narrower range of Wnt signaling at higher density, but their model also shows a later onset of Wnt signaling in Fig.4c which the experimental data in Fig.4e does not seem to show. Can this be explained? Can the model be improved to capture this? If it cannot be addressed, perhaps the authors can add it to the discussion.\n- As a separate point about the model, it is surprising that SOX2 is modeled essentially as “not GATA3 or BRA”, given that causally it is not true that SOX2 is directly downregulated by these factors (possibly the reverse) and that ultimately the authors emphasize SOX2 stability as the most important model feature. Is it fair to say SOX2 stability simply means degree of differentiation the way the model is set up? More importantly, if the authors truly think drugs directly affect SOX2 stability, it would be good if they would show this, similar to how they showed the density-dependent Wnt signaling.\n- The discussion illustrates the issue with thinking of all drugs as directly affecting SOX2 stability. It mentions in line 620 that ROCK inhibitors destabilize SOX2. However, what more likely happens is that they prevent tight junction formation, leading to increased BMP signaling, which is what actually destabilizes SOX2. In other words, there are many possible indirect of a drug that can affect differentiation.\n- In lines 487-496, the authors ask if the mechanistic model can capture the phenotype of drug treatment. They infer the SOX2 stability from interpolation of simulation for different SOX2 stabilities on the morphospace. Thus, it is not surprising that if one provides SOX2 stability and density one ends up roughly in the same place, the only way I can see this going wrong is if the experimental density used deviates strongly from the mean density coordinate on the morphospace. I guess this is what the authors may mean in lines 548-553. To summarize, this seems somewhat like circular reasoning and the authors may wish to make clearer whether this test really proves anything. Perhaps this could involve merging figures 5 and 6. Can the authors discuss more why the model cannot capture the effect of isoxazole and longdaysin? Should it not predict reduced/increased SOX2 stability for a Wnt agonist/inhibor?\n\nMinor comments:\n- The introduction line 60 states: “methods to assess the teratogenicity of compounds have not advanced beyond exposing model organisms”, but this is inaccurate, and the authors seem aware as they later cite two micropatterned approaches of teratogenicity testing. There are several other papers describing vitro approaches they may also wish to include,.e.g., PMID: 26387793, PMID: 31207406, and others.\n- Line 81: “3D models such as embryoids, blastoids, and SEMs” – embryoid is a general term that includes blastoids and SEMs. Whas something different meant?\n- Line 85: “cost of altering the topology of the developmental tissue” – which tissue? It is true that the topology of the amnion is altered, but the epiblast geometry and topology is a disc as in vivo, which it actually is not some 3D models. It is a misunderstanding that 3D is always more similar to in vivo than 2D.\n- Figure 1c shows classes of signal pattern perturbation but as far as I can tell two of these never occur: inverted and mixed, so why list them?\n- Clusters should be indicated in Fig.1e.\n- The authors use GELU for one network, RELU for another and make a point out of this in the main manuscript, however, it is not explained why. More generally data supporting the choice of network architecture is missing.\n- The shadows underneath the colony pictures throughout the manuscript, e.g., Fig.1g make the patterns harder to distinguish, perhaps a black outline would be better if contrast from the which background is needed.\n- Line 345: “cell seeding follows a Poisson distribution”. That is only true for very sparse seeding and typically not for micropatterns which are densely seeded. The authors should show this or could simply remove it since it does not really matter, their only point is there is naturally occurring density variation.\n- The regression of BRA peak vs. cell number in Supp. Fig. 5ac is quite nice. However, it would be good if in Fig. 1c not only the cluster averages but all datapoints were shown on the graph. Conversely, it would be helpful to see the cluster averages on the t-SNE plot.\n- In line 461-466, it is said density and SOX2 stability emerged as key parameters. How did they emerge? How were key parameters identified?\n- Line 566: “animal models, which are limited by species-specific differences and dose-dependent lethality”. This appears to suggest that dose-dependent effects are not an issue in this work but they are and the authors acknowledge this later in the discussion (everything done at 10uM), which is good, so they could consider rephrasing line 566.", "reviewer_3": "Summary:\nIn this manuscript, Rufo et al., develop an experimental and analytical screening platform, FATE-MAP, to assess phenotypic outcomes in human 2D stem cell models of early embryogenesis (human 2D gastruloids) treated with hundreds of chemical perturbations. FATE-MAP includes a deep learning framework to predict the effects of any molecule on 2D gastruloids development based on the molecule’s chemical structures. The authors first characterized “failure modes”, phenotypes that emerge when the 2D gastruloids are not properly patterned and then analyze over two thousand structures subjected to roughly 200 perturbations, including clinically relevant teratogens, eventually identifying several molecules that results in aberrant developmental programs. The authors also performed some validation experiments using the zebrafish model to corroborate the predicted molecule-induced phenotypic outcomes. They then model morphogen dynamics in 2D gastruloids, an approach that revealed how cell density and expression stability of the transcription factor SOX2 play a major role in the correct patterning of 2D gastruloids. Finally, they attempt to derive mechanistic principles of failure modes from their deep learning approach, implementing a mechanistic model to gain interpretability and biological meaning.\n\nOverall, this reviewer thinks the developed methodology is important and timely, given the little knowledge we have about the effect of environmental offences on embryonic development, especially humans. Moreover, cheap and accessible high-throughput screening opportunities are being developed, and I think this method and its future implementations (like others of its kind) are in need to the field. Science and experiments are solid and well executed. However, this reviewer thinks the biological findings are not groundbreaking, especially given the opportunity provided by the method. The manuscript is very well written and easy to understand for a broad audience. I generally support its publication at Nature Communications, pending some points that would need to be addressed.\n\nMajor comments:\n1) List of selected molecules:\n- It would be important if the author could comment more on the list of 210 drugs tested. How were they selected? What is known about their effect in other species? Are they already classified or stratified based on their chemical structures?\n- This becomes more relevant in the efforts to build a classifier based on chemical structures. Are all molecules causing one phenotypic cluster similar/dissimilar in structure? Do they share any chemical group/property (e.g. structure, domains, etc..)? Figure S3a looks at the similarity between seen and unseen drugs in the classifier but does not compare the general chemical structures in the drug list.\n- What is the inter cluster vs intra cluster chemical structure similarity? Do the tested molecules cluster based on their chemical structure?\n2) Screening design:\n- The authors should test different concentrations for the drugs they use. This reviewer understands well that this won’t be possible for all compounds, but using a fixed concentration for all drugs (10uM) limits the real assessment of their effects given that not every molecule has the same stability, processivity, and activity. The authors should titer different concentrations to understand dose-specific responses for a shortlisted number of compounds (e.g. the ones resulting in phenotypes belonging to C1, C11, and/or C12).\n3) Comparison with other screening approaches, deep learning methods and relevant literature:\n- A similar teratogen screen has been performed in 3D human gastruloids, which offer a more physiological and morphologically relevant model (PMID: 34425190). This work uses a way shorter list of compounds but various drug concentrations and combines morphometric analysis and gene expression patterning. Given that many molecules are present in both studies, the authors should compare their observed phenotypes in 2D in relation to the one in 3D. Also, if some molecules that have been used in PMID: 34425190 are not present in the authors’ list, they should predict the phenotype using FATE-MAP and compare the prediction with the observed phenotype in 3D.\n- Another relevant study has been recently published (PMID: 40245869) combining morphological variation in 3D mouse trunk-like-structures with gene expression profiles to train a classifier able to predict phenotypic outcomes from early morphological features. Even if this reviewer acknowledges that the experimental design is different (this work uses longitudinal imaging over time to identify predictive features), the authors should cite this work and comment on similarities and differences in their approaches.\n- A recent study in engineered mouse 3D gastruloids (PMID: 39358450) provides temporal recording measurements of signaling patterns and specifically studies the emergence of Wnt domains and its relationship to Nodal. The authors should comment and relate their results to the ones reported in this work, since this comparison could provide species-specific and shared features between mouse and human.\n4) Use of the term “morphospace”:\n- The authors use the term “morphospace” to define the landscape of phenotypic outcomes arising in the treated 2D micropatterns. Even if several colleagues in the field have argued about the morphological organization of micropatterns, including the 2.5D nature of the model with a z dimension also being relevant for proper cellular organization (which the authors do not use in their analysis), this reviewer thinks that 2D gastruloids lack morphological features of an embryo. The model is great to study tissue patterning, cell interactions and the effect of molecules on these features, but not so great to assess the acquisition of complex embryonic morphologies. I would recommend converting the wording “morphospace” into a “patterning space” (or similar) since this is what the model (and the readout used by the authors) can accurately inform.\n5) Predictive power:\n- The screening experiments performed by the authors offer the unique opportunity to inspect what was not “well predictive” (e.g inter cluster violations between prediction and observation) which could potentially result in unexpected, interesting observations/exceptions to follow up with. Did the authors observe any of these violations in their data?\n- In line with what this reviewer commented on before, are there any features of the molecules’ chemical structures with predictive power? Finding chemical similarities that results in similar phenotypic outcomes might be a key and novel aspect of the study.\n6) Cell death vs differentiation defects\n- Cell death measurements are not included in the experiments. For example, C11 (failure to express BRA) could result from cells failing to differentiate (as hypothesized by the authors) or compound lethality (which could still result in a similar phenotype). To pinpoint a real example, there is a steep difference between Cerdulatinib treatment (looks more like a lethal phenotype) and D4476, but both end up in C11. The authors should classify C11 further and pinpoint lethality vs differentiation failure phenotypes. Ideally, for a set of compounds, a live/dead dye should be used to determine the cause of the observed failure mode. This is extremely relevant to correctly interpret the phenotypic outcome and failure mode.\n- Cell density vs BRA expression: this relationship is established at the endpoint differentiation. It is known that exit from pluripotency to the PS is accompanied by cell death (and therefore might result in lower density). Is this hypothesis tested against the initial cell density, or the observed density reduction could be a consequence of enhanced differentiation?\n7) FDA drugs effect prediction:\n- This reviewer understands well the challenges associated with implementing the Pregnancy and Lactation Labeling Rule (PLLR) for FDA approved drugs given their non categorical classification, but the categorization used by the authors dates back from 1970s and lack rigorous testing. As specified by the authority and clinicians, this categorization is not a “grading system of risk” based on the actual severity of the teratogen, but rather a classification based on available data at the time. The authors should compare/consider including information from PLLR in their analysis. This could even improve their results since the current classification might suffer from lack of extensive testing (see the next two points regarding the A, B and X graded molecules).\n- On a similar line, the authors should show the embedding (or the statistical distance to clusters) for all FDA drugs, including class A and B compounds (not shown in the current version).\n- Where are all the other class X compounds ending up (71 out of 98 class X compounds are not shown). Are they part of the other aberrant clusters or do they end up in the canonical “morphospace”?\n- The authors should also comment on the known phenotypes (at the cellular and organismal levels) associated with Ulixertinib and Nirogacestat in other cellular systems or animals. Do their predicted and observed phenotypes relate to the literature knowledge of these molecules and to their molecular functions (ERK 1/2 and γ-secretase inibitors, respectively)?\n8) Zebrafish experiments:\n- This is, according to this reviewer, the most important technical aspect the authors should clarify. Based on both the observed and predicted phenotypes, both Ulixertinib and Nirogacestat treatments result in phenotypic cluster C11, where micropatterns fail to undergo gastrulation and to induce the expression of BRA, consequently lacking the primitive streak (PS). The timing of the treatment in zebrafish embryos is correct, matching the epiblast stage in human micropatterns. The obtained phenotype though does not seem to be entirely compatible with the lack of PS induction observed in C11 micropatterns. If it is true that the exposed animals have a short-tail phenotype, this does not resemble a scenario where BRA (ntl in fish) is not expressed (see PMID: 8402905). In fact, especially for Ulixertinib, axial elongation is still present and somitic rows and other PS derivatives are visible. This might have various sources of explanation, but the authors should clarify why. The authors could also check earlier stages (e.g. late gastrula) to assess the initial formation of early mesoderm (which should lack or be aberrant if the drug induces a failure in PS induction). Overall, the teratogen classification is correct, but the phenotypic prediction is not accurate (or may vary between species).\n9) WNT diffusion modeling\n- WNT diffusion as a function of tissue density is a very interesting hypothesis and might explain the evolutionary need for a correct PS size during gastrulation. This reviewer’s doubts related to cell death during differentiation (see comment above) remains in this assay as well, and this potentially confounding factor should be ruled out before interpreting the phenotype only on the ability of WNT to diffuse (especially because the lower diffusion of WNT was imposed in the mathematical model).\n- Density correlates with WNT waves: this represents the main finding of the manuscript, which also has the potential to inform some basic mechanism during early development. It is indeed an interesting observation, though not completely unexpected. The effect of structure sizes on WNT response has been studied before, also in 3D mouse gastruloids (PMID: 38491138). This work highlighted how varying the size of the initial aggregate does not result in differences in the PS-like domain induction (BRA, CDX, SOX2 patterning). The gastruloid protocol is entirely based on WNT induction using CHIR (small molecule). It could be that the difference in results observed by the two studies relates to the fact that CHIR permeability/diffusion might follow different rules than the WNT itself, being less sensitive to cell density/aggregate size. The authors should compare their results with these and/or try to model WNT vs CHIR diffusion.\n\nMinor comments:\n- The authors focus extensively on the undifferentiated SOX2+ and PS-like BRA+ domains and less on the trophectoderm-like GATA3+/CDX2+ outer ring. Are there molecules that cause specific aberrations in this domain (e.g. its expansion or complete loss)?\n- The authors do not check for endoderm induction, a domain observed in the original protocol (PMID: 24973948). While this reviewer understands the challenges in staining for an extra marker, the authors should comment on the expected outcomes or acknowledge this as a possible limitation of the study (apparently canonical phenotypes might mask defects in endoderm induction which were not measured in this study).\n- The authors should provide high-resolution imaging of some (or all in supplementary figures or repositories) conditions since the micropattern digital imaging panels are very hard to inspect.\n- The link to access the entire morphospace dataset seems to be broken and I could not inspect that closely (lines 199-200).\n- Fig S4 e-h are barely visible and with low image quality, making them quite hard to read.\n- The overall image quality is low and does not allow detailed inspection of the expression patterns and phenotypes. The authors should upload high-resolution figures.\n- Despite being published in ref 45, the authors should detail how the membrane bound β-catenin signal was masked in the live cell imaging experiments, since no dyes to label cell membranes were used (at least not specified in the methods nor in the legends). Clarification is needed given that this is a non-trivial segmentation.\n- Lines 512-513: the use of the term “malformation” to describe an aberrant phenotypic pattern is not ideal. I would recommend using “aberrant gene expression patterns” or similar expressions.\n- The name of the tool FATE-MAP is a bit misleading, since the current study mainly provides information about the ability to exit from pluripotency (one transition towards nascent mesoderm) and does not go beyond that measurement into fate decision-making.", "reviewer_1_attached": "In this work, the authors undertake a formidable effort to screen and quantify the effects of 210 compounds on the development and differentiation of 2D gastruloids from hPSCs, a system well established, primarily by the Warmflash lab. They quantify three fluorescent markers in 10 colonies per treatment, ultimately quantifying images from over 2,000 gastruloids at the single-cell level, a remarkable scale for this kind of analysis. The fluorescence quantification is converted to a vector describing the radial intensity profile for each marker.\n\nComment 1: The idea of generating a mean behaviour image for each cluster is nice, seems to work well in these radially symmetric colonies, and provides an intuitive understanding of each response type. The authors cluster the gastruloid vectors into 12 clusters, assigning 9 of them to normal behaviour and 3 clusters to distinct failure modes (loss of Bra/Sox2 or loss of radial symmetry). Note that converting a profile with a geometric meaning (expression along the radial axis) into a vector (simply by binning) loses all geometric information. For example, under standard metrics the vector 1000 will be equidistant from 01000 and from 00001, while in its geometric interpretation it is clearly closer to the first one. All the same, the clusters obtained seem to be meaningful. The authors may want to discuss this point, and perhaps point to limitations of this embedding. Of course, when we move to network learning in the next step, the network model can re-learn the geometric connections between output nodes.\n\nResponse 1: We agree that a Euclidean distance on binned vectors does not, by itself, encode geometric information (e.g., a peak at position 1 is Euclidean-equidistant from peaks at positions 2 and N). Importantly, our embedding retains the ordered radial sequence, so adjacent positions remain adjacent in the representation. Figure R1 compares (a) Euclidean distance between mean drug-treatment profiles, which separates drug conditions by their average radial patterns, and (b) mean pairwise cosine similarity across all individual drug treated colonies, which emphasizes pattern shape and provides finer separation when treatments produce more heterogeneous profiles. Consistent with the t-SNE embedding, which preserves local neighborhoods rather than absolute distances, the cosine similarity shows tight within-treatment neighborhoods and low similarity across treatments, while the Euclidean distances highlight separation between treatment means that induce distinct global patterns. We revised the manuscript to explicitly acknowledge this limitation and to explain that, while distances are not geometry-aware, the ordered radial sequence is retained.\n\nChanges to the manuscript:\n“Although this embedding does not explicitly preserve geometric distances between adjacent bins, it preserves the exact radial ordering, enabling spatially coherent comparison of radial expression patterns across colonies (Supplementary Fig. 1b-c).”\nFig. R1 has been added as Supplementary Fig. 1b-c.\nFigure R1: Similarity of ordered radial profiles across treatments. (a) Euclidean distance between mean control drug-treatment profiles. (b) Mean pairwise cosine similarity between all control colonies (higher = more similar; color bar inverted). Treatments: mTeSR1, BMP4, CHIR-98014.\n\nComment 2: In the next step, the authors encode the 210 chemical compounds as vectors, and train a network model to predict the gastruloid vector (output) from this encoded input. They also train an additional (smaller) network to embed the predicted vector in tSNE “morphospace embedding”. I assume the prediction errors were done on the original 150-dim vector space, correct? Why is it important to embed into the 2-dim tSNE space? It is only used for visualization purposes. As the name says, it is only a 2-dim embedding of a 150-dim vector, which itself is a representation of a 2-d circular pattern. If the purpose is mostly to evaluate the error (as depicted in Fig. 2b) for the 150-dim space, the embedding process only adds some more distortion to it. It is more informative to show some quantification of the prediction errors (as depicted in the red/green/blue plots), how they behave in training vs. test set, etc.\n\nResponse 2: We agree with the reviewer that quantitative error should be evaluated in the original 150-dimensional phenotype space. Accordingly, we added Supplementary Table 3 reporting RMSE and cosine similarity between predicted and experimental phenotypes for both training and test sets, which provides a rigorous assessment of prediction accuracy. However, learning the 150-dimensional phenotype vector before projecting into morphospace is essential for two conceptual and practical reasons:\n1. Biological interpretability and modularity (clarifying the reviewer’s “distortion” point). The 2D morphospace embedding is not used to evaluate error; all quantitative comparisons are performed in 150D. Its purpose is to co-locate experimental colonies, AI-predicted phenotypes, and PDE-based mechanistic simulations on a shared manifold where developmental trajectories can be visualized and interpreted (as in Fig. 6c). Because any projection can introduce distortion, we explicitly quantified the incremental distortion from projecting to 2D: the learned (parametric) embedding closely reproduces the original t-SNE density field (JSD ≈ 0.09; Supplementary Fig. 2c). Thus, while accuracy is judged in the high-dimensional space, the 2D projection enables direct visualization of how compounds shift along biologically interpretable axes (e.g., cell density, SOX2 stability) that emerge from our mechanistic simulations.\n2. Preservation of nonlinear developmental structure. Predicting the full 150-dimensional phenotype before embedding ensures the model captures biologically meaningful relationships between individual fate markers rather than overfitting to the arbitrary geometry of a particular 2D t-SNE map. The intermediate phenotype vector anchors predictions to measurable biological quantities: radial expression of GATA3, BRA, and SOX2, providing an interpretable format shared with our PDE simulations. Embedding those phenotypes afterward preserves the nonlinear topology of developmental trajectories and enables direct comparison of experimental, simulated, and predicted outcomes within a single interpretable coordinate system. This two-step approach ensures that the morphospace remains stable to retraining and extensible to new experimental data or mechanistic simulations.\n\nChanges to the manuscript:\nAdded Supplementary Table 3 which reports the RMSE and cosine similarity between all chemical structure–to–phenotype predictions and experimental results for both the training and test sets.\n\nComment 3: In the sample images in Fig. 3a we see some bald regions within the colony (which we do not see in the examples in Fig. 2. What are these? Why do we only see them here?\n\nResponse 3: We thank the reviewer for this observation. The “bald” regions visible in some colonies in Fig. 3a are not actual gaps in the cell layer, but rather out-of-focus vertical protrusions. These occur when certain drug treatments cause 3 dimensional cell growth out of the plane of focus of the confocal stack imaged for each colony. This phenomenon is particularly common in retinoids, which have been shown to promote stem cell proliferation.1,2 In the primary high-throughput screen, we acquired shallow z-stacks to keep scan time and data volume tractable across >2,000 colonies; full 3D stacks would have dramatically increased acquisition time and produced many blank z-slices for otherwise flat monolayers.\n\nFig. R2 shows three representative retinoid-treated colonies alongside their radial intensity profiles. The profiles clearly demonstrate that these protrusions are captured as a drop in SOX2 signal in the center of the colony. Nevertheless, these colonies still fall within Cluster 11, which is characterized by a lack of BRA expression. It should also be noted that the radial binning procedure partially mitigates this effect in some colonies, whereas in others the impact on the radial intensity profile is more pronounced. In Fig. R2d, we performed new imaging of tretinoin-treated colonies with the axial range increased to 20 μm; this resolves the apparent “holes” by capturing the out-of-plane cells while preserving the same BRA-loss phenotype\nFigure R2: Representative retinoid-treated colonies exhibiting vertical protrusions due to increased cell density. Example colonies treated with (a) isotretinoin, (b) TTNPB, and (c) bexarotene show lighter central regions that are not gaps in the cell layer but vertical protrusions arising from confinement within the micropattern. Corresponding normalized radial intensity profiles demonstrate loss of mesodermal identity (BRA, red), consistent with Cluster 11 phenotypes. (d) Tretinoin-treated colonies imaged with an increased focal range (20 μm) show filled centers, validating that the “bald” appearance reflects out-of-plane growth rather than absence of cells.\n\nChanges to the manuscript:\n“These colonies frequently displayed lighter central regions that reflect vertical protrusions caused by increased cell density under micropattern confinement (Supplementary Fig. 6)”\nFig. R2 has been added as Supplementary Fig. 6.\n\nComment 4: The authors then run 947 FDA approved drugs through their prediction pipeline, finding 100 to fall within the no-Brachyury cluster. The association statistics corresponds with pre-known risk levels (as teratogens) of those drugs. A few comments on this approach/result:\n1. Looking at Fig. 3c, it is not clear how association of the new compounds to cluster 11 was determined. Some of them (including several category X compounds) seem just as close (or closer) to clusters 8 and 9. Or is it because the association was computed in the 150-d space, and the 2-d projection is just misleading? The same is true for Ulixtertinib prediction.\n2. In this specific example, the association to cluster 11 describes failure in the form of no Brachyury expression. Did we have to go through the sophisticated experimental model (2D gastruloid) and heavy computational framework? Wouldn’t we get a similar (or better) prediction simply screening Bra-expression levels (i.e. a 1d readout instead of 150d readout) in e.g. simple diff protocol applied to a 2D culture? It is worth at least discussing this point.\n\nResponse 4:\nPoint 1: Previously, in the original submission, cluster 11 membership for predicted compounds was determined using a manually defined linear decision boundary in morphospace (Fig. R3a). As the reviewer noted, 2-dimensional t-SNE embeddings can visually make some compounds appear close to multiple clusters. To address this, we have replaced the manual boundary with an objective, data-driven assignment method that uses the same clustering algorithm originally used to identify clusters from the experimental dataset.\nFigure R3: Classification of cluster 11 compounds using a fixed linear decision boundary in morphospace. (a) Original linear classification boundary for cluster 11. (b) Revised morphospace boundaries computed using KDE + watershed clustering, applied consistently to experimental and predicted phenotypes. (c) Benchmarking against 33 FDA-approved compounds previously tested in TeraTox and mEST assays. (d) Ulixertinib and nirogacestat remain classified in cluster 11 under the revised method.\n\nIn brief, the original clustering pipeline first computes a kernel density estimate over the t-SNE space, detects local density maxima, and then applies watershed segmentation to partition morphospace into discrete cluster regions (see Methods). In the revised analysis, we apply these same watershed-defined morphospace boundaries to all predicted drug phenotypes, ensuring that experimental and predicted compounds are classified consistently and without any ad hoc criteria (Fig. R3b). This revision removes subjectivity and directly ties cluster assignment to the same unsupervised morphological structure learned from the experimental dataset. Corresponding language has been added to the Results and Methods sections.\n\nThe analysis comparing cluster 11 to FDA pregnancy-risk label categories has also been revised. While FDA labeling provides clinical guidance, it is influenced by heterogeneous data sources, risk communication practices, and historical regulatory decisions, making it an imperfect ground truth for mechanistically defined teratogenic potential. In the revised manuscript, we have removed the FDA-label analysis and instead benchmark predictions against a curated panel of 33 drugs with known teratogenicity potential, supported by either human data (as reported in FDA drug labeling) or in vivo embryo–fetal development (EFD) studies in rats and/or rabbits. Crucially, these 33 drugs have also been tested in 2 different stem-cell-based assays for teratogenicity assessment: TeraTox (3D iPSC multilineage differentiation) and the mouse embryonic stem cell test (mEST). Only 4 of these compounds were present in our experimental dataset; the remaining 29 were classified fully in silico using FATE-MAP.\n\nPredicted 150-dimensional phenotypes from 33 compounds were embedded into morphospace and classified using the fixed watershed-defined cluster boundaries derived from the experimental dataset (Fig. R3c). Embeddings falling within failure-mode clusters (C1, C11, C12) are assigned positive teratogen status, whereas embeddings within the canonical patterning region (C2–C10) are assigned non-teratogenic status. Across this set, FATE-MAP achieved 72.7% accuracy, outperforming the mEST assay (51.5%) and matching the TeraTox assay (72.7%) while exhibiting higher balanced accuracy than TeraTox (74.8% vs. 66.3%). These results illustrate the potential of FATE-MAP as a scalable screening framework, particularly because our approach does not require experimentally testing each compound, enabling scalable in silico teratogenicity prediction.\n\nFinally, the two compounds selected for follow-up zebrafish validation, ulixertinib and nirogacestat, are still robustly classified within cluster 11 under the revised clustering procedure (Fig. R3d, Revised Fig. 3), confirming that the conclusions of the study are unchanged by the improved classification approach.\n\nChanges to the manuscript:\nRevised manuscript to include benchmarking of FATE-MAP predictions against 33 FDA-approved compounds with experimentally established teratogenicity outcomes in stem-cell-based assays.\nFig. R3c-d has been added as Revised Fig. 3c and Fig. 3e.\n“Of note, drugs were assigned to clusters not by visual inspection of their position in morphospace, which can be ambiguous, but using the same watershed-derived morphospace boundaries defined during the original clustering, ensuring consistent, data-driven classification across experimental and predicted phenotypes.”\nSupplementary Note 1 has been revised to clearly explain the cluster classification process.\n\nPoint 2: This is a good point, and if our sole goal was teratogen identification, a 2D culture differentiation protocol would be much simpler to implement. In fact, prior work has demonstrated monolayer-directed differentiation of hESCs for teratogen detection, on 71 drug-like compounds with known in vivo effects.3 Rather than using BRA as a marker, the authors used the absence of the mesendodermal marker SOX17 as the readout. However, while highly effective for classification, this single-marker framework does not provide insight into how developmental signaling pathways are perturbed, nor does it offer a systematic way to resolve diverse failure modes.\n\nIn contrast, FATE-MAP was designed not only as a teratogen detection platform but as a framework to understand how chemical perturbations reshape developmental trajectories. There are several key differences between our approach and the monolayer system:\n Controlled starting conditions via micropatterning\nFATE-MAP uses micropatterned colonies, which provide reproducible starting conditions and well-defined boundary constraints. This level of control was essential to our discovery of the role of cell density in modulating WNT signaling and mesoderm formation. In the monolayer study, boundary conditions were not controlled, and the authors themselves note that cultures contained both “dense” and “sparse” regions by day 3. They further observed that fate specification varied drastically between these regions, with high-density areas often staining positive for mesendodermal markers and sparse regions lacking them. While this heterogeneity did not undermine their assay because of the highly sensitive marker chosen (SOX17), it limits the ability to extract mechanistic insights about the contribution of density or other developmental parameters.\n Robustness across markers\nThe authors of the monolayer study also acknowledge that “in our hands, the expression levels of T-brachyury varied greatly among cell populations within a single well,” highlighting the variability of monoculture-based platforms. By contrast, FATE-MAP leverages multiparametric readouts from three lineage markers (BRA, SOX2, and GATA3) and encodes their spatial distributions, allowing us to resolve failure modes beyond the absence of a single marker.\n Embryo-relevance and morphological context\nUnlike monolayer differentiation, 2D gastruloids are a validated in vitro model of the gastrulating human embryo, displaying the correct radial organization of germ layers and recapitulating single-cell RNA-sequencing signatures that map to in vivo human gastrulation. This morphological and transcriptional fidelity provides biological context for interpreting compound effects, allowing us to ask how perturbations alter embryonic patterning itself. By contrast, monolayer assays, although highly effective for teratogen classification, lack morphological significance and thus cannot inform how compounds impact the spatial organization of developmental programs.\n\nNevertheless, we thank the reviewer for raising this point. The monolayer-directed differentiation approach was not discussed in our initial manuscript and represents an important contribution to the field. We also acknowledge that its simplicity in both experimental setup and single-marker readout is compelling for applications focused purely on teratogen classification. However, our goal with FATE-MAP is broader: to develop a mechanistically interpretable platform capable of resolving diverse failure modes and revealing how drug perturbations reshape developmental signaling pathways. To highlight this distinction, we have added a discussion of the monolayer approach to the revised manuscript.\n\nChanges to the manuscript:\n“In parallel, monolayer-directed differentiation protocols have also been implemented for teratogen detection, achieving high classification accuracy from single-marker readouts.57 While compelling in their simplicity, these assays lack morphological context and mechanistic resolution, underscoring the added value of FATE-MAP in resolving developmental failure modes.”\n\nComment 5: The authors then devise a mechanistic model, that models the spatio-temporal dynamics of Wnt, BMP and Nodal activity, and tie it with the measureable marker expressions (that serve as a proxy for cell fate) through a simple logic-based model. They demonstrate how such a model can predict the spatial patterns observed, as well as their dependence on cell density. They experimentally validate one implication from the model simulations, which is the inverse dependence of Wnt signaling on cell density.\n\nThey then parametrize morphospace by projecting gastruloids predicted by the mechanistic model while varying two parameters: cell density and SOX2 stability.\n1. It is not clear (or well motivated) why SOX2 stability was chosen – this is not a parameter that can be directly controlled experimentally, and it is not clear to what degree it naturally varies between the gastruloids.\n2. The authors claim that Sox2 stability (Cs) “emerged as a key regulator of fate boundary positions in our model.“ This claim is a bit misleading – the model was devised in a way that gives Cs a high weight, determining the resulting patterns. For example, the stabilities of Gata3 or Brachyury do not appear is parameters in the model at all. This should be better explained and motivated.\n\nResponse 5: We thank the reviewer for raising this important point. To address this, we expanded our analysis to systematically vary all parameters (including the addition of parameters such as GATA3 stability and BRA stability) in our cell fate rules across biologically plausible ranges. A schematic of the expanded rules (Fig. R4a) and a table of parameter definitions (Fig. R4b) are provided below. Figure R4c–d shows the results from our initial manuscript when SOX2 stability (Cₛ) was modulated across multiple cell densities, producing a set of trajectories that spanned orthogonal axes of variation in morphospace. As noted by the reviewer, Cₛ is not directly controllable experimentally. However, it serves as a mechanistically interpretable parameter that scales the persistence of SOX2 against competing lineage programs (GATA3 and BRA). This parameter effectively represents the balance between pluripotency maintenance and lineage commitment over time, and thus provides a natural axis for linking drug-induced shifts to signaling perturbations. While not directly measurable, Cₛ likely reflects integrated regulatory processes such as post-translational modifications (e.g., phosphorylation, methylation) and chromatin accessibility that influence SOX2 degradation and stability.\nFigure R4: Expanded parameterization of mechanistic cell fate model to evaluate regulatory contributions. (a) Equations defining GATA3, BRA, and SOX2 as functions of all possible stability terms, repression weights, and activation thresholds. (b) Table of model parameters with biological interpretation. (c-d) Simulated colony embeddings as a function of cell density and SOX2 stability (Cₛ), showing orthogonal axes of variation.\n\nTo evaluate whether this interpretability was unique to Cₛ, we extended the parameter sweep to include all stability terms, repression weights, and activation thresholds (Fig. R5). Interestingly, none of these other parameters produced trajectories that aligned with major axes of morphospace variation, with the partial exception of C₁. \nFigure R5: Systematic variation of all model parameters across morphospace. Simulated colony embeddings in morphospace while varying each parameter across a biologically plausible range (0.25–1.75) and at multiple cell densities. Among all parameters tested, only Cₛ (Fig. R4) and cell density generate trajectories that align with the major axes of developmental variation, highlighting their interpretability for linking drug-induced shifts to developmental mechanisms.\n\nC₁ (baseline SOX2 level) does generate a trajectory, but this reflects a trivial shift in the initial condition of the system rather than a regulatory process. Mathematically, C₁ sets SOX2(0) at t=0 in the absence of any signaling input. While this rescaling can alter final outcomes, it does not capture the dynamic regulation of boundary positions during gastrulation-like patterning. Furthermore, all stem cells were harvested from the same culture well before being seeded into micropatterned plates, making it unlikely that they possessed significantly different baseline levels of C₁. In contrast, Cₛ multiplies the entire SOX2 term, meaning it governs the effective persistence of SOX2 in the presence of antagonistic inputs (GATA3, BRA). This makes Cₛ a more appropriate analog for drug action, since most perturbations alter stability or turnover rates of transcriptional regulators rather than their absolute initial levels.\n\nIt is also worth noting that parameters like Wnt do produce effects, but these are qualitatively different: for example, lowering or raising Wnt shifts BRA induction to become less density-dependent, resulting in all colonies at a given density to have higher or lower levels of BRA expression, respectively. However, these changes do not project along a clear developmental axis in morphospace and thus are less useful as organizing principles.\n\nIn summary, while we acknowledge that our initial emphasis on Cₛ may have overstated its uniqueness, our expanded parameter sweeps support the conclusion that Cₛ provides the clearest axis for interpreting morphospace displacements. This choice is mechanistically plausible from a drug MOA standpoint: compounds can modulate the effective stability of SOX2 via epigenetic or post-translational control (e.g., phosphorylation creating degron motifs that recruit ubiquitin–proteasome turnover, and lysine/arginine methylation altering DNA binding, nuclear retention, and degradation kinetics). In parallel, shifts in SOX2 abundance have been shown to reconfigure WNT/β-catenin target engagement, biasing β-catenin toward pluripotency enhancers when SOX2 is high and toward mesendodermal gene programs when SOX2 is reduced, thereby switching the net output of WNT signaling (ref. 71).\n\nWe appreciate the reviewer’s point and have tempered our language accordingly; we now present as a mechanistically interpretable axis rather than a unique determinant. We have clarified this rationale in the revised text and included the expanded parameter sweep in the revised Supplementary Information.\n\nChanges to the manuscript:\nAdded Figs. R4 and R5 as the revised Supplementary Fig. 11 and Supplementary Fig. 12.\n“To assess whether other coefficients could generate comparable trajectories, we systematically swept all model parameters over biologically plausible ranges; only density and SOX2 stability consistently produced extended, coherent axes in morphospace (Supplementary Figs. 11-12)”\nAbstract: Changed “Using FATE-MAP, we flagged two clinical molecules as potential teratogens and also revealed cell density and SOX2 stability as determinants of gastruloid patterning” to “Applying this framework, we flagged two clinical molecules as potential teratogens and identified two parameters, cell density and SOX2 stability, that form orthogonal morphospace axes along which canonically patterned gastruloids systematically vary.”\nRevised section header “Parameterizing gastruloid morphospace: cell density and SOX2 stability define developmental axes” to “Parameterizing gastruloid morphospace: cell density and SOX2 stability provide interpretable axes”.\nChanged “Although not directly measurable, this parameter (SOX2 stability) emerged as a key regulator of fate boundary positions in our model” to “Although not directly measurable, this parameter provides a mechanistically interpretable readout that broadly reflects drug mechanism and influences a colony’s position within morphospace”.\nChanged “The resulting structure suggests that morphospace can be decomposed into two interpretable axes defined by these parameters” to “The resulting structure reveals two interpretable axes aligned with these parameters, while not precluding additional contributing factors.”\nChanged “This integration facilitates the identification of key signaling parameters, such as cell density-modulated Wnt signaling and SOX2 stability, that organize developmental trajectories” to “This integration facilitates the identification of key signaling parameters, such as cell density-modulated Wnt signaling and SOX2 stability, that help organize developmental trajectories.”\n\nComment 6: The idea of parametrizing the abstract UMAP embedding by varying interpretable parameters is nice, as it adds to the interpretability of the morphospace.\n\nResponse 6: We thank the reviewer for this positive comment.\n\nComment 7: The authors present the AI and mechanistic models as if they were two alternatives to choose from. This is a bit misleading, as the mechanistic model does not incorporate the potential effects of teratogens. If I understand the rationale of the authors, in integrating (and comparing) the two, they assume that the effect of the teratogens can be summed up through their effect on just 2 parameters (cell density and SOX2 stability).\n1. What is the rationale for that?\n2. Are the effects of the compounds on density simply due to toxicity? (this could be cross-checked against simple growth curves in 2D cell culture for these compounds).\n3. And if they affect Sox2 stability (as rationalized in the Discussion), don’t they have similar effects on GATA3 and Brachyury? Or on Wnt, Nodal and BMP proteins?\n\nResponse 7:\nPoint 1: We do not assume that the effects of all teratogens can be fully captured by modulating just two parameters. Rather, our key insight is that canonical patterning outcomes, those in which the spatial organization of germ layers remains intact, can be largely explained by variation along two axes: cell density and SOX2 stability. These parameters reflect continuous, tunable aspects of developmental regulation that define the structure of the canonical morphospace (clusters C2-C10).\n\nIn contrast, teratogenic compounds fundamentally change the topology of morphospace, pushing colonies out of the canonical density–SOX2 stability continuum and into failure-mode clusters (C1, C11, or C12). Such drugs perturb the rules that underlie pattern formation. For example, they can alter morphogen production or transmission (for example, by directly inhibiting BMP, Wnt, or NODAL). They can disrupt cell–cell contacts that regulate boundary formation and receptor localziation, including effects secondary to drug-induced cell loss. They can also rewire fate-specification rules (for example, by changing competence thresholds or the cross-repression between SOX2, GATA3, and BRA through epigenetic or post-translational mechanisms) so that identical inputs are decoded into different fates.\n\nTo further illustrate this point, we selected two morphogen-pathway inhibitors from our screen, WIKI4 and LY3200882, that produced colonies in cluster 11 for follow-up analysis. WIKI4, a tankyrase inhibitor that suppresses WNT/β-catenin,4 was modeled by eliminating WNT auto-activation (setting s = 0). LY3200882, a selective ALK5 (TGF-βRI) inhibitor that attenuates the uSMAD2/3 arm shared by TGF-β/Activin/Nodal,5 was modeled by eliminating the NODAL/Activin production term (setting s = 0). With all cell-fate rules unchanged, both perturbations yielded BRA-null simulated colonies that aligned with the cluster 11 experimental phenotypes at comparable densities (Fig. R6).\nFigure R6: Mechanistic modeling of teratogen-induced failure modes. (a)–(b) WNT inhibition modeled by eliminating WNT auto-activation (s = 0) recapitulates the effect of WIKI4, yielding uBRA-null colonies consistent with experimental cluster 11 phenotypes. (c)–(d) NODAL inhibition modeled by eliminating Activin/NODAL production (s = 0) reproduces the effect of LY3200882, likewise producing BRA-null colonies that align with experimental cluster 11 outcomes.\n\nThese examples help illustrate that the mechanistic modeling framework is not a universal model of drug action. Rather, it provides a compact description of canonical variation and, when supplied with specific mechanistic hypotheses (e.g., WNT or ALK5/SMAD2/3 inhibition), can also reproduce failure modes, though this requires prior knowledge of pathway targeting.\n\nOverall, the mechanistic model codifies the rules that govern the structured, canonical region of morphospace; it tells us that when a compound drives colonies outside this region, it is operating under a different set of rules. The AI model, in turn, learns patterning directly from chemical structure, predicting where a compound will land in morphospace. Together, the two are complementary: the AI model flags and localizes compound effects, and the mechanistic model interprets whether those effects reflect canonical displacements (density/C s) or topology-altering failures, yielding insight into how a teratogen may act (i.e., inhibition of Wnt signaling). In the revised manuscript, we have rewritten the relevant section in the Results section to articulate this distinction more clearly.\n\nChanges to the manuscript:\n“While our AI-based structure-to-phenotype model predicts 2D gastruloid phenotypes with high accuracy (Fig. 6a), it functions as a black box and does not directly reveal which developmental processes are perturbed. To provide mechanistic insight, we developed a complementary model grounded in morphogen diffusion and fate transitions (Fig. 6b). By projecting simulated phenotypes into the same learned manifold as experimental data, this model assigns interpretable structure to morphospace and enables causal reasoning about how perturbations modulate underlying parameters. In practice, the mechanistic model codifies the rules that govern the structured, canonical region of morphospace and, when supplied with prior mechanistic hypotheses (e.g., morphogen inhibition; Supplementary Fig. 10c), can also reproduce specific failure modes. The AI model, in turn, predicts where a compound will land in morphospace directly from chemical structure. Together, the two are complementary: the AI localizes compound effects, and the mechanistic model interprets whether they reflect canonical displacements or topology-altering failures.”\nWe have added Fig. R6d to Supplementary Fig. 10c.\n\nPoint 2: We appreciate the reviewer’s thoughtful comment. A central finding of this work is that cell density defines a continuous axis of morphospace variation, influencing fate patterning through biophysical and signaling-mediated mechanisms rather than simply reflecting toxicity. This relationship was first observed from the embedding of experimental colonies (Supplementary Fig. 9), where the slope and fit of the BMP4-control trendlines, together with mean cluster representations, revealed systematic, density-dependent shifts within the canonical patterning region. We then developed and validated a density-dependent model of Wnt signaling (Fig. 4) and subsequently identified cell density as a principal morphospace axis (Fig. 5).\n\nWhile conventional 2D growth-curve assays could, in principle, provide complementary information about toxicity, our screen captures this in a physiologically relevant format that reflects edge-localized receptor dynamics and BMP→Wnt→Nodal feedbacks unique to the 2D gastruloid model. By counting every nucleus within each colony, we directly quantify cell number at single-colony resolution and can identify gross reductions in viable cell density without assuming a mechanism. Because all wells start at the same seeding density (50,000 cells per well), large drops in final nuclear counts primarily reflect drug-induced toxicity and exceed what would be expected from growth arrest alone.\nFigure R7: Density-based toxicity threshold. (a) Mean nuclei per colony for each BMP4+drug condition. Toxicity threshold (red, ≤50% mTeSR mean) and hyperproliferation threshold (green, >mTeSR mean). (b) Toxic compounds map predominantly to failure modes C1, C11, and C12; insets show representative colonies. The embedding illustrates that overt cytotoxic collapse can be distinguished from density-dependent modulation within the canonical patterning region.\n\nTo objectively define toxicity, we used a control-anchored cutoff: treatments with mean nuclear counts <50% of the mTeSR control mean were classified as toxic. As shown in Figure R7a, seven compounds met this criterion, while seven others exceeded the mTeSR1 mean and formed a high-density (hyperproliferative) group. Notably, 96.1% of colonies exposed to toxic compounds mapped to failure-mode clusters C1, C11, and C12 (Fig. R7b), demonstrating that our assay distinguishes toxicity-driven collapse from density-dependent modulation of canonical patterning.\n\nThese analyses show that while drug-specific perturbations can alter proliferation or morphogenesis, the underlying patterning program remains remarkably robust, preserving canonical organization across a broad range of cell densities.\n\nChanges to the manuscript:\nAdded Fig. R7 as Supplementary Fig. 3.\n“Because failure modes can arise from either differentiation failure or compound-induced cytotoxicity, we quantified total nuclear counts for every colony and identified a small subset of drugs that reduce viable cell number below 50% of the mTeSR-only control (Supplementary Fig. 3). This analysis separates genuine patterning failures from cytotoxicity-driven collapse within the failure-mode regions.”\n\nPoint 3: We agree that drug perturbations are unlikely to act exclusively on SOX2. Indeed, in our expanded parameter sweeps (Fig. R4-5; Supplementary Fig. 10-11), we systematically varied all model coefficients, including GATA3 stability, BRA stability, and morphogen thresholds. Some of these perturbations did shift colonies within morphospace; however, none produced smooth, monotonic trajectories that traversed the canonical region of morphospace.\n\nThis suggests that while drugs may act at many levels (e.g., modulating WNT, NODAL, or BMP production, or altering GATA3/BRA turnover), their downstream effects converge on parameters that influence the persistence of SOX2 relative to competing fates. In our framework, SOX2 stability (C s) captures this convergence: it reflects the effective stability of SOX2 expression in the presence of GATA3 and BRA, integrating diverse upstream influences into a single, interpretable readout.\n\nWhile further work will be needed to dissect the precise molecular mechanisms by which specific compounds alter SOX2 stability, the fact that a single parameter can organize heterogeneous perturbations within morphospace is powerful for interpreting developmental outcomes. In the revised manuscript, we temper claims of uniqueness and explicitly acknowledge that SOX2 stability serves as a useful interpretable axis rather than an exclusive determinant.\n\nChanges to the manuscript:\nAbstract: Changed “Using FATE-MAP, we flagged two clinical molecules as potential teratogens and also revealed cell density and SOX2 stability as determinants of gastruloid patterning” to “Using FATE-MAP, we flagged two clinical molecules as potential teratogens and identified two parameters, cell density and SOX2 stability, that form orthogonal morphospace axes along which canonically patterned gastruloids systematically vary.”\nChanged “Although not directly measurable, this parameter (SOX2 stability) emerged as a key regulator of fate boundary positions in our model” to “Although not directly measurable, this parameter provides a mechanistically interpretable readout that broadly reflects drug mechanism and influences a colony’s position within morphospace”.\nChanged “The resulting structure suggests that morphospace can be decomposed into two interpretable axes defined by these parameters” to “The resulting structure reveals two interpretable axes aligned with these parameters, while not precluding additional contributing factors.”\nChanged “This integration facilitates the identification of key signaling parameters, such as cell density-modulated Wnt signaling and SOX2 stability, that organize developmental trajectories” to “This integration facilitates the identification of key signaling parameters, such as cell density-modulated Wnt signaling and SOX2 stability, that help organize developmental trajectories.”\n\nComment 8: In summary, the effort to generate a prediction system for the effects of compounds on human embryonic development is commendable. The AI-generated predictions seem to generate useful predictions. The integration with the mechanistic model, while adding mathematical gravity to the manuscript, is not very clear currently.\n\nResponse 8: In the revised manuscript, we have clarified the rationale and scope of the mechanistic model. Specifically, we emphasize that the two-parameter framework (cell density and SOX2 stability) provides a compact, interpretable description of canonical variation in gastruloid patterning, while not serving as a universal model of drug action. Through expanded parameter sweeps (Fig. R4-5), we show that other coefficients can shift colony outcomes but do not generate coherent developmental trajectories across morphospace. We also demonstrate, with mechanistic simulations of WNT and NODAL inhibition (Fig. R6), that the model can reproduce teratogen-induced failure modes when supplied with prior pathway knowledge, complementing the AI model’s ability to predict phenotypes directly from chemical structure. Together, these revisions make clear that the AI and mechanistic approaches are not alternatives but complementary: the AI localizes compound effects, and the mechanistic model interprets whether they reflect canonical displacements or topology-altering failures, thereby providing a coherent framework for understanding how drugs reshape developmental trajectories."}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have satisfactorily addressed all my concerns. I believe the MS is now ready for publication.", "reviewer_2": "The authors have sufficiently addressed nearly all my concerns.\nMy only remaining major concern is the modeling of the density dependence of Wnt. In the first round of feedback, I mentioned their model appears to predict a delay in signaling that is not visible in the data. In their reply, the authors argue this is a matter of normalization of the data and changed the normalization. However, looking at the new normalization (Fig.R14), there is still a qualitative difference between simulation and experiment, now in the amplitude: while simulations all approach the same maximal signaling level over different time scales, the amplitudes are clearly density dependent in the data. Overall, it seems to me the disagreement can’t be normalized away and should be discussed in the manuscript.\nI also had a few more very minor comments:\n• In Figure R8, the different blues for the two different BMP4 groups become indistinguishable when I print it.\n• Why did the authors choose not to show error bars on radial profiles in fig 1?\n• The authors state in their rebuttal that GATA3 follows rapid switch-like dynamics citing Gunne-Braden et al, however this paper disagrees with a large body of literature. Teague et al Nat Comm 2024 observed different GATA3 dynamics and discuss the disagreement.\n• In the revised discussion on ROCK inhibition the authors cite Maldonado et al. In the same place they should also cite Etoc et al 2016 since they used ROCK inhibitor to disrupt tight junctions and increase differentiation as discussed.", "reviewer_3": "The authors have done a very good job to address all the concerned raised. This reviewer thanks them for their great effort.\nThis reviewer still thinks that the translational aspect of the tool prediction should be less emphasized and acknowledged as a limitation of the study and further point of improvement in the Discussion section. This comment comes mainly from the results of the zebrafish experiments (where the predictive power of the tool might not be directly translatable in other systems or real embryos) and from the inaccuracy in the prediction of a subset of compounds (including very well-known teratogens).\nApart from this comment, this reviewer supports the publication of the study without further revision."}, "author_response": "Reviewer 2 (Remarks to the Author)\nThis paper develops a framework for high throughput drug screening and phenotype prediction for 2D gastruloids. Improving on previous work describing teratogenicity screening in vitro, the authors develop an integrated approach of deep learning and mechanistic modeling to predict and interpret phenotypes. The paper contains several nice contributions to the field, including 1) the machine learning model mapping phenotypes into a (somewhat) interpretable 2D space, 2) a large dataset made public of the effect of different drug treatments on 2D gastruloids, 3) an extension of the model from Chhabra et al 2019 to include cell density and explicit modeling of cell fate. However, a major weakness is a lack of rigor and detail, including poorly motivated choices in approach. For a paper about what is essentially a statistical testing and prediction framework, rigor seems particularly important.\n\nMajor comments:\nComment 1: The clarity and accuracy of the abstract could be improved. First, it claims to elucidate “failure modes” but neither the abstract nor the rest of the paper clearly defines what a failure “mode” is. If what is meant by this is that the phenotypes can be captured by cell density and SOX2 stability, perhaps that should be made more explicitly. It also says the paper “revealed” cell density and SOX2 stability as determinants of patterning, but both were already known, the author do in fact cite previous work showing the importance of density (e.g. Etoc et al) and SOX2 stability (Camacho-Aguilar et al, Teague et al). Perhaps corroborate is a better word and the discussion would ideally more explicitly incorporate this past work.\nResponse 1: We thank the reviewer for this constructive feedback. In the revised abstract, we now provide an explicit definition of failure modes, stating “we mapped a phenotypic morphospace that separates canonical patterning, in which primitive-streak fates are correctly specified and radially organized, from failure modes, defined as departures from this organization and marked by a loss of a required fate and/or radial symmetry.” This clarification distinguishes canonical patterning outcomes (colonies in which germ-layer topology is preserved and variation can be captured by variations in cell density and SOX2 stability) from failure modes, which fall outside this density–SOX2 stability continuum and reflect topology-breaking disruptions (e.g., impaired morphogen production or signaling, disrupted cell–cell contacts, or rewired fate-decision rules).\nWe also agree with the reviewer that prior work has highlighted the importance of both density and SOX2 stability. Our contribution is not in re-identifying these regulators, but in showing that they emerge as orthogonal, interpretable axes that structure canonical morphospace and can be used prospectively to parameterize developmental trajectories.\nWe have updated the Discussion to more explicitly integrate these foundational studies and to clarify how our work extends them by embedding density and SOX2 stability as orthogonal morphospace axes, thereby providing a framework for distinguishing canonical variation from failure modes.\nChanges to the manuscript:\nRevised abstract: “…we mapped a phenotypic morphospace that separates canonical patterning, in which primitive-streak fates are correctly specified and radially organized, from failure modes, defined as departures from this organization and marked by a loss of a required fate and/or radial symmetry. … Applying this framework, we flagged two clinical molecules as potential teratogens and identified two parameters, cell density and SOX2 stability, that form orthogonal morphospace axes along which canonically patterned gastruloids systematically vary.\nAdded paragraph to the Discussion: “Previous studies have established the importance of density and SOX2 dynamics in fate control in human 2D gastruloids. Etoc et al. showed that colony density both modulates morphogen accessibility (receptor relocalization and NOGGIN restrict TGF-β signaling to the colony edge) and regulates BRA patterning (as density rises, a BRA band emerges and shifts toward the edge).40 Camacho-Aguilar et al. showed that SOX2 must decline for WNT to drive mesoderm and linked SOX2 loss to BMP/SMAD4 dynamics,46 while Teague et al. demonstrated that SOX2 acts as an integrator of BMP/SMAD4 input, with its rate of decrease proportional to the signaling integral, highlighting its role as a memory of morphogen exposure.43 Building on this foundation, we show that WNT signaling is inversely proportional to cell density and that cell density and SOX2 stability emerge as orthogonal, interpretable axes spanning the canonical morphospace. This framework not only corroborates earlier findings but also supports prospective mapping of developmental trajectories and separation of canonical variation from topology-breaking failure modes.”\n\nComment 2: Rigor should be improved:\n The paper states “10 colonies per drug condition were imaged”. However, it is not clear whether these all came from the same experiment and if so if they came from the same well. Inter-colony variation within a well is very small, between wells on the same plate small, between experiments or different plates within an experiment can be significant. To understand if phenotypes of two drugs are significantly different, or whether the prediction is correct within the experimental uncertainty, it is essential to quantify the inter-experiment phenotypic variation under the same drug perturbations and to make clear whether all drug treatment data was collected in the same experiment and if not to show if there are differences in the control condition between experiments and determine if there are any batch effects.\nResponse 2: We agree that assessing day/plate (“batch”) effects is essential. The 211 drug conditions were run on four 96-well plates on different days; for each condition, fluorescence profiles (GATA3, BRA, SOX2) were robustly normalized per plate by subtracting the plate mean and scaling by that plate’s interquartile range (IQR; 25th–75th percentile), which reduces day-to-day intensity drift while preserving relative differences across drugs. After this normalization, colonies do not segregate by plate in morphospace (Fig. R8a–d), controls evaluated across plates and across wells occupy overlapping regions (Fig. R8e). Finally, a PERMANOVA on the 150-dimensional phenotype vectors showed a large drug effect (R² = 0.836, p = 0.001) and a small plate effect (R² = 0.025, p = 0.001) (Fig. R8f), indicating that residual plate/day variation is minimal relative to drug effects and does not drive the observed clustering or predictions.\nFigure R8. Plate/day effects are minimal after per-plate robust normalization. (a–d) Morphospace of all colonies (gray) with colonies from each plate overlaid and colored by cell density. (e) Control treatments (BMP4-only and mTeSR-only) from the indicated plates/wells occupy overlapping regions of morphospace. (f) PERMANOVA on the 150-dimensional phenotype vectors shows a large Drug effect (R² = 0.836, p = 0.001) and a small Plate effect (R² = 0.025, p = 0.001), indicating residual plate/day variation is minimal relative to drug effects.\nChanges to the manuscript:\nFig. R8f has been added to the revised Supplementary Fig. 1f.\n“Across the broader compound screen, plate-to-plate variability was minimal compared to drug-induced effects (Supplementary Fig. 1f), confirming that the 2D gastruloid model analyzed in this way meets the criteria for constructing a morphospace map.”\n\nComment 3: The untreated (BMP only) control conditions is not shown in most figures but should be for comparison. For example, in which cluster does if fall in Fig 1g and where is it in 1f, and in Fig 6c – how does Raloxifene HCl compare to untreated? Can it be said to be significantly different? The same for all other figures. If I understand it correctly, Supplemental Fig. 5a shows the BMP only control clusters and their natural variation. This covers a large part of the map and the majority of the “Canonical Cell Patterning” clusters in Fig.1g, suggesting that those clusters should perhaps be a single cluster and that most drug treatments within that cluster cannot be distinguished from control.\nResponse 3: We thank the reviewer for emphasizing the importance of visualizing the untreated (BMP-only) controls alongside drug perturbations. Our figures are already dense, and plotting all control colonies in every panel would obscure treatment-specific patterns. To balance clarity with completeness, we have:\n1. Made control placement explicit in the text for the global maps (Fig. 1f–g): BMP4-only controls fall within canonical clusters C2, C3, C4, C5, C8, and C9, consistent with density-driven variation in BRA width/position.\n2. In Fig. 6c, we now overlay the BMP4-only control centroid so the reader can directly compare model predictions with the control.\n3. Pointed to consolidated control maps that show their global distribution and density dependence (Supplementary Fig. 2a for overall location; Supplementary Fig. 9a for the density/BRA continuum).\nRegarding the reviewer’s suggestion that canonical clusters might collapse into a single control cluster: BMP4-only controls indeed span a continuous canonical region because cell density systematically modulates BRA width/position. We keep these as separate clusters in the main analysis because they capture interpretable regimes along the canonical manifold (rather than a single undifferentiated bin), while referring to them collectively as the “canonical” region to avoid over-interpretation. Statistical analysis of cluster significance is addressed in our response to Comment 4.\nChanges to the manuscript:\n“The canonical patterning region (C2–C10) contains colonies that retain all three fate markers in concentric radial order (SOX2 central, BRA intermediate, GATA3 peripheral) and differ primarily in the relative sizes of these domains. BMP4-only controls map exclusively within this region (C2, C3, C4, C5, C8, C9; Supplementary Fig. 2a). In contrast, clusters C1, C11, and C12 lie outside the canonical region and exhibit either the loss of a required cell-fate marker and/or loss of radial symmetry.”\n“Because cell seeding produces inherent variation in cell density across colonies, BMP4-treated controls provide a natural spectrum for analyzing density-dependent effects. Indeed, the distribution of these control colonies in morphospace is explained by their density (Supplementary Fig. 9a),…”\n\nComment 4: More generally, it is not clear the clusters are meaningful. The authors use a custom method based on watershed segmentation. It is not clear why they need this, and the method is not described in sufficient detail to be reproducible. To motivate the approach, it is stated that an advantage is that the number of clusters is not predetermined but there are many methods for which that is true, such as Leiden/Louvain. Furthermore, one can optimize cluster number for methods with fixed cluster number using cluster analysis. Regardless, cluster analysis should be performed here to determine if the obtained clusters are statistically meaningful, e.g. calculating Silhouette coefficients, etc (https://en.wikipedia.org/wiki/Cluster_analysis). This reviewer would not be surprised if the number of statistically and biologically meaningful clusters is a lot smaller than the 12 that are discussed.\nResponse 4: We thank the reviewer for this important point. The choice of watershed segmentation was motivated by a key requirement of our downstream analyses: new predicted phenotypes must be assigned to clusters in a consistent, reproducible, and parameter-free manner. Our density-based watershed segmentation produces spatially contiguous cluster regions in morphospace that act as stable decision regions. This allows post-hoc classification of new embeddings simply by determining which region they fall into, without re-running clustering or re-estimating parameters. In contrast, methods such as Leiden, Louvain, or Spectral clustering generate partition assignments only for the specific dataset used during clustering; when new points are introduced, the clustering must be re-computed, and the solution can vary depending on initialization, neighborhood size, and resolution parameters. For our application (evaluating model predictions and comparing experimental and in silico phenotypes) consistent post-hoc assignment was essential. We have clarified this rationale in the revised Methods.\nFigure R9. Alternative clustering methods support a similar partitioning of morphospace.\n(a) Silhouette score heatmap for Spectral Clustering across combinations of nearest-neighbor graph size (kNN = 10–25) and cluster numbers (k = 2–15). (b) Silhouette score heatmap for Leiden Clustering across kNN values (10–25) and resolution parameters, showing comparable optima. (c) Representative Spectral clustering result at kNN = 15 (k = 9) corresponding to its top Silhouette score (Sil = 0.400), yielding ~9 coherent clusters. (d) Representative Leiden clustering result at kNN = 15 (resolution = 0.2; C = 10) corresponding to its top Silhouette score (Sil = 0.403), yielding ~10 clusters.\nTo evaluate whether the clusters identified by watershed segmentation reflect statistically meaningful structure, we conducted a complementary cluster validity analysis using both Spectral and Leiden clustering across a range of k-nearest-neighbor graph sizes and resolution/cluster-number settings, and computed Silhouette coefficients (Fig. R9). Across both clustering frameworks, the optimal solutions consistently yielded ~9–12 clusters, depending on hyperparameters, with Silhouette scores comparable to the watershed-derived clusters. These results indicate that the embedded morphospace contains multiple separable and biologically coherent phenotypic states and that the number of clusters identified is not artificially inflated by the watershed method.\nA direct quantitative comparison is provided below:\nClustering Method kNN # of Silhouette Davies-Bouldin\nused clusters Score Index\nFATE-MAP (KDE+watershed) N/A 12 0.289 0.895\nSpectral (nearest-neighbors) 15 9 .400 .700\nLeiden (graph clustering) 15 10 .403 .767\nThese results show that while the watershed approach yields slightly lower Silhouette scores than the best-parameterized Spectral/Leiden solutions, this reflects a trade-off between marginal gains in cluster compactness and the substantial practical advantage of having fixed, reusable cluster boundaries that enable consistent classification of new predicted phenotypes.\nChanges to the manuscript:\n“The normalized single-colony profiles were embedded into two dimensions with t-SNE. To define morphological phenotypes, we applied a Gaussian kernel density estimate (KDE) over the t-SNE space of all experimental colonies, identified local density maxima, and then performed watershed segmentation to partition the morphospace into discrete, contiguous regions (Supplementary Fig. 1d). This unsupervised procedure does not require pre-specifying the number of clusters and yields cluster boundaries that reflect the empirical structure of the embedded data. While alternative clustering methods (e.g., Spectral or Leiden) produce similar numbers of clusters, the watershed-defined boundaries offer a key practical advantage in that they can be directly reused to assign new predicted phenotypes without re-running clustering, ensuring consistent and unbiased classification across experimental and in silico data.”\n\nComment 5: Error bars on radial profiles. The analysis throughout the paper only considers the mean radial profile for each colony, but there are error bars on these radial profiles that may also change by perturbation, e.g. the variance in expression for some distance from the edge goes up when the patterns become less symmetric, and generally one could imagine propagating these error bars to the morphospace to get a sense of how significantly different the profiles represented by dots there are.\nResponse 5: We thank the reviewer for this thoughtful suggestion. We have now added error bars representing the standard deviation across all experimental colonies for each drug treatment when plotting the radial intensity profiles (revised Figs. 2, 5, and 6). These error bars convey the extent of colony-to-colony variability and allow direct visual comparison between the predicted phenotype and the range of experimentally observed outcomes. In addition, when evaluating prediction accuracy, we now report whether the model-derived radial profiles fall within the experimental variability margins (Supplementary Table 6), enabling a clearer interpretation of the biological significance of prediction errors.\n\nComment 6: In testing the prediction framework, the separation into training and testing data is essential. How was the 10% of compounds held out for test selected? Randomly? It would be good to see N-fold cross-validation to see how much performance varies as different compounds are held out for validation.\nResponse 6: We thank the reviewer for this helpful suggestion. In the revised manuscript, we have strengthened both the validation dataset and the model evaluation strategy. First, we revised the external validation approach. In the original submission, model predictions were compared to FDA pregnancy-risk labeling. However, FDA labeling aggregates heterogeneous evidence sources (human case reports, post-marketing surveillance, historical risk communication practices, and regulatory conventions) and therefore does not provide a consistent mechanistic ground truth for teratogenic potential.\nTo address this, we now benchmark FATE-MAP against a curated panel of 33 FDA-approved compounds with well-characterized teratogenicity outcomes, supported by either human clinical evidence or standardized embryo–fetal development (EFD) studies in rats and rabbits. Importantly, all 33 compounds have also been evaluated in two stem-cell–based assays designed for teratogenicity testing: TeraTox (3D iPSC multilineage differentiation) and the mEST (mouse embryonic stem cell test) assay.\nOur structure-to-phenotype model first generates a 150-dimensional predicted phenotype vector for each compound and then projects this vector into the morphospace, where cluster identity is assigned using the same watershed-derived decision boundaries defined in the original clustering. Compounds whose predicted embeddings fall into the failure-mode clusters C1, C11, or C12 are labeled teratogenic, while those in canonical patterning clusters C2–C10 are labeled non-teratogenic (Fig. R10c). \nOf the 33 benchmark compounds, only four were present in our experimental dataset; the remaining 29 were predicted fully in silico. Across this set, FATE-MAP achieved 72.7% accuracy, outperforming mEST (51.5%) and matching TeraTox (72.7%), while exhibiting higher balanced accuracy than TeraTox (74.8% vs. 66.3%; Fig. R10d). These results demonstrate that FATE-MAP can serve as a structurally grounded, scalable framework for in silico teratogenicity screening.\nRegarding train/test partitioning:\nThe held-out compounds in the main analysis were selected at random (10% of compounds). To assess whether performance depends on the particular composition of the hold-out set, we performed additional resampling analyses in which we repeatedly re-trained the model while holding out different randomly selected subsets of compounds. Representative results are shown in Fig. R11. Across resamples, performance on the 33-compound benchmark panel remained within a similar range, indicating that model accuracy is not driven by a specific train/test split. Furthermore, within each training split, the model is trained using 5-fold cross-validation, and predictions are generated by ensembling the best-performing fold models. This procedure improves robustness in a limited data regime and reduces variance arising from any single training subset. Together, the 5-fold cross-validation and the resampling analysis demonstrate that the predictive performance of FATE-MAP is stable across different train/test partitions and is not dependent on the particular held-out compounds.\nFigure R11. Predictive performance is stable across different training/testing splits. (a–f) Representative evaluation results on the 33-compound benchmark panel, each corresponding to a separate model trained with a different randomly selected 10% of compounds held out. In each split, predicted embeddings were classified in morphospace using the fixed watershed-defined cluster boundaries, and accuracy reflects agreement between predicted and known teratogenicity outcomes. Across splits, overall accuracy remains within a similar range (≈66–73%), indicating that model performance is robust to variation in train/test partitioning.\n\nComment 7: The authors use RMSE as a measure of goodness of fit, but these numbers have no absolute meaning. Cosine similarity is also hard to interpret as good or bad. More meaningful would be to see if prediction and data agree within the experimental error.\nResponse 7: We agree that RMSE and cosine similarity by themselves do not convey whether prediction accuracy falls within the range of natural experimental variability. To address this, we have now incorporated a coverage metric that quantifies, for each compound, the fraction of the 150 phenotypic features for which the model’s predicted value lies within one standard deviation of the empirical replicate distribution. This measure provides an interpretable, data-driven benchmark: if predictions fall within one standard deviation, they are statistically indistinguishable from replicate-to-replicate biological variation.\nUsing this approach, we find that on average 74% of features per compound in the training set fall within experimental variability, and 54% of features in the blinded set fall within variability. Across all 209 compounds, the global mean coverage is 72%, demonstrating that the model routinely predicts phenotypes at a level comparable to experimental noise. We have added these statistics to the Results and summarize them alongside RMSE, cosine similarity, and Pearson correlation coefficients. These values allow for a more interpretable assessment of predictive accuracy and directly address the concern that model performance should be evaluated relative to biological variation rather than abstract error metrics. As noted in Response 5, we also now include error bars and variability ranges directly on the radial profile plots, allowing readers to visually assess whether predicted profiles fall within the empirically observed variability.\nChanges to the manuscript:\n“The model accurately predicts 150-dimensional phenotype vectors and corresponding morphospace positions for both the training and blinded compounds, capturing canonical and failure-mode clusters. Across 188 training compounds, the mean RMSE was 0.17, with a mean cosine similarity of 0.90 (mean Pearson r = 0.79). Importantly, on average 74% of features per compound fell within one standard deviation of the empirical replicate distribution, indicating that most predicted feature values lie inside the observed biological variability. On the blinded set (n = 21), performance remained strong (mean RMSE = 0.24; mean cosine similarity = 0.85; mean Pearson r = 0.70), with 54% of features per compound within one standard deviation of replicate variation. Aggregating all 209 compounds, the overall means were RMSE = 0.17, cosine similarity = 0.90, Pearson r = 0.78, and 72% within one standard deviation. Representative predictions for all 12 clusters are shown in Fig. 2b, including both phenotype profiles and projected morphospace coordinates.”\n\nComment 8: In figure 3, the authors compare their prediction to a “ground truth” drug classification. They should show at least a confusion matrix with accuracy, % false positive, false negative, perhaps also check other measures like precision. Systematically showing %enrichment of drug classes in clusters would also be good. Looking at Figure 3b it looks well mixed.\nResponse 8: In the revised manuscript, we now provide a full confusion matrix and report standard binary classification metrics, including accuracy, precision, sensitivity, specificity, and balanced accuracy, for the curated 33-compound validation panel. These results are shown in the revised Fig. 3d (Response Fig. R10). This provides a transparent summary of false positives and false negatives and allows readers to quantitatively interpret model performance beyond the overall accuracy value.\nAdditionally, to evaluate the model’s ability to generalize beyond the training distribution and to benchmark against other computational toxicity predictors, we carried out an independent comparative analysis across 20 compounds held out from training. We compared FATE-MAP to two established in silico developmental toxicity models: the CAESAR developmental toxicity predictor and the Procter & Gamble (P&G) teratogenicity model. For CAESAR and P&G, we used the originally published model outputs, ensuring a fair comparison without retraining or parameter adjustments. Across this set, FATE-MAP achieved an accuracy of 75.0%, matching CAESAR (75.0%) and substantially outperforming the P&G model (35.0%), while also demonstrating the highest balanced accuracy (85.7% for FATE-MAP vs. 81.2% for CAESAR) (Supplementary Fig. 7a–b; Fig. R12). These results indicate that FATE-MAP not only performs competitively with existing experimental assays but also surpasses or matches leading computational developmental toxicity models.\nFigure R12 (Revised Supplementary Fig. 4a-b in manuscript). Comparative performance of FATE-MAP and existing in silico teratogenicity models. (a) Predicted morphospace embeddings for the 20 compounds held out from training, with points labeled according to whether the model’s teratogenicity classification matched the known outcome. Boundaries correspond to the failure-mode (pink) and canonical (blue) morphospace regions. (b) Confusion matrices and performance metrics comparing FATE-MAP to two widely used in silico developmental toxicity predictors (CAESAR and P&G) on the same 20-compound benchmark. FATE-MAP matches CAESAR in overall accuracy (75%) and achieves the highest balanced accuracy (85.7%), substantially outperforming the P&G model across all metrics.\n\nComment 9: For figure 3, the authors focus on ulixertinib and nirogacestat for verification. Why are these not shown among the data points in Figure 2? Looking at the Supp. Fig. 4ab, the colony images show dark blue lines streaking along and across the colonies. These do not look like nuclear stains for the expected markers. What are they? If they are not cells, how do they affect analysis? are they excluded? If they are cells, why does their pattern look so strange?\nResponse 9: Ulixertinib and nirogacestat were selected for follow-up validation because both compounds were predicted by FATE-MAP to induce distinct and reproducible gastrulation failure phenotypes representative of cluster 11 (loss of BRA). These compounds were not explicitly labeled in Fig. 2 to maintain visual clarity given the large number of compounds displayed; however, both are included among the data points.\nRegarding the reviewer’s observation of dark blue streaking lines in Supplementary Fig. 4a–b, the newly included Fig. R13a–c provides a clearer view of these features. As shown in Fig. R13a, the blue streaks appear as non-nuclear artifacts likely arising from non-specific binding of the secondary antibody rather than true biological signal, since they are not co-localized with DAPI-stained nuclei. Our quantitative analyses are based exclusively on nuclear intensities of the three lineage markers (GATA3, BRA, and SOX2), and thus these signals are effectively excluded from downstream analysis.\nSpecifically, our image analysis pipeline first performs DAPI-based nuclear segmentation and computes the radial position of each nucleus relative to the colony edge (Fig. R13b). Using the resulting masked colony image (Fig. R13c), we then calculate the mean fluorescence intensity for GATA3, BRA, and SOX2 across 50 concentric annular bins (~5 µm width) spanning from the periphery to the center. This yields a 150-dimensional vector capturing the azimuthally averaged nuclear expression pattern of the three markers as a function of radial position, ensuring that any non-nuclear background signal (including the blue streaks) does not affect quantification.\nFigure R13: Image analysis pipeline for nirogacestat-treated colonies. (a) Full colony image showing non-nuclear blue streaking, likely due to non-specific secondary antibody binding. (b) DAPI-based nuclear segmentation and radial mapping of each nucleus relative to the colony edge. (c) Masked colony image used to quantify mean nuclear intensities for GATA3, BRA, and SOX2 across 50 concentric annular bins (~5 µm width), generating a 150-dimensional vector describing spatial marker expression profiles.\n\nComment 10: As stated above, the modeling in Fig. 4 is one of my favorite parts of this paper, but it does raise some questions that if addressed would further strengthen the manuscript. The authors model Wnt production as density-dependent and verify the predicted Wnt signaling dynamics. This nicely captures the narrower range of Wnt signaling at higher density, but their model also shows a later onset of Wnt signaling in Fig.4c which the experimental data in Fig.4e does not seem to show. Can this be explained? Can the model be improved to capture this? If it cannot be addressed, perhaps the authors can add it to the discussion.\nResponse 10: We thank the reviewer for this insightful comment. The apparent “later onset” of Wnt signaling in the original Fig. 4c arises from a visual thresholding effect in the heatmap representation rather than from the model dynamics themselves. As shown in the revised Response Fig. R14, the simulated Wnt profiles increase continuously from the earliest time points across all densities (Fig. R14a). However, because the heatmaps use a perceptually nonlinear colormap, small early increases remain within the dark-blue range and only become visually apparent once they exceed a certain fraction of the full dynamic range. This creates the appearance of a delayed rise despite a smooth, monotonic increase in the underlying trajectories.\nTo directly compare model and experiment in the same radial region, Fig. R14 shows quantitative traces computed from bins 8–12, corresponding to 40–60 μm from the colony edge. These traces show that both the simulations (Fig. R14a) and the experimental data (Fig. R14b) rise from the outset. The heatmap-based visual delay is therefore a contrast-rendering artifact rather than a biological discrepancy.\nFigure R14: Simulated and experimental Wnt activation dynamics across colony densities.\n(a) Simulation kymographs and mean Wnt trajectories (bins 8–12; ~40–60 µm from the colony edge) show continuous early activation with density-dependent amplitude. (b) Experimental kymographs and bin-8–12 mean traces using the initial global percentile normalization similarly show early activation; any impression of delayed onset arises from colormap thresholding rather than the underlying data. (c) Revised experimental kymographs and bin-8–12 traces using baseline-scalar normalization (mean intensity at t = 0) better represent the true dynamic range across densities and align more directly with the simulation output.\nImportantly, we also identified that our original normalization approach contributed to this apparent mismatch. In the initial submission, we normalized experimental intensities using global 25th–75th percentiles computed across all colonies. This suppressed genuine baseline differences (particularly the higher starting Wnt levels present in high-density colonies) thereby exaggerating the visual thresholding effect. In the revised manuscript, we instead normalize each colony to its mean intensity at t = 0, which preserves biologically meaningful differences in initial signal levels and yields time-series traces that align closely with the model predictions (Fig. R14c).\nChanges to the manuscript:\nWe have revised Fig. 4f to include the baseline normalized kymographs shown in Fig. R14c.\n\nComment 11: As a separate point about the model, it is surprising that SOX2 is modeled essentially as “not GATA3 or BRA”, given that causally it is not true that SOX2 is directly downregulated by these factors (possibly the reverse) and that ultimately the authors emphasize SOX2 stability as the most important model feature. Is it fair to say SOX2 stability simply means degree of differentiation the way the model is set up? More importantly, if the authors truly think drugs directly affect SOX2 stability, it would be good if they would show this, similar to how they showed the density-dependent Wnt signaling.\nResponse 11: We appreciate the reviewer’s point and clarify that while the rule SOX2 = Cs·(c₁ + c₂·GATA3 + c₃·BRA) can be interpreted as a form of “not GATA3 or BRA,” our intent was to construct a minimal and interpretable model that reflects experimentally observed relationships. Specifically, SOX2 and BRA display mutually exclusive expression domains during early human gastrulation. In human embryos, SOX2 expression is maintained in epiblast-like regions, whereas BRA expression emerges in the primitive streak and nascent mesoderm, with minimal spatial overlap between the two domains, highlighting their mutually exclusive expression during germ layer segregation.6 Single-cell RNA-seq from gastrulating human embryos shows that as BRA expression emerges, SOX2 is sharply downregulated in the same cells, suggesting an antagonistic transition in cell state.7\nWe extend this reasoning to GATA3 as well. In 2D gastruloid systems, GATA3 acts as a rapidly induced commitment gene during BMP4-driven differentiation of human ESCs. Its expression follows bistable, switch-like dynamics regulated by a positive feedback loop involving both autoregulation and BMP4-induced signaling.8 Additionally, spatial analyses of 2D gastruloids show GATA3 is excluded from SOX2+ territories and enriched in mesoderm- and trophectoderm-like cell types, reinforcing mutual exclusivity with the pluripotency network.9\nOur model does not assert direct repression of SOX2 by GATA3 or BRA. Indeed, evidence suggests that SOX2 loss is a precondition for robust BRA expression and mesodermal differentiation, and that SOX2 itself may antagonize BRA transcriptional targets.10-12 Our model is consistent with this framework: the SOX2 stability rule serves as an abstraction of the pluripotency–lineage opposition, not as a mechanistic repression map. Finally, we fully agree that clarifying the causal control of SOX2 is important, and in future work, we plan to investigate this using human ESCs engineered with live SOX2 reporters.\n\nComment 12: The discussion illustrates the issue with thinking of all drugs as directly affecting SOX2 stability. It mentions in line 620 that ROCK inhibitors destabilize SOX2. However, what more likely happens is that they prevent tight junction formation, leading to increased BMP signaling, which is what actually destabilizes SOX2. In other words, there are many possible indirect of a drug that can affect differentiation.\nResponse 12: We agree that many compounds in our screen likely influence SOX2 levels indirectly, by modulating upstream signaling pathways or epithelial organization rather than targeting SOX2 directly. In particular, we appreciate the point that ROCK inhibitors may destabilize SOX2 not by acting on SOX2 itself, but by disrupting apical junctions, thereby enhancing endogenous BMP signaling. We have revised the manuscript to reflect this change.\nMore broadly, our model treats SOX2 stability as a parameter that integrates diverse molecular and structural influences on pluripotency maintenance. Rather than implying direct targeting of SOX2, our analysis highlights how drugs that shift the pluripotency–lineage balance are consistently captured by this parameter in the model.\nChanges to the manuscript:\n“In our dataset (Supplementary Table 4), compounds associated with increased SOX2 stability included valproic acid, a histone deacetylase (HDAC) inhibitor widely used to enhance pluripotency and reprogramming efficiency.71 In contrast, ROCK inhibitors such as hydroxyfasudil and ripasudil were among the strongest SOX2 destabilizers. While these effects may not stem from direct SOX2 targeting, ROCK inhibition disrupts apical junctions and epithelial architecture, indirectly enhancing BMP signaling, which promotes differentiation.72 Recent studies have further shown that SOX2 levels reconfigure Wnt/β-catenin binding, thereby controlling whether Wnt signals promote pluripotency or mesoderm.73 Together, these findings reinforce the idea that SOX2 stability encodes a pluripotency threshold, linking morphospace position to chromatin state, environmental responsiveness, and the broader logic of fate decisions during gastrulation.”\n\nComment 13: In lines 487-496, the authors ask if the mechanistic model can capture the phenotype of drug treatment. They infer the SOX2 stability from interpolation of simulation for different SOX2 stabilities on the morphospace. Thus, it is not surprising that if one provides SOX2 stability and density one ends up roughly in the same place, the only way I can see this going wrong is if the experimental density used deviates strongly from the mean density coordinate on the morphospace. I guess this is what the authors may mean in lines 548-553. To summarize, this seems somewhat like circular reasoning and the authors may wish to make clearer whether this test really proves anything. Perhaps this could involve merging figures 5 and 6. Can the authors discuss more why the model cannot capture the effect of isoxazole and longdaysin? Should it not predict reduced/increased SOX2 stability for a Wnt agonist/inhibor?\nResponse 13: We thank the reviewer for this insightful comment. We agree that if SOX2 stability and density are used as inputs, and SOX2 stability is derived from the embedding, then the resulting morphospace location will often match the original, especially if the density does not deviate significantly from the mean. However, the key utility of this approach is not in reconfirming expected matches, but in identifying cases where simulated outcomes diverge significantly from experimental results, thereby revealing perturbations to the underlying developmental logic.\nSpecifically, the morphogen spread model assumes fixed biophysical and regulatory rules for how BMP, Wnt, and Nodal diffuse and react (Fig. 4a), including BMP receptor localization, Wnt auto-activation, and density-dependent modulation of Wnt signaling. These assumptions are built into the reaction-diffusion equations and held constant across all simulated conditions. However, small molecules can violate these assumptions by directly targeting key nodes in the signaling architecture.\nThis is precisely what we observe with the small molecules, isoxazole and longdaysin, in which their predicted location in morphospace deviated strongly from the experimental embedding. In the case of isoxazole 9, the model predicts moderate BRA levels given the drug’s moderate cell density (mean cell count: 1020 cells), but the actual experimental colonies display strong BRA expression. This deviation signals a mechanism outside of cell density or SOX2 stability. Indeed, literature shows that isoxazole 9 potentiates Wnt/β-catenin signaling by covalently binding Axin1 and stabilizing β-catenin, promoting mesodermal gene expression (Jang et al., 2015). Conversely, longdaysin–treated colonies have a low cell density and as a result, the model predicts high BRA expression; however, experimental colonies show low BRA expression. Again, this is consistent with its known Wnt-inhibitory effects via disruption of LRP6 phosphorylation and reduced nuclear β-catenin accumulation (Xiong et al., 2019). In both cases, the drugs appear to modify the core rules of morphogen spread or cell fate determination, resulting in model deviation.\nThus, rather than being circular, this approach validates the model’s predictive power while revealing biologically meaningful deviations that highlight candidate mechanisms of action. We have updated the relevant text to better clarify this distinction.\nChanges to the manuscript:\n“Even deviations between mechanistic model predictions and experimental results revealed valuable insights. The morphogen spread model assumes fixed rules for BMP, Wnt, and Nodal signaling. However, certain compounds appear to override these assumptions. For example, isoxazole 9 (ISX-9)–treated colonies had relatively high cell counts (mean = 1,020 cells), which would predict moderate BRA expression in simulation; experimentally, however, they exhibited strong BRA expression (Supplementary Fig. 14a). This suggests an additional mechanism at play; indeed, isoxazole 9 has been shown to activate Wnt/β-catenin signaling by covalently binding Axin1 and stabilizing β-catenin,50 thereby enhancing mesodermal gene expression. Conversely, longdaysin-treated colonies had lower cell counts (mean = 892 cells), which would predict high BRA expression, but showed reduced BRA experimentally (Supplementary Fig. 14b), consistent with its reported Wnt-inhibitory action via disruption of LRP6 phosphorylation and reduced nuclear β-catenin accumulation.51 These discrepancies reveal when a perturbation alters the fundamental rules of morphogen signaling, enabling the model to identify biologically meaningful outliers rather than merely recapitulating the input.\n\nMinor comments:\nComment 14: The introduction line 60 states: “methods to assess the teratogenicity of compounds have not advanced beyond exposing model organisms”, but this is inaccurate, and the authors seem aware as they later cite two micropatterned approaches of teratogenicity testing. There are several other papers describing vitro approaches they may also wish to include,.e.g., PMID: 26387793, PMID: 31207406, and others.\nResponse 14: In the revised manuscript, we have removed this statement and now explicitly benchmark our model against existing human stem cell–based in vitro approaches. This comparison validates our method’s ability to capture teratogenic effects across a shared set of compounds and better positions our work within the context of ongoing developments in the field. PMID: 26387793 and PMID: 31207406 have been added as references in the Discussion section.\n\nComment 15: Line 81: “3D models such as embryoids, blastoids, and SEMs” – embryoid is a general term that includes blastoids and SEMs. Was something different meant?\nResponse 15: We agree that “embryoid” is a general term encompassing blastoids and SEMs. We have corrected the text to reference only blastoids and SEMs.\n\nComment 16: Line 85: “cost of altering the topology of the developmental tissue” – which tissue? It is true that the topology of the amnion is altered, but the epiblast geometry and topology is a disc as in vivo, which it actually is not some 3D models. It is a misunderstanding that 3D is always more similar to in vivo than 2D.\nResponse 16: We thank the reviewer for this clarification. We have revised the text to specify that the topological limitation applies primarily to extraembryonic tissues such as the amnion, not the epiblast, whose disc-like geometry is faithfully recapitulated in 2D systems. We also agree that 2D is not inherently less physiological than 3D and have updated the text to reflect this nuance.\n\nComment 17: Figure 1c shows classes of signal pattern perturbation but as far as I can tell two of these never occur: inverted and mixed, so why list them?\nResponse 17: We thank the reviewer for pointing this out. The inverted and mixed classes in Fig. 1c were intended as illustrative of possible signal perturbation modes, rather than empirically observed classes. However, we do observe a form of “mixed” patterning in cluster C12, which is characterized by a loss of radial symmetry and interspersed cell fates. To prevent confusion, we have updated the figure caption to clarify that these classes represent potential rather than strictly observed outcomes.\n\nComment 18: Clusters should be indicated in Fig.1e.\nResponse 18: Figure 1e was designed to present an unbiased visualization of all colony phenotypes prior to introducing any clustering or categorization, serving as a reference for the full dataset’s structure. We intentionally omitted cluster labels here to avoid prematurely implying discrete groupings before they are formally defined in Figure 1f. This approach allows readers to appreciate the diversity and distribution of patterning behaviors without presupposing outcomes from downstream analysis.\n\nComment 19: The authors use GELU for one network, RELU for another and make a point out of this in the main manuscript, however, it is not explained why. More generally data supporting the choice of network architecture is missing.\nResponse 19: In the original version, the two models used different activation functions (GELU and ReLU), though this reflected a historical implementation detail rather than a deliberate design choice. In the revised manuscript, we have standardized both models to use ReLU activations for consistency and clarity. The manuscript and code have been updated accordingly.\n\nComment 20: The shadows underneath the colony pictures throughout the manuscript, e.g., Fig.1g make the patterns harder to distinguish, perhaps a black outline would be better if contrast from the which background is needed.\nResponse 20: We thank the reviewer for pointing this out. In the revised manuscript, we have removed the colony shadows in the figures.\n\nComment 21: Line 345: “cell seeding follows a Poisson distribution”. That is only true for very sparse seeding and typically not for micropatterns which are densely seeded. The authors should show this or could simply remove it since it does not really matter, their only point is there is naturally occurring density variation.\nResponse 21: We thank the reviewer for this clarification. Our intent was simply to convey that natural density variation arises across colonies, which enables analysis. We have revised the text accordingly to avoid making an inaccurate statistical claim.\n\nComment 22: The regression of BRA peak vs. cell number in Supp. Fig. 5ac is quite nice. However, it would be good if in Fig. 1c not only the cluster averages but all datapoints were shown on the graph. Conversely, it would be helpful to see the cluster averages on the t-SNE plot.\nResponse 22: We thank the reviewer for this helpful suggestion. In the revised version (Fig. R15), we now display all individual data points colored by cell density alongside the regression of cluster averages. This update provides a more comprehensive view of the relationship between cell density and BRA peak position. We have not added cluster averages to the t-SNE plot to maintain visual clarity (and because this would require mixing the projected mean embedding with experimental data), but would be happy to do so if the reviewer feels it is important.\nFigure R15: Revised Supplementary Fig. 6c with additional data points shown.\n\nComment 23: In line 461-466, it is said density and SOX2 stability emerged as key parameters. How did they emerge? How were key parameters identified?\nResponse 23: We appreciate the reviewer’s request for clarification. To address this, we performed an expanded analysis where all parameters in the cell fate model were systematically varied across biologically plausible ranges. As shown in Fig. R16, we defined each regulatory term mathematically and scanned each parameter at multiple cell densities (Fig. R17).\nFigure R16: Expanded parameterization of mechanistic cell fate model to evaluate regulatory contributions. (a) Equations defining GATA3, BRA, and SOX2 as functions of morphogen inputs and regulatory coefficients. (b) Table of model parameters with biological interpretation. (c-d) Simulated colony embeddings as a function of cell density and SOX2 stability (Cₛ), showing orthogonal axes of variation.\nFigure R17: Systematic variation of all model parameters across morphospace. Simulated colony embeddings in morphospace while varying each parameter across a biologically plausible range (0.25–1.75) and at multiple cell densities. Among all parameters tested, only Cₛ (Fig. R4) and cell density generate trajectories that align with the major axes of developmental variation, highlighting their interpretability for linking drug-induced shifts to developmental mechanisms.\nThis parameter scan revealed that only cell density and SOX2 stability (Cₛ) produced trajectories that aligned with major axes of morphospace variation and span the entire canonical patterning region. We have added these figures to the revised manuscript and revised the text accordingly.\nChanges to the manuscript:\nAdded Figure R16-17 as the revised Supplementary Fig. 11-12.\n“To assess whether other coefficients could generate comparable trajectories, we systematically swept all model parameters over biologically plausible ranges; only density and SOX2 stability consistently produced extended, coherent axes in morphospace (Supplementary Figs. 11-12).\n\nComment 24: Line 566: “animal models, which are limited by species-specific differences and dose-dependent lethality”. This appears to suggest that dose-dependent effects are not an issue in this work but they are and the authors acknowledge this later in the discussion (everything done at 10uM), which is good, so they could consider rephrasing line 566.\nResponse 24: We thank the reviewer for this comment. Our original statement was meant to highlight the challenge of identifying embryo-specific teratogenic effects in animal models, where high doses may be lethal to the mother before revealing developmental toxicity. However, we agree this could be misread as implying that FATE-MAP avoids all issues related to dose-dependent effects, which it does not. We have revised the sentence to focus instead on practical limitations of animal models, such as species differences and the need for in vivo dosing and animal care.\n\nReviewer 3 (Remarks to the Author)\nIn this manuscript, Rufo et al., develop an experimental and analytical screening platform, FATE-MAP, to assess phenotypic outcomes in human 2D stem cell models of early embryogenesis (human 2D gastruloids) treated with hundreds of chemical perturbations. FATE-MAP includes a deep learning framework to predict the effects of any molecule on 2D gastruloids development based on the molecule’s chemical structures. The authors first characterized “failure modes”, phenotypes that emerge when the 2D gastruloids are not properly patterned and then analyze over two thousand structures subjected to roughly 200 perturbations, including clinically relevant teratogens, eventually identifying several molecules that results in aberrant developmental programs. The authors also performed some validation experiments using the zebrafish model to corroborate the predicted molecule-induced phenotypic outcomes. They then model morphogen dynamics in 2D gastruloids, an approach that revealed how cell density and expression stability of the transcription factor SOX2 play a major role in the correct patterning of 2D gastruloids. Finally, they attempt to derive mechanistic principles of failure modes from their deep learning approach, implementing a mechanistic model to gain interpretability and biological meaning.\n\nOverall, this reviewer thinks the developed methodology is important and timely, given the little knowledge we have about the effect of environmental offences on embryonic development, especially humans. Moreover, cheap and accessible high-throughput screening opportunities are being developed, and I think this method and its future implementations (like others of its kind) are in need to the field. Science and experiments are solid and well executed. However, this reviewer thinks the biological findings are not groundbreaking, especially given the opportunity provided by the method. The manuscript is very well written and easy to understand for a broad audience. I generally support its publication at Nature Communications, pending some points that would need to be addressed.\n\nMajor comments:\n1) List of selected molecules:\nComment 1: It would be important if the author could comment more on the list of 210 drugs tested. How were they selected? What is known about their effect in other species? Are they already classified or stratified based on their chemical structures?\nResponse 1: We thank the reviewer for raising this point. The 210 compounds tested were selected as a subset from the commercially available Stem Cell Signaling Compound Library (Selleck Chemicals, Cat #: L2100), with the goal of probing perturbations to major developmental signaling pathways. Since the initial aim was not specifically teratogen screening, many of these compounds do not have annotated teratogenicity data in other species.\n\nComment 2: This becomes more relevant in the efforts to build a classifier based on chemical structures. Are all molecules causing one phenotypic cluster similar/dissimilar in structure? Do they share any chemical group/property (e.g. structure, domains, etc..)? Figure S3a looks at the similarity between seen and unseen drugs in the classifier but does not compare the general chemical structures in the drug list.\nResponse 2: To address the reviewer’s suggestion regarding structural stratification, we clustered all experimentally screened compounds based on their chemical structure. We computed MACCS chemical fingerprints for each compound, projected them into 2D space using t-SNE, and applied KMeans clustering to identify major chemical clusters. We also extracted the most common Bemis–Murcko scaffold within each cluster and visualized representative drugs that share these scaffolds. These results are now summarized in Fig. R18.\nFigure R18: Structure-Based Clustering of Experimentally Screened Compounds. (a) t-SNE projection of MACCS chemical fingerprints for all experimentally screened compounds, colored by KMeans-derived chemical cluster. (b–d) Representative drugs from Clusters 6, 7, and 8, each shown alongside the most common Bemis–Murcko scaffold for the cluster.\n\nIn panel a, all experimental compounds are visualized in t-SNE space and colored by their assigned chemical cluster. Panels b–d highlight clusters 6, 7, and 8, respectively. For each cluster, we show the most common scaffold (left) and representative drugs that contain this scaffold (right). For example, Cluster 6 includes PD98059, 2-D08, and Tangeretin, all of which share a core aromatic ether scaffold. Cluster 7 contains sulfonamide-linked tricyclic inhibitors such as BP-1-102 and SH-4-54, which target the STAT3 pathway. Cluster 8 includes PD169316 and SB203580, which share a pyrazolopyridine core and are known p38 MAPK inhibitors. This structural clustering highlights the presence of chemically and functionally coherent drug groupings within our experimental library.\n\nComment 3: What is the inter cluster vs intra cluster chemical structure similarity? Do the tested molecules cluster based on their chemical structure?\nResponse 3: To further assess whether chemical similarity aligns with phenotypic outcomes, we compared the structure-based groupings from Fig. R18 to the morphology-derived clusters obtained from our morphospace embedding.\nFigure R19: Comparison of Phenotypic and Structural Clustering of Experimental Compounds. (a) Experimental compound centroids (n=210) visualized in the phenotypic embedding and colored by morphology-derived cluster assignment. (b) The same compounds colored by chemical structure clusters (derived from MACCS fingerprints and KMeans clustering). The structural clusters show extensive overlap across phenotypic space, indicating weak correspondence between chemical similarity and phenotypic behavior.\n\nThis comparison is shown in Fig. R19. While the left panel shows the original phenotypic clusters, the right panel displays the same compounds colored by their chemical cluster. The chemical clusters are broadly dispersed across morphospace, with substantial overlap and little alignment to the phenotype-defined groupings. This discrepancy suggests that structurally similar molecules do not necessarily induce similar morphological effects, highlighting the challenges of relying solely on chemical structure for phenotype prediction. Together, these results underscore the added value of phenotypic profiling in stratifying small molecule effects beyond what structural similarity alone can capture.\nChanges to the manuscript:\nFigures R18 and R19 have been added to the revised Supplementary Fig. 5.\n“To explore this further, we clustered all experimentally screened compounds based on structural similarity and compared these groupings to morphology-defined phenotypic clusters (Supplementary Fig. 5). Despite the presence of chemically coherent scaffold-based clusters, we observed minimal spatial overlap between structure- and phenotype-defined groupings. This suggests that even compounds with similar scaffolds can produce distinct developmental effects, reinforcing the need for neural network–based approaches to uncover structure–function relationships that scaffold similarity alone may obscure.”\n\n2) Screening design:\nComment 4: The authors should test different concentrations for the drugs they use. This reviewer understands well that this won’t be possible for all compounds, but using a fixed concentration for all drugs (10uM) limits the real assessment of their effects given that not every molecule has the same stability, processivity, and activity. The authors should titer different concentrations to understand dose-specific responses for a shortlisted number of compounds (e.g. the ones resulting in phenotypes belonging to C1, C11, and/or C12).\nResponse 4: We agree that assessing compounds at multiple concentrations provides critical insight into dose-specific effects, cytotoxicity, and mechanistic thresholds. This approach enables the identification of metrics such as the lowest-observed-adverse-effect concentration (LOAEC) and offers a more nuanced characterization of teratogenic potential. While our initial screen was conducted at a fixed 10 µM dose to ensure coverage of a wide compound library, we recognize the importance of titration experiments for selected hits.\nTo address this, we followed up with dose–response testing for a subset of retinoid compounds identified in our initial screen. Retinoids in our initial screen, including isotretinoin, bexarotene, TTNPB, and tretinoin (Fig. R20a) reproducibly induced the morphological C11 phenotype, marked by overgrowth and suppression of BRA expression. Notably, tretinoin also induced a secondary phenotype associated with overexpression of GATA3, aligning with C12 (Fig. R20b).\nIn addition to these annotated retinoids, we included an unannotated compound (IB-2744; Specs Cat. #AK-968/15360521) in our follow-up experiments. This allowed us to demonstrate that the same titration workflow can also reveal mechanistic similarity for compounds lacking prior annotation, based solely on their morphological phenotype (Fig. R20e).\nTo explore the dose dependence of these effects, we titrated tretinoin and isotretinoin across 10 µM, 1 µM, and 0.1 µM concentrations (Fig. R20c–d). Both compounds recapitulated the original failure mode phenotypes observed at 10 µM, with dose-specific trends.\nFigure R20: Dose–Response Effects of Retinoids on Patterning and Morphology. (a) Retinoid compounds from the initial screen, including isotretinoin, bexarotene, TTNPB, tretinoin, and AM580, visualized in morphospace. (b) Representative colony phenotypes induced by each compound at 10 µM. (c–d) Follow-up titration experiments for tretinoin and isotretinoin at 10 µM, 1 µM, and 0.1 µM. (e) Titration of an unannotated compound (IB-2744), demonstrating the platform’s ability to identify compounds with phenotypes suggestive of shared mechanisms of action.\n\nFor tretinoin, a clear dose-dependent shift in phenotype was apparent. At 1 µM, most colonies continued to exhibit disrupted patterning and GATA3 overexpression, though with some colonies beginning to show patches of BRA expression. By 0.1 µM, a substantial number of colonies showed restoration of canonical BRA ring formation, suggesting a concentration threshold below which teratogenic disruption is mitigated.\nIn contrast, isotretinoin exhibited a more consistently potent effect across all doses. At 10 µM and 1 µM, colonies showed complete absence of BRA expression, with dense overgrown SOX2+ cores. At 0.1 µM, some of the colonies began to express BRA; however, none of the colonies recapitulated canonical patterning, suggesting isotretinoin’s teratogenic effects are more potent and less reversible within the concentration range tested. Interestingly, the unannotated compound IB-2744 produced a highly similar dose–response profile, suggesting that it may act through a related retinoid-like mechanism despite lacking prior annotation.\nThese findings suggest that isotretinoin and tretinoin may exert their teratogenic effects through distinct mechanisms, even though isotretinoin is a known metabolic precursor to tretinoin in vivo. While some conversion to tretinoin may occur in our stem cell platform, recent studies indicate that isotretinoin itself can directly perturb key developmental signaling pathways in human pluripotent stem cells, including Wnt/β-catenin, TGF-β/Nodal, Notch, and Hedgehog. These pathways are critical for early germ layer specification, and isotretinoin has been shown to suppress BRA expression and inhibit mesoderm formation at nanomolar concentrations. Thus, our observation of strong and persistent morphological disruption by isotretinoin likely reflects a combination of direct and metabolite-driven activity, and points to a distinct teratogenic mechanism that may not be fully recapitulated by tretinoin alone. Our platform offers a unique opportunity to dissect such isoform-specific effects in a controlled, human stem cell context.\nIn future work, we aim to scale our assay to include systematic dose–response curves across a broader range of drug classes. Additionally, incorporating metabolic activation systems, such as enzyme mixes commonly used in genotoxicity assays like the Ames test, could provide further insight into the effects of metabolic conversion. Nevertheless, these proof-of-concept experiments demonstrate the utility of our platform for quantitative phenotypic profiling and highlight the value of dose stratification in characterizing compound-specific effects.\nChanges to the manuscript:\nAdded Figure R20 as the revised Supplementary Fig. 15.\n“As a proof-of-principle, we performed titration experiments for two representative cluster 11 compounds (tretinoin and isotretinoin) and for an unannotated compound (IB-2744), observing clear dose-dependent phenotypes that mirrored their activity at 10 µM (Supplementary Fig. 15). These results demonstrate that our platform can resolve graded responses across a concentration range and highlight its potential for identifying lowest-observed-adverse-effect thresholds and for detecting mechanistic similarity among previously uncharacterized compounds.”\n\n3) Comparison with other screening approaches, deep learning methods and relevant literature:\nComment 5: A similar teratogen screen has been performed in 3D human gastruloids, which offer a more physiological and morphologically relevant model (PMID: 34425190). This work uses a way shorter list of compounds but various drug concentrations and combines morphometric analysis and gene expression patterning. Given that many molecules are present in both studies, the authors should compare their observed phenotypes in 2D in relation to the one in 3D. Also, if some molecules that have been used in PMID: 34425190 are not present in the authors’ list, they should predict the phenotype using FATE-MAP and compare the prediction with the observed phenotype in 3D.\nResponse 5: We thank the reviewer for this helpful suggestion. In the revised manuscript, we now benchmark FATE-MAP structure-to-phenotype predictions against a curated panel of 33 FDA-approved compounds with well-characterized teratogenicity outcomes, supported by either human clinical evidence or standardized embryo–fetal development (EFD) studies in rats and rabbits (Response Fig. R21b; revised Fig. 3b). Importantly, all 33 compounds have also been evaluated in two stem-cell–based assays designed for teratogenicity testing: TeraTox (3D human iPSC multilineage differentiation) and the mEST (mouse embryonic stem cell test) assay.\nOur structure-to-phenotype model first generates a 150-dimensional predicted phenotype vector for each compound and then projects this vector into the morphospace, where cluster identity is assigned using the same watershed-derived decision boundaries defined in the original clustering (Response Fig. R21c). Compounds whose predicted embeddings fall into the failure-mode clusters C1, C11, or C12 are labeled teratogenic, while those in canonical patterning clusters C2–C10 are labeled non-teratogenic.\nOf the 33 benchmark compounds, only four were present in our experimental dataset; the remaining 29 were predicted fully in silico. Across this set, FATE-MAP achieved 72.7% accuracy, outperforming mEST (51.5%) and matching TeraTox (72.7%), while exhibiting higher balanced accuracy than TeraTox (74.8% vs. 66.3%; response Fig. R21d)). These results demonstrate that FATE-MAP can serve as a structurally grounded, scalable framework for in silico teratogenicity screening.\nFigure R21 (Revised Figure 3 in manuscript). FATE-MAP predicts teratogenic risk from chemical structure and identifies novel human-relevant teratogens. (a) Representative 2D gastruloids treated with known teratogens isotretinoin, TTNPB, and bexarotene, all converging on cluster 11 and exhibiting a shared failure mode marked by BRA loss. (b) A benchmark panel of 33 compounds with known teratogenicity outcomes (positive and negative), previously evaluated in stem cell-based assays of developmental toxicity (TeraTox and mEST; refs), was used to assess predictive performance. Predicted phenotypes were embedded into morphospace and classified using the fixed watershed-defined cluster boundaries derived from the experimental dataset. (c) Morphospace-based classification of the benchmark set: embeddings falling within failure-mode clusters (C1, C11, C12) are assigned positive teratogen status, whereas embeddings within the canonical patterning region (C2–C10) are assigned non-teratogenic status. Points are labeled as true positive, true negative, false positive, or false negative based on known outcomes. (d) Confusion matrices and comparative performance metrics for FATE-MAP, TeraTox, and mEST on the benchmark panel.(e) Two blinded compounds without prior human pregnancy data—ulixertinib and nirogacestat—were predicted and experimentally confirmed to fall within cluster 11. (f) Zebrafish teratogenicity assays show that both compounds, like the positive control isotretinoin, induce qualitative and quantitative defects consistent with failure modes in embryonic development, including abnormal body curvature (BA), yolk edema (YE), necrosis (NC), craniofacial edema (CE), scoliosis (SC), snout and jaw defects (SJ), reduced lateral heart area (LH), and decreased lateral body length (BL), as indicated by the color-coded bars.\nRegarding the teratogen screen performed in 3D human gastruloids (PMID: 34425190) referenced by the reviewer, the only overlapping experimental compounds were retinoic acid and valproic acid. Both approaches (2D and 3D) detected clear morphological changes to retinoic acid treated gastruloids (although we only test at 10 µM and the 3D approach noted morphological changes at 333 µM). In our approach, we did not note any morphological changes to valproic acid treated 2D gastruloid at 10 µM. The 3D approach reports minimal effect at 4 µM and morphological effects, gene expression effects, and size reduction at 333 µM. For comprehensive comparison, we used our FATE-MAP structure-to-phenotype model to predict the phenotype of all 7 compounds tested (Fig. R22). The model correctly predicted the teratogen status of retinoic acid and bosentan as well as the non-teratogenic status of penicillin and ibuprofen (referred to as “less teratogenic” in PMID: 34425190); however, our model incorrectly classified valproic acid, thalidomide, and phenytoin. We note the valproic acid and thalidomide are both included in our 33-compound validation set.\nFigure R22: Benchmarking Structure-Based Predictions Against 3D Gastruloid Assay Outcomes. Predicted phenotype embeddings for 7 compounds previously evaluated in a 3D human gastruloid teratogenicity screen (PMID: 34425190), overlaid on the experimental morphospace. The background shows the fixed cluster boundaries derived from experimental data (as in Fig. R19c), with light grey dots representing the distribution of all experimentally screened compounds. While FATE-MAP correctly identifies the teratogenicity of retinoic acid and bosentan, and the non-teratogenicity of penicillin and ibuprofen, it misclassifies valproic acid, phenytoin, and thalidomide.\nChanges to the manuscript:\nResponse Fig. 21 has been added as the revised Fig. 3 to include benchmarking of FATE-MAP predictions against 33 FDA-approved compounds with experimentally established teratogenicity outcomes in stem-cell-based assays.\n\nComment 6: Another relevant study has been recently published (PMID: 40245869) combining morphological variation in 3D mouse trunk-like-structures with gene expression profiles to train a classifier able to predict phenotypic outcomes from early morphological features. Even if this reviewer acknowledges that the experimental design is different (this work uses longitudinal imaging over time to identify predictive features), the authors should cite this work and comment on similarities and differences in their approaches.\nResponse 6: The work by Chhabra et al. (PMID: 40245869) represents a powerful example of integrating live imaging, transcriptomics, and machine learning to predict morphogenetic outcomes in 3D trunk-like structures derived from mouse stem cells. While their study focuses on early morphological predictors of axial elongation and segmentation in 3D contexts, our work centers on linking chemical structure to phenotypic outcomes in a human 2D gastruloid model, with an emphasis on deciphering morphogen-driven failure modes during early germ layer patterning. Despite these differences in biological system and input features, both approaches share a common goal of using computational models to derive mechanistic insight and improve predictive inference in developmental systems. In particular, their use of multimodal data integration and explainable machine learning provides a valuable model for future work aiming to relate latent parameters (such as SOX2 stability in our model) to underlying molecular and morphological processes. We have now cited this study in the Discussion.\nChanges to the manuscript:\n“A recent study integrating time-lapse imaging of morphology with single-cell transcriptomics in 3D mouse trunk-like structures used machine learning to identify early features predictive of fate and found that metabolic state plays a key role in regulating differentiation outcomes.74 Similar multimodal approaches could be applied to investigate the molecular determinants of SOX2 stability and its role in patterning outcomes in human gastruloids.”\n\nComment 7: A recent study in engineered mouse 3D gastruloids (PMID: 39358450) provides temporal recording measurements of signaling patterns and specifically studies the emergence of Wnt domains and its relationship to Nodal. The authors should comment and relate their results to the ones reported in this work, since this comparison could provide species-specific and shared features between mouse and human.\nResponse 7: We thank the reviewer for highlighting this study of Wnt and Nodal patterning dynamics in mouse 3D gastruloids (PMID: 39358450). While our approach employs 2D human gastruloids and theirs uses 3D mouse aggregates, both systems interrogate early symmetry breaking using live-cell Wnt reporters and experimentally probe BMP, Wnt, and Nodal signaling pathways.\nThere are several notable points of alignment. Both studies show that Nodal signaling is required for proper morphogenesis: in the mouse gastruloids, inhibition of Nodal abolishes elongation and symmetry breaking, consistent with existing literature that shows that Nodal inhibition prevents BRA expression in 2D gastruloids. Furthermore, both models explore BMP–Wnt–Nodal crosstalk and find that modulation of one pathway can influence downstream responses in the others. In particular, the mouse system shows that early Activin A (Nodal) stimulation can dampen Wnt signaling later in development, echoing the reciprocal feedbacks modeled in our system.\nHowever, the mechanistic basis of symmetry breaking differs substantially. The mouse gastruloid model emphasizes cell sorting and spatial segregation as primary drivers of axis formation, whereas our system exhibits robust symmetry breaking and patterning based on dynamic, spatially coupled signaling waves and inferred reaction–diffusion–like behavior. Nevertheless, we now include a sentence in the Discussion acknowledging this complementary approach and emphasizing the broader utility of live-cell signaling reporters in dissecting human and mouse gastrulation models.\nChanges to the manuscript:\n“Recent work in engineered mouse 3D gastruloids has similarly used live-cell signaling reporters to dissect the timing and interplay of BMP, Wnt, and Nodal signaling during early axis formation,69 highlighting the utility of these tools across species and model systems.”\n\n4) Use of the term “morphospace”:\nComment 8: The authors use the term “morphospace” to define the landscape of phenotypic outcomes arising in the treated 2D micropatterns. Even if several colleagues in the field have argued about the morphological organization of micropatterns, including the 2.5D nature of the model with a z dimension also being relevant for proper cellular organization (which the authors do not use in their analysis), this reviewer thinks that 2D gastruloids lack morphological features of an embryo. The model is great to study tissue patterning, cell interactions and the effect of molecules on these features, but not so great to assess the acquisition of complex embryonic morphologies. I would recommend converting the wording “morphospace” into a “patterning space” (or similar) since this is what the model (and the readout used by the authors) can accurately inform.\nResponse 8: We thank the reviewer for this thoughtful and constructive suggestion. We agree that 2D gastruloids do not fully capture the complex morphological architecture of the embryo and that our model is best suited for studying fate patterning and signaling interactions rather than three-dimensional morphogenesis. Our use of the term morphospace was not intended to imply embryonic-level morphological complexity, but rather follows prior usage in the field (ref. 16), where morphospace refers to a configuration space in which individual systems are positioned based on defined phenotypic properties. In our case, these properties correspond to spatial patterning features in 2D micropatterned colonies.\nTo avoid confusion, we now include a clarifying statement when the term is first introduced in the manuscript, emphasizing that our use of morphospace refers specifically to 2D patterning phenotypes in micropatterned colonies and is not meant to suggest a direct morphological analogy to embryonic structure.\nChanges to the manuscript:\n“Central to this approach is the construction of a morphospace16,17—a low-dimensional, quantitative map of morphological phenotypes—that captures both canonical and disrupted trajectories of human gastrulation. Here, we use morphospace specifically to describe the diversity of 2D patterning outcomes in micropatterned stem cell colonies, not to imply full embryonic morphology or organization.”\n\n5) Predictive power:\nComment 9: The screening experiments performed by the authors offer the unique opportunity to inspect what was not “well predictive” (e.g inter cluster violations between prediction and observation) which could potentially result in unexpected, interesting observations/exceptions to follow up with. Did the authors observe any of these violations in their data?\nResponse 9: One of the central findings of our study is that most compounds were accurately predicted by our reaction-diffusion–based morphogen spread model, with simulated phenotypes closely matching experimentally observed fate patterns when informed by cell density and inferred SOX2 stability. However, the most biologically revealing cases were those where prediction and experimental data diverged. These deviations were not random; rather, they reflected meaningful violations of the model’s core assumptions, providing mechanistic insights into how certain compounds disrupt normal signaling dynamics.\nAs described in Fig. 4a, our model encodes fixed biophysical and regulatory rules for BMP, Wnt, and Nodal signaling, including edge receptor localization, Wnt auto-activation, and cell density–dependent modulation of morphogen availability. These parameters remain constant across all simulations. Thus, any discrepancy between predicted and experimental embeddings indicates that a compound has altered the underlying signaling logic (i.e., targeting morphogen production, degradation, diffusion, or pathway response in a way not captured by the model).\nSupplementary Fig. 9 (Response Fig. R23) highlights two such examples. In the case of isoxazole 9 (ISX-9), the predicted BRA pattern (based on low cell density and moderate SOX2 destabilization) was weaker than observed, whereas the experimental data showed robust BRA induction (Supplementary Fig. 9a). This discrepancy aligns with prior reports showing that ISX-9 directly stabilizes β-catenin by covalently binding Axin1, thereby enhancing Wnt signaling and mesodermal gene expression.13 Conversely, longdaysin was predicted to induce strong BRA expression due to its low-density profile, but the actual colonies lacked BRA (Supplementary Fig. 9b). This too is consistent with literature showing that longdaysin inhibits Wnt signaling by blocking LRP6 phosphorylation and β-catenin nuclear translocation.14\nFigure R23 (Supplementary Fig. 9): Deviations between model predictions and experimental phenotypes reveal mechanistic insights into drug action. (a) Isoxazole 9 (ISX-9)–treated colonies had high cell counts (mean = 1,020), which the morphogen model predicted would yield moderate BRA expression; however, experimental data showed strong BRA induction, consistent with reports that ISX-9 activates Wnt/β-catenin signaling via Axin targeting. (b) Longdaysin–treated colonies had relatively low cell counts (mean = 892), which would predict high BRA levels in simulation, but instead showed reduced BRA experimentally, aligning with evidence that Longdaysin inhibits Wnt/β-catenin signaling. \nThese examples underscore how deviations from model predictions serve as a diagnostic tool: identifying compounds that rewire the signaling architecture in unanticipated ways. We have updated the text in the revised manuscript to emphasize that such outliers not only validate the utility of the model but also point to novel mechanisms of action for individual compounds.\n\nComment 10: In line with what this reviewer commented on before, are there any features of the molecules’ chemical structures with predictive power? Finding chemical similarities that results in similar phenotypic outcomes might be a key and novel aspect of the study.\nResponse 10: To investigate whether shared chemical features underlie similar morphological effects, we examined associations between chemical and morphological clusters. As shown in Fig. R19b, we identified a modest enrichment between morphological cluster 11 (characterized by failure to induce BRA expression) and chemical cluster 9, with 12 of 36 compounds in this morphological cluster belonging to chemical cluster 9 (p = 0.0016). This was the only chemical cluster that showed significant enrichment (p<0.05) among the morphological clusters. However, further analysis of the 12 enriched compounds revealed that no generic Bemis–Murcko scaffold was shared across more than one compound. This suggests that even when abstracting chemical structures to core scaffolds, the enriched compounds are structurally diverse.\nThis result underscores the challenge of inferring predictive structure–activity relationships solely from scaffold similarity in phenotypic screening contexts. While we observe modest enrichment at the chemical cluster level, the absence of recurring core scaffolds across enriched compounds emphasizes the complexity of morphological outcomes and suggests that phenotypic convergence can arise from structurally diverse agents, possibly via shared functional or mechanistic motifs not easily captured through traditional scaffold analysis.\n\n6) Cell death vs differentiation defects\nComment 11: Cell death measurements are not included in the experiments. For example, C11 (failure to express BRA) could result from cells failing to differentiate (as hypothesized by the authors) or compound lethality (which could still result in a similar phenotype). To pinpoint a real example, there is a steep difference between Cerdulatinib treatment (looks more like a lethal phenotype) and D4476, but both end up in C11. The authors should classify C11 further and pinpoint lethality vs differentiation failure phenotypes. Ideally, for a set of compounds, a live/dead dye should be used to determine the cause of the observed failure mode. This is extremely relevant to correctly interpret the phenotypic outcome and failure mode.\nResponse 11: We thank the reviewer for raising this important point regarding the distinction between cytotoxicity and true differentiation failure, particularly within cluster 11. We fully agree that determining whether BRA loss reflects impaired fate induction or reduced viability is essential for correctly interpreting failure modes.\nFigure R24: Density-based toxicity threshold. (a) Mean nuclei per colony for each BMP4+drug condition. Toxicity threshold (red, ≤50% mTeSR1 mean) and hyperproliferation threshold (green, >mTeSR1 mean). (b) Toxic compounds map predominantly to failure modes C1, C11, and C12; insets show representative colonies. The embedding illustrates that overt cytotoxic collapse can be distinguished from density-dependent modulation within the canonical patterning region.\nAlthough we did not perform live/dead staining in the initial screen, our assay incorporates a built-in, physiologically relevant measure of viability by quantifying every nucleus in every colony. Because all wells are seeded at the same initial density (50,000 cells per well), large reductions in total nuclear count directly reflect compound-induced toxicity rather than density-dependent variation or growth arrest. Using this approach, we defined an objective, control-anchored toxicity threshold, classifying treatments whose mean nuclear count fell to ≤50% of the mTeSR-only control as toxic. As shown in Fig. R24a, only a small subset of compounds fall below this threshold, while several others exceed the mTeSR control mean and form a hyperproliferative group.\nMapping these conditions into morphospace revealed that 96.1% of colonies exposed to toxic compounds reside within failure-mode regions (C1, C11, or C12; Fig. R22b). This analysis allows us to distinguish lethal phenotypes from patterning failures within the same cluster, such as the different phenotypes pointed out by the reviewer between Cerdulatinib and D4476-treated colonies. Cerdulatinib-treated colonies fall far below the toxicity threshold and exhibit sparse, collapsed phenotypes consistent with pronounced cell loss, whereas D4476-treated colonies retain near-normal nuclear counts yet fail to induce BRA. Although both map to C11, their underlying causes differ: cerdulatinib produces a toxicity-driven collapse, while D4476 produces a bona fide patterning failure with preserved viability. Thus, within the existing dataset, we can already distinguish cytotoxic and non-cytotoxic subtypes of C11 based on quantitative viability metrics.\nWe have revised the Results section to clarify how toxicity-associated subclusters were identified and how these data distinguish lethal phenotypes from differentiation failures.\nChanges to the manuscript:\nAdded Fig. R24 as Supplementary Fig. 3.\n“Because failure modes can arise from either differentiation failure or compound-induced cytotoxicity, we quantified total nuclear counts for every colony and identified a small subset of drugs that reduce viable cell number below 50% of the mTeSR-only control (Supplementary Fig. 3). This analysis separates genuine patterning failures from cytotoxicity-driven collapse within the failure-mode regions.”\n\nComment 12: Cell density vs BRA expression: this relationship is established at the endpoint differentiation. It is known that exit from pluripotency to the PS is accompanied by cell death (and therefore might result in lower density). Is this hypothesis tested against the initial cell density, or the observed density reduction could be a consequence of enhanced differentiation?\nResponse 12: We thank the reviewer for raising an important point regarding the potential role of differentiation-associated cell death in shaping endpoint cell density and BRA expression. While it is true that exit from pluripotency can involve apoptotic events, we believe this is not the primary cause of the inverse relationship we observe between cell density and BRA expression.\nTo directly disentangle this, we performed live-cell time-lapse imaging of endogenous Wnt signaling in colonies seeded at deliberately different initial cell densities (Fig. R25; Fig. 4e–f in the manuscript). Since these colonies began with controlled and measured cell counts, this allowed us to assess the impact of initial density (rather than differentiation-induced cell loss) on BRA induction. In these experiments, lower-density colonies consistently exhibited earlier and more extensive propagation of the Wnt signaling wave, followed by stronger BRA expression at the endpoint (Fig. R25g). Importantly, these colonies maintained their relative density rankings throughout the experiment.\nThus, while we acknowledge that some degree of cell death may accompany differentiation, our results indicate that the observed density-dependent effects on BRA are largely attributable to initial cell seeding differences rather than a consequence of mesoderm formation. These findings further support our modeling assumption that Wnt signaling is modulated by cell density.\nFigure R25 (Fig. 5e-f in manuscript): Experimental validation of the impact of cell density on Wnt signaling and BRA expression. (e) Experimental investigation of density-dependent Wnt dynamics using CRISPR-tagged β-catenin (β-cat) hESC lines. Live-cell time-lapse imaging over 48 hours shows the effect of cell density on Wnt signaling. Scale bars = 100 μm. (f) Mean kymographs of Wnt signaling activity across 8 colonies per cell density range, demonstrating that lower cell densities result in a faster and broader spread of the Wnt signal. (g) Quantification of BRA staining intensity profiles across different densities, taken from the same colonies used to generate the kymographs. The area under the curve (AUC) analysis confirms increased total BRA signal at lower densities, validating our model’s prediction of density-induced variations in mesoderm dynamics.\n\n7) FDA drugs effect prediction:\nComment 13: This reviewer understands well the challenges associated with implementing the Pregnancy and Lactation Labeling Rule (PLLR) for FDA approved drugs given their non categorical classification, but the categorization used by the authors dates back from 1970s and lack rigorous testing. As specified by the authority and clinicians, this categorization is not a “grading system of risk” based on the actual severity of the teratogen, but rather a classification based on available data at the time. The authors should compare/consider including information from PLLR in their analysis. This could even improve their results since the current classification might suffer from lack of extensive testing (see the next two points regarding the A, B and X graded molecules).\nResponse 13: We thank the reviewer for raising this important point. As noted in our response to Comment 5, we have removed the PLLR- and pregnancy-category–based validation from the manuscript because these classifications are not mechanistic and were never intended to function as quantitative measures of teratogenic risk. Instead, we now validate FATE-MAP using a curated benchmark of 33 FDA-approved compounds with well-characterized teratogenicity outcomes based on human clinical evidence or standardized embryo–fetal development (EFD) studies, all of which have also been evaluated in stem cell–based teratogenicity assays (TeraTox and mEST).\nFor each compound, the structure-to-phenotype model generates a predicted phenotype vector that is embedded into the experimental morphospace, and classification is assigned using the fixed watershed-derived cluster boundaries (Fig. R21). Compounds mapping to failure-mode clusters (C1, C11, C12) are labeled teratogenic, while those in canonical clusters (C2–C10) are labeled non-teratogenic. Across the full 33-compound panel (only four of which appear in our experimental dataset) FATE-MAP achieves 72.7% accuracy and higher balanced accuracy than TeraTox (74.8% vs. 66.3%), providing a more rigorous and mechanistically grounded evaluation than PLLR-based labels.\n\nComment 14: On a similar line, the authors should show the embedding (or the statistical distance to clusters) for all FDA drugs, including class A and B compounds (not shown in the current version).\nResponse 14: We appreciate the reviewer’s suggestion. In the revised manuscript, Figure 3 (Fig. R21c) now displays the full morphospace embedding for all 33 validation compounds, including both non-teratogenic agents (corresponding to the previously categorized class A and B compounds) and known teratogens. By embedding all compounds directly into the experimental morphospace, rather than showing only a subset, we provide a complete visualization of where each drug falls relative to canonical and failure-mode clusters, addressing the reviewer’s request for comprehensive representation of the validation set.\n\nComment 15: Where are all the other class X compounds ending up (71 out of 98 class X compounds are not shown). Are they part of the other aberrant clusters or do they end up in the canonical “morphospace”?\nResponse 15: We thank the reviewer for this question. In the initial submission, the predicted embeddings for all FDA class X compounds, including the majority not individually shown, were broadly distributed across both aberrant and canonical regions of the morphospace, reflecting the expected heterogeneity of the historical pregnancy-category system rather than a coherent mechanistic pattern.\nFigure R26 (Revised Supplementary Fig. 7a-b in manuscript). Comparative performance of FATE-MAP and existing in silico teratogenicity models. (a) Predicted morphospace embeddings for the 20 compounds held out from training, with points labeled according to whether the model’s teratogenicity classification matched the known outcome. Boundaries correspond to the failure-mode (pink) and canonical (blue) morphospace regions. (b) Confusion matrices and performance metrics comparing FATE-MAP to two widely used in silico developmental toxicity predictors (CAESAR and P&G) on the same 20-compound benchmark. FATE-MAP matches CAESAR in overall accuracy (75%) and achieves the highest balanced accuracy (85.7%), substantially outperforming the P&G model across all metrics.\nTo strengthen the evaluation and provide a more rigorous benchmark, we replaced the class-based analysis with two complementary validation strategies. First, as described in Response 5, we evaluated FATE-MAP against a curated 33-compound panel with well-characterized teratogenicity outcomes. Second, to assess generalization and compare performance with established in silico developmental toxicity predictors, we conducted an independent analysis using 20 compounds held out from training. Across this set, FATE-MAP achieved 75% accuracy, matching the CAESAR developmental toxicity model and substantially outperforming the Procter & Gamble model (35% accuracy), while also obtaining the highest balanced accuracy (85.7%) (Fig. R26; Revised Supplementary Fig. 7a–b). These results demonstrate that the revised validation framework provides a more reliable and mechanistically grounded assessment than pregnancy-category labels and that FATE-MAP performs competitively with leading computational toxicity models.\n\nComment 16: The authors should also comment on the known phenotypes (at the cellular and organismal levels) associated with Ulixertinib and Nirogacestat in other cellular systems or animals. Do their predicted and observed phenotypes relate to the literature knowledge of these molecules and to their molecular functions (ERK 1/2 and γ-secretase inibitors, respectively)?\nResponse 16: To address the reviewer’s request, we summarize below the known developmental and teratogenic phenotypes associated with ulixertinib and nirogacestat and relate them to their molecular mechanisms of action.\nUlixertinib (ERK1/2 Inhibitor): Ulixertinib targets ERK1/2 kinases, which lie at the core of the MAPK pathway controlling cell proliferation and differentiation.15 Disrupting ERK signaling is known to have severe developmental consequences. For example, genetic ablation of ERK2 (MAPK1) in mice leads to early embryonic lethality, with studies indicating failures in mesoderm formation16 and placental development.17 Consistently, pharmacologic inhibitors upstream in this pathway show pronounced teratogenic effects: the MEK inhibitor trametinib was embryotoxic and abortifacient in animals at exposures as low as ~0.3× the human therapeutic dose (\"Product Information. Mekinist (trametinib).\" GlaxoSmithKline. 2013). These findings align with our observed phenotypes for ulixertinib, as blocking ERK1/2 would be expected to derail essential embryogenic processes, resulting in developmental defects and embryo-fetal loss.\nNirogacestat (γ-Secretase Inhibitor): Nirogacestat is a γ-secretase inhibitor that blocks Notch receptor activation,18 and its known phenotypes in vivo underscore the critical role of Notch signaling in development. While no human pregnancy data are available, animal studies reveal clear teratogenicity: in pregnant rats, nirogacestat treatment during organogenesis caused significant embryo-fetal toxicity, including increased embryo loss, early resorptions, and reduced fetal weights in surviving pups, at maternal exposures below the clinical dose (\"Product Information. Ogsiveo (nirogacestat).\" SpringWorks Therapeutics, Inc. 2023). These adverse developmental outcomes were anticipated given that Notch signaling is indispensable for normal embryogenesis (Notch pathway loss in transgenic mice is embryonically lethal).19,20\nOverall, the phenotypes we predicted and observed for both ulixertinib and nirogacestat are strongly supported by prior knowledge of their molecular targets. The developmental defects identified (in our assays and in the literature) are concordant with the established roles of ERK and Notch pathways in embryonic development, reinforcing the relevance of our findings to these molecules’ known mechanism-of-action and teratogenic risk profiles.\nChanges to the manuscript:\n“Notably, the predicted failure-mode phenotypes for both compounds are consistent with the known essential roles of ERK and Notch signaling in embryogenesis, where disruption of either pathway leads to severe developmental defects and embryo–fetal loss in animal models.”\n\n8) Zebrafish experiments:\nComment 17: This is, according to this reviewer, the most important technical aspect the authors should clarify. Based on both the observed and predicted phenotypes, both Ulixertinib and Nirogacestat treatments result in phenotypic cluster C11, where micropatterns fail to undergo gastrulation and to induce the expression of BRA, consequently lacking the primitive streak (PS). The timing of the treatment in zebrafish embryos is correct, matching the epiblast stage in human micropatterns. The obtained phenotype though does not seem to be entirely compatible with the lack of PS induction observed in C11 micropatterns. If it is true that the exposed animals have a short-tail phenotype, this does not resemble a scenario where BRA (ntl in fish) is not expressed (see PMID: 8402905). In fact, especially for Ulixertinib, axial elongation is still present and somitic rows and other PS derivatives are visible. This might have various sources of explanation, but the authors should clarify why. The authors could also check earlier stages (e.g. late gastrula) to assess the initial formation of early mesoderm (which should lack or be aberrant if the drug induces a failure in PS induction). Overall, the teratogen classification is correct, but the phenotypic prediction is not accurate (or may vary between species).\nResponse 17: Thank you for pointing out this important aspect. We agree that the Ulixertinib (ERK pathway inhibitor) and Nirogacestat (γ-secretase/Notch pathway inhibitor) treatments in zebrafish do not phenocopy a complete absence of primitive streak induction, even though our human stem cell micropattern model predicted a Cluster C11 phenotype (failure of gastrulation/BRA expression).\nUlixertinib: In contrast to a genetic null mutation, pharmacological inhibition by Ulixertinib likely resulted in partial suppression of the pathways needed for PS formation. Consequently, zebrafish embryos treated at the epiblast stage exhibited short-tail phenotypes rather than a complete loss of mesodermal structures. This suggests that some Brachyury/Ntl expression still occurs, enabling residual axial elongation. Indeed, literature on FGF/MAPK pathway interference supports this interpretation: Blocking FGF/ERK signaling in early embryos drastically truncates posterior development but does not entirely abolish it.21 For example, expressing a dominant-negative FGF receptor in zebrafish eliminates most structures posterior to the trunk (severe axis truncation), yet even under these extreme conditions some ntl (Brachyury) expression remains detectable in patchy domains.21 Importantly, general mesoderm formation still initiates via other cues (e.g., Nodal signaling), as evidenced by unaffected expression of other mesodermal markers like Snail1 in FGF-inhibited embryos. In our case, Ulixertinib (an ERK1/2 inhibitor) likely attenuates but does not completely block mesoderm/PS induction, leading to a truncated axis with somites and notochord still present, as observed by the reviewer. This aligns with the notion that partial reduction of BRA function yields a shortened posterior body rather than a total failure.\nNirogacestat: Nirogacestat, a γ-secretase inhibitor, blocks Notch signaling. Notch is known to be crucial for somite segmentation and certain aspects of neurogenesis, but it is not required to initiate the primitive streak in the early embryo. Thus, inhibiting Notch in zebrafish would not be expected to abolish BRA/ntl expression at gastrulation outright; instead, it leads to defects in how the mesoderm is patterned and segmented. Consistent with this, treating zebrafish embryos with γ-secretase inhibitors like DAPT produces phenotypes characteristic of Notch loss-of-function: e.g., irregular or fused somites and a curved, shortened trunk/tail due to segmentation failures, while mesodermal tissues are still present.22 This study showed that γ-secretase inhibitors cause dose-dependent curvature of the tail/trunk without eliminating axial tissues. In our observations, Nirogacestat-treated fish embryos had somites and partial axial structures (short tail), indicating that the primitive streak was induced but subsequent development was aberrant. We will clarify that our “Cluster C11” prediction in the human model (lack of BRA expression) does not perfectly translate to zebrafish’s morphological outcome for a Notch inhibitor; rather, the compound’s teratogenic effect lies in disrupting posterior patterning (leading to truncation) more so than preventing PS initiation.\nOverall, we agree that the teratogen classification of Ulixertinib and Nirogacestat is correct, but the exact phenotypic manifestations differ between the human in vitro model and the zebrafish embryo. Species-specific differences and the degree of pathway inhibition likely explain why we see a short-tail instead of a complete PS failure in zebrafish. We clarify this point in the revised manuscript. In particular, we note that BRA/primitive streak induction is only partially impaired by these treatments in zebrafish, resulting in a truncated body axis (short-tail) with some primitive streak derivatives still forming (e.g. somites and notochord), whereas the human stem cell micropattern prediction suggested a more absolute loss of PS. This variance is consistent with literature reports across species, where complete loss of Bra/T results in no tail, but transient or partial pathway inhibition allows limited axis formation.\nChanges to the manuscript:\n“Although ulixertinib and nirogacestat produced the expected posterior truncation in zebrafish, axial tissues such as somites and notochord remained present, consistent with partial, rather than complete, disruption of primitive streak–derived mesoderm and reflecting species-specific differences in pathway sensitivity.”\n\n9) WNT diffusion modeling\nComment 18: WNT diffusion as a function of tissue density is a very interesting hypothesis and might explain the evolutionary need for a correct PS size during gastrulation. This reviewer’s doubts related to cell death during differentiation (see comment above) remains in this assay as well, and this potentially confounding factor should be ruled out before interpreting the phenotype only on the ability of WNT to diffuse (especially because the lower diffusion of WNT was imposed in the mathematical model).\nResponse 18: We appreciate the reviewer’s thoughtful comment and agree that it is important to distinguish density-dependent WNT signaling from potential confounding effects of differentiation-associated cell death. We took several steps to rule out toxicity-driven artifacts:\nFirst, as detailed in our responses to Comments 11–12 and in Supplementary Fig. 3, we quantified nuclei counts for every colony and defined an explicit toxicity threshold (≤50% of the BMP4-only control mean). Only a small set of strongly cytotoxic compounds fell below this threshold, and these mapped almost exclusively to failure-mode clusters (C1, C11, C12). All analyses of WNT dynamics (Fig. 4) and density–patterning relationships were restricted to the non-cytotoxic regime, ensuring that the density variation examined reflects viable differences in colony size rather than collapse due to cell death.\nSecond, the experimental validation of the WNT model does not depend on the imposed reduction of WNT diffusion in the reaction–diffusion equations. The core result that WNT activity spreads faster and deeper in low-density colonies, and more slowly in high-density colonies, was measured directly in BMP4-only controls using a live β-catenin reporter (Fig. 4e–g). This effect is bidirectional: decreasing initial seeding density accelerates WNT target activation, while increasing initial seeding density attenuates it. A unidirectional confound such as differentiation-associated cell death cannot generate this symmetric relationship.\nThird, this density-dependent modulation of patterning is not restricted to BMP4-only cultures. Across the full canonical patterning region (C2–C10), which contains diverse, non-cytotoxic drug treatments that modulate proliferation and growth rates, BRA expression remains inversely correlated with final cell density. This conserved trend across drugs with different mechanisms further supports that cell density itself, rather than toxicity or death, is the primary driver of the observed variation in WNT/BRA dynamics.\nTogether, the toxicity-filtered cell counts, the live β-catenin reporter measurements, and the consistency of the effect across the canonical morphospace demonstrate that density-dependent WNT signaling is an intrinsic property of the 2D gastruloid system and not an artifact of cell death.\n\nComment 19: Density correlates with WNT waves: this represents the main finding of the manuscript, which also has the potential to inform some basic mechanism during early development. It is indeed an interesting observation, though not completely unexpected. The effect of structure sizes on WNT response has been studied before, also in 3D mouse gastruloids (PMID: 38491138). This work highlighted how varying the size of the initial aggregate does not result in differences in the PS-like domain induction (BRA, CDX, SOX2 patterning). The gastruloid protocol is entirely based on WNT induction using CHIR (small molecule). It could be that the difference in results observed by the two studies relates to the fact that CHIR permeability/diffusion might follow different rules than the WNT itself, being less sensitive to cell density/aggregate size. The authors should compare their results with these and/or try to model WNT vs CHIR diffusion.\nResponse 19: We thank the reviewer for this valuable comment. The reviewer is correct that size-dependent effects on WNT signaling have been examined in 3D mouse gastruloids, including the study cited (PMID: 38491138). Importantly, that work used CHIR99021, a small-molecule GSK-3 inhibitor, as the sole driver of WNT pathway activation. CHIR acts downstream of the receptor by stabilizing β-catenin directly and therefore bypasses WNT ligand production, secretion, and extracellular diffusion entirely. Because CHIR is a small molecule with well-defined permeability and rapid intracellular access, its activity is expected to be insensitive to tissue size, density, or extracellular-space constraints, which aligns with the authors’ conclusion that primitive streak–like patterning scales robustly with gastruloid size.\nThis is directly compatible with our observations. In our 2D gastruloid model, WNT signaling is not induced pharmacologically but instead arises endogenously downstream of BMP4, making it subject to ligand secretion, short-range distribution, receptor availability, and local physical constraints. Under these conditions, we observe a bidirectional density dependence: low-density colonies show faster and deeper WNT propagation, whereas high-density colonies show attenuated and spatially restricted signaling (Fig. 4e–g). Indeed, in our chemical screen, CHIR-treated 2D colonies violate these diffusion-dependent rules: they uniformly express BRA and lose SOX2 across the colony (resulting in colonies exclusively in C1), precisely because CHIR circumvents the need for WNT ligand transport. Thus, the differences between the two studies can be explained by the fundamental mechanistic distinction between a diffusing morphogen (WNT) and a cell-permeant pathway agonist (CHIR).\nMore broadly, both results support a unifying principle: the endogenous WNT system is sensitive to tissue geometry and cell density, whereas synthetic pathway activation masks these dependencies. As highlighted in our Discussion, our findings are consistent with in vivo data showing that cell density is actively regulated during gastrulation. Mouse studies have shown that gastrulation is delayed when early embryos contain fewer cells, whereas aggregated double-sized embryos still initiate gastrulation on schedule despite having more cells at early stages. Moreover, whole-embryo imaging demonstrates minimized density variance at the primitive streak. These observations, together with our 2D results, suggest that endogenous WNT signaling may function as a mechano-sensitive density detector, helping ensure robust initiation of the primitive streak even when starting conditions vary.\n\nMinor comments:\nComment 20: The authors focus extensively on the undifferentiated SOX2+ and PS-like BRA+ domains and less on the trophectoderm-like GATA3+/CDX2+ outer ring. Are there molecules that cause specific aberrations in this domain (e.g. its expansion or complete loss)?\nResponse 20: In our system, GATA3/CDX2 expression is directly driven by BMP4 and is normally confined to the colony edge because BMP receptors are enriched and exposed at the perimeter, while receptors in the center are shielded by cell–cell contacts. Since all compounds in our screen are administered together with BMP4, many drug treatments do not disrupt the outer GATA3 ring simply because the initiating BMP4 cue remains present and the edge-localized receptor geometry remains intact.\nFigure R27. Perturbations that selectively expand or eliminate the trophectoderm-like GATA3 domain. (a) ROCK inhibition with Y-39983 HCl disrupts junctional tension and exposes BMP receptors across the colony, leading to uniform BMP4 signaling and expansion of GATA3 expression throughout the colony (Cluster C12). (b) BMP pathway inhibition with LDN-214117 blocks BMP4-mediated signaling and nearly abolishes GATA3 expression, producing colonies that resemble mTeSR-only controls (Cluster C11). BMP4-only and mTeSR-only conditions are shown as positive and negative controls, respectively.\nHowever, several classes of molecules do generate clear and mechanistically interpretable perturbations of this domain, as shown in Fig. R24. ROCK inhibition (Fig. R27a) provides a prominent example: treatment with the ROCK inhibitor Y-39983 HCl disrupts junctional tension and eliminates receptor shielding, allowing BMP receptors throughout the colony to become uniformly accessible. Under these conditions, BMP4 signaling becomes global rather than edge-restricted, resulting in a nearly uniform expansion of GATA3 expression and producing the asymmetric, TE-expanded phenotype characteristic of Cluster C12.\nConversely, inhibition of BMP signaling using LDN-214117 (Fig. R27b) nearly abolishes GATA3 expression. In this case, the upstream BMP→WNT→Nodal cascade fails to initiate, yielding colonies that resemble the mTeSR-only negative control and corresponding to Cluster C11, where both trophectoderm-like identity and primitive streak induction are lost. Thus, while most compounds preserve the GATA3 ring because BMP4 is supplied exogenously, drugs that modulate receptor accessibility (ROCK inhibitors) or block pathway initiation (BMP inhibitors) produce robust and interpretable alterations of the trophectoderm-like domain consistent with our model’s assumptions.\n\nComment 21: The authors do not check for endoderm induction, a domain observed in the original protocol (PMID: 24973948). While this reviewer understands the challenges in staining for an extra marker, the authors should comment on the expected outcomes or acknowledge this as a possible limitation of the study (apparently canonical phenotypes might mask defects in endoderm induction which were not measured in this study).\nResponse 21: We agree that the absence of an endoderm marker limits our ability to assess defects in definitive endoderm induction within otherwise canonical phenotypes. While prior work has shown that SOX17⁺ endoderm emerges in this system, our assay focused on SOX2, BRA, and GATA3 as primary readouts of ectoderm, mesoderm/PS, and trophectoderm-like identities. We now acknowledge in the Discussion that not measuring endoderm markers represents a limitation of the study, and we note that future extensions of FATE-MAP will incorporate additional markers, including SOX17 and FOXA2, or spatial transcriptomics approaches to more fully capture endoderm specification.\nChanges to the manuscript:\n“We also note that the absence of a definitive endoderm marker (e.g., SOX17 or FOXA2) may mask potential defects in endoderm induction within otherwise canonical phenotypes.”\n\nComment 22: The authors should provide high-resolution imaging of some (or all in supplementary figures or repositories) conditions since the micropattern digital imaging panels are very hard to inspect.\nResponse 22: We thank the reviewer for this comment. In the revised submission, we have uploaded high-resolution versions of all representative images, including both main-text figures and full supplementary panels.\n\nComment 23: The link to access the entire morphospace dataset seems to be broken and I could not inspect that closely (lines 199-200).\nResponse 23: We thank the reviewer for bringing this to our attention. We have confirmed that the link is fully functional and now directs to the complete morphospace dataset, which includes all colonies as high-resolution images: https://max-wilson.mcdb.ucsb.edu/research/gastruloid-morphospace\n\nComment 24: Fig S4 e-h are barely visible and with low image quality, making them quite hard to read.\nResponse 24: In the revised submission, Fig. S4e–h has been replaced with higher-resolution panels and improved contrast to ensure the features are clearly visible.\n\nComment 25: The overall image quality is low and does not allow detailed inspection of the expression patterns and phenotypes. The authors should upload high-resolution figures.\nResponse 25: In the revised submission, we have uploaded high-resolution versions of all main-text figures and full supplementary panels.\n\nComment 26: Despite being published in ref 45, the authors should detail how the membrane bound β-catenin signal was masked in the live cell imaging experiments, since no dyes to label cell membranes were used (at least not specified in the methods nor in the legends). Clarification is needed given that this is a non-trivial segmentation.\nResponse 26: The procedure for masking membrane-bound β-catenin is fully described in the Methods section; however, the main text did not previously direct the reader there. We have now revised the manuscript to explicitly reference the Methods section in this context, ensuring that readers can easily locate the detailed description of how membrane and non-membrane β-catenin were computationally separated.\n\nComment 27: Lines 512-513: the use of the term “malformation” to describe an aberrant phenotypic pattern is not ideal. I would recommend using “aberrant gene expression patterns” or similar expressions.\nResponse 27: We thank the reviewer for this suggestion. We have updated the text accordingly and removed the inadvertent use of the word “malformation” to describe an abhorrent phenotypic pattern.\n\nComment 28: The name of the tool FATE-MAP is a bit misleading, since the current study mainly provides information about the ability to exit from pluripotency (one transition towards nascent mesoderm) and does not go beyond that measurement into fate decision-making.\nResponse 28: We appreciate the reviewer’s perspective. The acronym FATE-MAP was chosen as a compact descriptor of the framework rather than as a literal enumeration of all fate decisions assayed in the present study. We agree that the current implementation primarily captures early fate transitions (particularly exit from pluripotency and progression toward nascent mesoderm) and we now clarify this explicitly in the revised manuscript. While we retain the acronym for consistency with prior use and code infrastructure, we note in the text that future extensions of FATE-MAP will incorporate additional markers and modalities to broaden the range of developmental decisions captured by the platform.\nChanges to the manuscript:\n“For clarity, we note that FATE-MAP is used as an acronym for the framework and does not imply comprehensive mapping of all downstream lineage fates in its current form.”\n\nREFERENCES"}]

Decision letter