Engineering a Computable Epiblast for in silico Modeling of Developmental Toxicity

1. Introduction

New approach methods (NAMs) that generate bioactivity data from automated high throughput screening (HTS) in vitro assays offer the potential to accelerate chemical hazard assessment while reducing the reliance on animal testing [1]. Vast collections of HTS data from ToxCast/Tox21 [https://comptox.epa.gov/dashboard] provide an open resource to examine cellular and molecular determinants of toxicity for drugs and chemicals. Pluripotent stem cell (PSC) culture is one of the most promising in vitro alternatives to pregnant animal testing for assessing developmental hazard potential of drugs and chemicals [2]. Their capacity to form most cell types in the embryo, together with their potential for self-renewal and self-organization into rudimentary tissues, has motivated the design of engineered microsystems that recapitulate various aspects of anatomical development [3,4,5,6,7,8,9,10,11].

Computational agent-based models (ABMs) that mobilize cellular dynamics to recapitulate self-organizing and emergent phenotypes can be used to predict systems-level outcomes linked to biomolecular lesions [12,13,14,15,16,17]. ABMs engineered with sufficient biological veracity allow a user to simulate limitless ‘what-if’ scenarios to localize and analyze critical phenomena and develop hypothesis-based testing strategies [18,19,20,21,22,23]. In these cases, nature-inspired agents (cells) and rules (behaviors) were set into motion using an open-source modeling environment (e.g., compucell3d.org). Information is specified from the known biology to enable bottom-up reconstruction of key morphogenetic events, and to facilitate understanding the global behavior of complex systems following perturbation. This may include effects at multiple biological levels, for example: molecular and physical properties of the system; cellular fate, behavior, and signaling dynamics; and systems-level properties such as emergence, self-regulation, and buffering. Phenotype-based protocols for multiscale modeling and simulation can employ mathematical algorithmic selection or fuzzy logic to fill in for missing biological information. Although these microsystems lack all essential features of in vivo development, they harbor enough information to pinpoint critical phenomenon for probabilistic rendering of where, when, and how a particular insult might lead to an adverse developmental outcome (cybermorph) [16].

In the present study, a cellular ABM was built for mesodermal specification during gastrulation in mouse and human development. This process establishes the anterior-posterior axis within an embryonic disc composed of epiblast and hypoblast layers, where inductive cell signaling centers pattern molecular specification of mesodermal linages [24,25,26,27]. In mammals, epiblast cells internalize at the posterior midline by epithelial-mesenchymal transition (EMT), marking the position of a primitive streak as a transitory hallmark of mesodermal specification [28]. The molecular profile of the epiblast resembles the primed pluripotent state of epiblast stem cells prior to entering primitive streak, and a determined state committed to regional mesodermal specification once they pass through the streak [29]. There is little mixing of the different mesodermal fields once cells exit the primitive streak. This suggests that when and where a cell exits the streak is of critical importance to determining its anatomical fate [30]. As such, cell position in the epiblast is a physical determinant of mesodermal patterning by setting the distance an individual epiblast cell must travel to reach the primitive streak [31].

Nascent mesodermal populations are patterned in part by the nested expression of homeobox (Hox) genes activated in an anterior-posterior sequence over time [32]. A conceptual ‘Hox clock’ runs autonomously to activate this colinear sequence [33]. Some of the key signals regulating this progression include, for example, fibroblast growth factors (FGF4, FGF8), bone morphogenetic protein 4 (BMP4), Wnt Family Member 3 (WNT3), and all-trans retinoic acid (ATRA) [7,34]. To computationally model quantitative impacts of genetic errors or chemical disruption on epiblast development, we built an ABM, referred to hereafter as ESABM (epiblast stem cell agent-based model). The model is driven by an ontology map of the signaling network of a bilaminar embryo and implements a physical progression of epiblast dynamics to simulate mesodermal fields and their disruption by biochemical alterations that can be tied to NAMs data from ToxCast/Tox21.

2. Methods and Implementation

2.1. Model scope

The morphogenetic series covered in the ESABM prototype focused on epiboly of pluripotent epiblast cells to the midline, their ingression by EMT forming nascent mesoderm carrying nested Hox gene addresses, and regional specification of four regional domains: chordamesoderm, paraxial mesoderm, lateral plate mesoderm, and posterior mesoderm (explained in the Results). The simulation begins with an initial tissue configuration corresponding to the bilayered embryonic disc at the egg cylinder stage, just prior to the onset of gastrulation, and progresses through formation of a trilaminar embryo at mid-gastrulation. This approximates a timeline from embryonic (E) days E6.25 – E7.25 in the mouse and gestational days 14-17, Carnegie Stage (CS) CS6 – CS7 in human development. Figure 1 shows the workflow for parameterization, representation, and operationalization.

2.2. Model construction.

ESABM was engineered as a three-dimensional (3D) model in CompuCell3D [https://compucell3d.org/], an open-source environment for simulating the spatial-temporal dynamics of multicellular systems given specified cell behaviors, regulatory signals, and physical properties [12,14,16]. CompuCell3D uses a cellular Potts model (CPM) that treats individual cells as autonomous agents interacting in a shared physical environment, modeled as a discrete lattice, in which each cell occupies multiple lattice sites. The version accessed for this working prototype is CompuCell3D v4.3.2, scripted with Python v3.9.2. The CompuCell3D lattice was 240 x 240 x 48 (2,764,800) voxels. The bilaminar embryonic disc was arbitrarily set to a diameter of 200 voxels filling 54.5% of the lattice from a dorsal view of the epiblast surface (31,400 voxels) and a total volume of 331,650 voxels (12% of the lattice) for the full embryonic disc.

Although most information is based on the mouse, initial rendering was designed to be disc-shaped rather than cup-shaped to better reflect the human embryonic disc [26]. We do not know the exact number of cells in a human embryo at the time of gastrulation; however, a parallel increase in growth expands the epiblast in both mouse [30] and human [35]. Therefore, species equivalence is assumed for the purposes of this model. About 120 cells make up the epiblast of the mouse embryo at E5.5, and about 660 cells at E6.5 [36]. With gastrulation, cell numbers increase in the embryonic region to about 15,000 at E7.5 [37]. Our initial rendering consisted of two epithelial layers, one for the epiblast and the other for hypoblast composed of, respectively, 3483 and 3150 cells, with gastrulation forming a trilaminar embryonic disc as pluripotent epiblast cells ingress through the primitive streak to form the middle mesodermal layer. The target parameter in CompuCell3D was set to maintain individual cell sizes between 33-67 voxels (average 50 voxels), corresponding to biological cell sizes of 9 – 11 µm (average 10 µm). With little extracellular space in the epiblast layer, this range allowed flexibility of cell shape, movement, and growth in the simulation within a discrete size range for epithelial and mesenchymal cells scaled to individual cells in the mouse embryo. An appropriately abstracted geometric boundary condition was fixed to scale the ESABM to an intermediate composition within the timeline for mouse (E6.25 – E7.25) or human (CS6 – CS7) development.

Model dynamics were discretized in Monte-Carlo steps (MCS), each corresponding to 240 × 240 x 48 = 2.8M attempted changes following a Metropolis Monte Carlo algorithm that repeatedly attempts to change which cell-voxels occupy randomly chosen lattice sites. The success of attempted changes is governed by corresponding changes in a pseudoenergy function (ΔE), where favorable changes (ΔE ≤ 0) are always successful; unfavorable changes (ΔE > 0) are successful with probability: P = e^−ΔE/T_eff where T_eff is the cell-specific fluctuation amplitude or effective temperature. This adds a cell-level stochastic element to the model framework (see [21] for details). Each simulation ran for 7500 MCS. Video tiles were imaged each 25 MCS and constructed into movies using Windows video editor. In addition, we queried CompuCell3D for morphodynamic meta data (e.g., cell counts, migration, differentiation) output to flat files for quantitative analysis.

2.3. Cell types and signaling network:

Specifications were developed from the literature, focusing on several key signaling pathways: FGF4, FGF8, BMP4, WNT3, and ATRA. We used Abstract Sifter v7 for entity identification [38] and mined the PubMed literature for the following terms relating to gastrulation in the mouse/human: cell types to be incorporated into the model, key signaling pathways regulating gastrulation, homeobox gene expression during epiblast differentiation, temporality and cellular localization of expression to define parameters for our model. The ontologies, supporting literature, and color schema for the various cell types included in the ESABM prototype are synopsized in Table 1.

Information from literature mining elucidated the flow of molecular regulatory information in a qualitative control network visually represented with BioTapestry v7.1.2 [https://biotapestry.systemsbiology.net/]. To address information gaps in the downstream pathways we used the Mouse Genome Informatics (MGI) resource [https://www.informatics.jax.org/] and the Reactome Pathway Browser for human systems [https://reactome.org/userguide/pathway-browser]. Although the ontology map is gene-centric, the individual relationships represent ligand-receptor interactions and downstream targets for this cell-oriented model. Each ligand is placed in a subregion representing the cell type from which it is secreted. Receptors were input as intercellular nodes within the cell types as they are known to be expressed in the BioTapestry. Downstream targets between the binding of one gene to its cognate receptor and expression of the next downstream gene (such as protein complexes) were input to BioTapestry as ‘bubble nodes’ to illustrate the path by which gene expression evolves.

Mesoderm induction and patterning is largely regulated by NODAL activity through Activin receptors and SMAD2/3 transducers. Other signals patterning the epiblast include NODAL from the epiblast and BMP4 from the posterior visceral endoderm (PVE). The anterior visceral endoderm (AVE) is a major source of LEFTY1 (and by inference, LEFTY2) signals that antagonize NODAL and BMP4 once gastrulation has started. We recognize some key regulators not included in this model implementation will be left for future development; for example, CRIPTO and DRAP1 inhibit NODAL during pre-gastrulation stages in the mouse and have an impact on the posterior mesoderm [55].

(i)

Start FGF4: FGF is required for derivation and maintenance of a primed pluripotent state for human embryonic stem cells (hESCs) in vitro, and for post-implantation growth of the embryonic disc in vivo [25]. FGFR activation preserves hESC pluripotency by transducing information to the nucleus leading to negative regulation of FOXO1 phosphorylation [56,57]. FRS2a (FGFR substrate 2) is an adapter protein linking activated receptor to downstream transducers that regulate the epiblast-stem cell fate decision switch between pluripotency (eg, FOXO1) or differentiation (e.g., RAF/ERK). PI3K-AKT-FOXO1 signaling emerged as the top pathway domain of developmental toxicity from the ToxCast hESC data [58]. It is thus relevant to the epiblast which depends on FGF signals from the endoderm [59]. FGF/FGFR signaling is mediated by the activation of RAS-mitogen-activated protein kinase (MAPK), phosphatidylinositol-4,5-bisphosphate 3-kinase (PI3K)-AKT, Phospholipase C Gamma (PLCγ), and signal transducers and activators of transcription (STAT), which intersect and synergize with other signaling pathways such as WNT, ATRA and transforming growth factor (TGF)-β [60]. Cerebrus (CER1) as an antagonist of NODAL/BMP/WNT is expressed near the AVE in co-expression with LEFTY1/2 [25]. The model represents broad expression of FGF4 and FGF8 expression localized to the primitive streak as a transient function of WNT signaling.

(ii)

Start ATRA: ATRA signaling is potentially active from early PSC stages [61], but the embryo does not make its own ATRA until mid-gastrulation (~E7.5) [62,63]. RAR/RXR liganding suppresses FGF-signaling and triggers Nodal/LEFTY1/CER1 antagonists in the AVE, ultimately influencing mesodermal Hox gene specification through CDX. The pluripotent epiblast is protected from premature differentiation until exposure to inductive cues in strictly controlled spatially and temporally organized patterns guiding fetal formation [31]. This includes precocious ATRA, an endogenous regulator of the patterned Hox gene expression. The removal of ATRA is required for correct Nodal expression during early embryonic patterning and is primarily due to Cytochrome P450 CYP26 enzymatic expression (primarily Cyp26a1) in the extraembryonic tissues surrounding the embryo E5.5 – E6.5 [39,40,54]. Cyp26-deficient mouse embryos show defects that resemble mutants lacking Lefty1 and Cer1 due to up-regulation of Nodal activity [64]. Similarly, exogenous maternal ATRA (oral, 50 mg/kg on E5.5) has similar effects suggesting the embryo is normally protected from premature ATRA exposure but may succumb to defects with retinoid drugs (e.g., tretinoin, isotretinoin) or CP26 inhibitors (e.g., triazole antifungal agents) [40]. TCF7 (gene for TCF-1, the WNT nuclear effector) has overlapping expression with SMADs.

3. Results

3.1. Morphodynamics:

Executing ESABM produced a quasi-normal simulation of the transformation of the embryo from bilaminar (epiblast, hypoblast) to trilaminar disc, representing the sequence of morphogenetic events from E6.25 – E7.25 (mouse) and 3rd week human gestation (CS6-CS7). The normal run encompasses 7500 MCS from start to finish. Epiblast cells are drawn to the midline by cell movement and cell growth (epiboly). This behavior is stimulated by FGF4 released from the hypoblast (visceral endoderm) across the entire simulation (MCS 0 to 7500). Biologically, the posterior visceral endoderm (PVE) releases BMP4 that positions the primitive streak and the anterior visceral endoderm (AVE) releases BMP4 inhibitors to polarize the antero-posterior axis [46,47]. AVE and PVE signaling centers are specified in the model code and thus release their signals that form gradients from the start of the simulation (Figure 2). AVE position in vivo is influenced by mechanical stress from the uterus which is partly responsible for the cup versus disc shape of the bilaminar embryo in mouse versus human, respectively [25]. This is not, however, represented in ESABM where the epiblast borders are circumferentially fixed. Epiboly of the in silico epiblast is therefore limited to inward expansion, represented by a midline primitive streak coded into the model that represents a restricted location of EMT. The timing, shape, and extension of this feature could be specified at different times and positions; however, for the purposes of this small working prototype an ’I-shaped’ primitive streak was specified at 100 MCS to represent the onset of gastrulation at E6.25 (CS6).

3.2. Primitive streak formation:

Specification of the primitive streak was modeled as an open slot for epiblast cells undergoing EMT to migrate through. This reflects a simple abstraction of the biological condition, where EMT in the posterior epiblast is initiated by a localized breakdown of the basal lamina and loss of E-cadherin expression [65]. The inward flux and movement of pluripotent epiblast cells in ESABM was, however, primed by position and timing [66]. Here, we focus on computational analysis of nascent mesodermal populations as the primary emergent output.

3.3. Mesodermal lineages:

Nascent mesoderm cells acquired a nested Hox address as they internalized (Figure 3). The cell field in the simulation (Figure 3A and Supplemental video) shows the bilaminar embryo viewed from above, with the epiblast rendered transparent for quantitative viewing of mesodermal populations. Mammals have 39 Hox genes arranged in 4 paralogous clusters numbered from 1 to 13 within each group. Pluripotent stem cells are usually devoid of Hox expression, but these genes are activated in spatial and temporal patterns defining the development of various body regions [32]. For simplicity, the ESABM prototype simulated only four Hox genes: Hoxd4, Hoxd9, Hoxd11, and Hoxd13 reported in the self-organizing mouse PSC-derived in vitro gastruloid [67].

Epiblast cells adopt different fates based on positional information determined by an autonomous HOX clock encoded into the model. As each epiblast stem cell passes through the streak, its mesodermal identity is locked in time sync to the combinatorial Hox profile [33], modeled in ESABM for colinear expression of the four Hox genes. Each MCS advances the HOX clock by one tick that is counted as different mesodermal lineages form. Each tick of the autonomous clock in ESABM computes to 0.087 sec biological time in the mouse and conversely each minute of biological time in the mouse epiblast is represented in the model by 5.2 MCS.

3.4. Dynamical regulation:

During early development, Cdx2/4 genes follow a similar expression pattern to Hox genes in conferring positional identity to mesodermal tissues derived from posterior epiblast and posterior primitive streak [68]. Biologically, the Hox clock is internally regulated by a small network of signals acting through Cdx genes [45]. The ATRA-WNT-FGF-CDX-HOX signaling axis in ESABM therefore regulates mesodermal fate specification based on levels of CDX activity. Timing in the model is a stochastic feature dependent on a cell’s spatial position in the epiblast that determines how long it took an individual cell to reach the primitive streak as well as where in the primitive streak ingression occurred, forming chordamesoderm at the streak-node region early in the simulation to posterior mesoderm late in the simulation as combinatorial Hox addresses are fixed at the time mesodermal cells release from the epiblast (Figure 4). ESABM showed little mixing of the different mesodermal domains once cells exit the primitive streak so when and where a cell exits the streak is paramount to its anatomical fate [30]. A rudimentary caudal growth zone forms but falls short of reaching a tail bud stage (E7.5) [33]. This process in silico differs from the in vivo condition but recapitulates synchronicity between the acquisition of genetic positional information over time (temporal colinearity) and rostral-caudal axis (spatial colinearity) during the initiation phase of primitive streak formation. For this reason, we next tested the model for predicted effects on disruption of the signaling network.

3.5. Perturbing the control network:

Whereas mesoderm formation in silico is driven by FGF4, FGF8 signaling from the primitive streak stimulates cell migration of nascent mesoderm into the central and peripheral mesodermal wings. Mesodermal cell migration speed is set as a function of FGF8 concentration, which is expressed at the primitive streak to disperse cells following their ingression (Figure 4).

(A)

FGF4 is a positive determinant of CDX-dependent regulation of the HOX clock. Progressive activation of Cdx2/4 specifies progressive mesodermal cell fates (green, chordamesoderm; red, paraxial mesoderm; purple, lateral plate mesoderm; orange, posterior mesoderm), and endoderm (yellow). Trajectories of individual cells in each mesodermal field can be traced in the simulation (see supplemental video S2).

(B)

Pace of the synthetic HOX clock was slowed by FGF down-regulation, thereby affecting colinear Hox activation through Cdx2/4. A 25% reduction in FGF4 had a critical effect on posterior mesodermal specification, with rarefication of cell numbers in this field foreshadowing caudal deficiencies observed in the mouse with 25% to 50% loss of FGFR1 function [53]. The broader effect reflected anteriorization of the mesodermal field shown here at 50% FGF4 reduction with substantial rarefication of the posterior mesodermal population.

(C)

FGF8 signaling stimulates cell migration of nascent mesoderm away from the primitive streak [51]. Mesodermal cell migration speed is set as a function of FGF8 concentration, which is expressed at the primitive streak. A 50% reduction of FGF8 is shown. In the complete absence of FGF8 (not shown), we have no migration away from the primitive streak and the mesodermal cells pile up and slowly expand outwards only due to contact forces from more cells undergoing EMT.

(D)

Retinoid signaling: the embryo acquires its potential for de novo ATRA synthesis around E7.5 [62,63] and is protected from maternal ATRA exposure by the local expression of Cyp26a1 [69]. Precocious ATRA exposure was modeled as an inhibition of CYP26a1 activity at a level that yielded a 50% reduction in Cdx2/4 expression.

(E)

Quantitative mesodermal fields and cell counts computed for 5 replicates (mean + SD for n = 5 runs).

(A)

Regulatory network implemented for abstracted Hox genes in the model, limited for simplicity to sequential activation of Hoxd4, Hoxd9, Hoxd11, and Hoxd13 (3’-anterior, 5’-posterior) as reported in a mouse gastruloid system [67]. An autonomous HOX clock in ESABM is paced by FGF4 and ATRA signaling coordinated by CDX. FGF activation of Cdx2/4 (assisted by WNT, BMP4) speeds the HOX clock for co-linear expression of more posterior (5’) cell fates, according to color addresses in the abstracted model. Responsiveness to ATRA is highest for Hox genes at the 3’ end in a HOX cluster while genes at the 5’ end are more responsive to FGF signaling [54].

(B)

Hox states in the model cell field visualized across the 7500 MCS for a normal simulation. Mid-longitudinal sections of the embryonic disc shown as a snapshot at the designated MCS, with the trilaminar disc oriented anterior (left) to posterior (right). Quantitative mesodermal fields were scored at MCS 7500 for n = 5 replicates (mean + compucell SD). Anteriorization of Hox gene patterning is predicted for FGF4 hypomorphs or ATRA hypermorphs.

Dysregulation of ATRA signaling was modeled in ESABM (Figure 5). Although de novo ATRA synthesis begins circa E7.5 in the mouse, precocious accumulation in the implantation chamber disrupts axial elongation and anteroposterior patterning, especially in the posterior mesoderm [40,70,71]. The embryo is protected from precocious ATRA exposure in the implantation chamber by FGF-induced Cyp26a1 expression in extraembryonic tissues surrounding the epiblast (E5.5 – E6.0) and later in the primitive streak and posterior mesoderm [69]. ESABM does not currently have homeostatic feedback loops for FGF-dependent recovery of Cyp26a1 or Cdx2/4 expression so kinematics of ATRA exposure were modeled as a 50% alteration of Cdx2/4 function. Overall, synthetic perturbation of FGF and ATRA signaling had pronounced consequences on the timing of posterior mesoderm formation consistent with experimental mouse embryology. Reduced FGF4, FGF8 and BMP4 signaling resulted in 86%, 98% and 59% loss, respectively in specification of the posterior mesodermal population. Elevated ATRA exposure resulted in a 78% decrease in this mesodermal population. These effects were highly significant by ANOVA (P < 0.001) across five replicates (n = 5), with the critical effect on FGF4 at or below 25% reduction.

4. Discussion

ESABM provides an in silico platform for the integration of in vitro chemical bioactivity profiles on specific molecular targets with known function in embryology for mesodermal specification. The model is driven by a biomolecular ontology (Figure 2) of the bilaminar embryo and implements a quasi-normal physical progression of epiblast dynamics to simulate nascent mesoderm formation. Although based largely on biological information from the mouse, the initial rendering was physically structured to represent the progression from bilaminar to trilaminar disc stages in the human. ESABM showed little mixing of the different mesodermal layers once cells exit the primitive streak. This is consistent with en toto imaging of the mouse and with the concept that when and where a cell exits the streak is critical for its final anatomical fate [30]. The in silico mesodermal fields produced in the 3D simulation were a quasi-normal facsimile to the pattern from en toto imaging of the early mouse embryo [30], albeit limited by a fixed boundary and absence of the heart field in our model. The in situ epiblast was rendered transparent as epiblast stem cells underwent EMT somewhat akin to the out-of-focal plane of the epiblast layer in the en toto mouse imaging study [30].

Computational dynamics of the ESABM simulation can be used to model countless ‘what-if’ scenarios for predictive biology. Scenarios can be tuned for primitive streak formation, colinear activation of homeobox patterning genes, and cell migration to anterior, paraxial, lateral, and posterior mesoderm fields. The current ESABM prototype arbitrarily implements only four (Hoxd4/d9/d11/d13) of 39 Hox genes in the mammalian genome. They were selected to recapitulate colinear expression during self-organization of the anterior-posterior axis in a murine in vitro gastruloid derived from PSCs [67]. Future expansion of ESABM can incoporate the full complement of Hox genes based on their known expression patterns from the relevant literature and emerging cell atlases of the early embryo [24,25] (see Supplemental Figure S3). This information would enhance resolution of chemical injury predicted on regional Hox codes determined for different mesodermal lineages.

Dual regulation of Hox gene activation by ATRA and FGF signaling is a common theme for shaping the early embryo [63]. ATRA regulates Hox gene activation directly through DR5 elements (primarily anterior Hox gene cluster) [54] or indirectly through Cdx expression in the primitive streak and surrounding epiblast [68,72,73]. The early peri-implantation embryo is protected from premature ATRA exposure (prior to E7.5) by locally expressed CYP26A1 enzymatic activity [52,53]; however, maternal exposure to retinoids during early gastrulation or agents that block CYP26 (e.g, triazole antifungal agents) may be expected to alter Hox patterning at a critical time during early mesodermal fate specification [54,74]. In contrast to ATRA overload, loss of FGF in early mouse embryos is embryolethal at pre-primitive streak stages (E5) for FGF4-deficiency, and early gastrulation (E7) for FGF8 deficiency [75]. Although it may be difficult to define a threshold of ‘lethal’ outcome in the model, the reduction of FGF8 signaling resulted in the progressive accumulation of cells at the primitive streak failing to complete EMT and dispersing from the primitive streak. The current ESABM prototype does not engage apoptosis, which remains for future model development; however, dispersal of nascent mesoderm from the primitive streak is an FGF8-driven cell behavior [75]. This pathway involves down-regulation of E-cadherin in the mouse epiblast [76] and RAS/MAPK and PI3K/AKT pathways in the chick epiblast [51]. FGF8-dependent dispersion of nacsent mesoderm is coarsely represented in ESABM at the level of FGFR1 activity, which is influenced by FGF4 from the hypoblast before primitive streak formation and then by FGF8 in the primitive streak.

Finally, ESABM predicted strong disruption in the timing of posterior mesodermal formation with FGF-deficiency or ATRA-excess, consistent with empirical observations for Cdx-dependent growth of the posterior mesoderm [53,73]. The small working prototype is proof-of-concept for predictive toxicology, particularly for how the embryo forms blood which is substantially derived from the posterior mesoderm [24,70,71]. In this way, the in vitro chemical effects data from largely human cell-based assays in ToxCast/Tox21 [https://comptox.epa.gov/dashboard/] can be imported to a dynamical, self-organizing in silico toxicodynamics platform for dose-time exposure scenarios. As a use case, we queried USEPA’s CompTox Chemicals Dashboard for ToxCast hits with an activity concentration of 50% (AC50) < 1 µM on the in vitro bioactivity profiles for inhibition of FGFR1 biochemical activity in a cell-free assay (NVS_ENZ_hFGFR1) and/or reporter gene activation of the Direct Repeat 5 retinoic acid response element in a cis-activation assay (ATG_DR5_CIS_up). For example, PharmaGSID_48519 (CASRN 686756-87-6) had an AC50 (concentration that altered bioactivity by 50%) of 0.005 µM on hFGFR1 liganding. At 0.36 µM, this kinase directed drug, identified as a SMAD binding element inhibitor at 0.57 µM [77], invoked a gene signature of EMT in the MCF7 cell line. At AC50 = 1.22 µM, GSID_48519 caused cis-activation of DR5 inferring an effect on ATRA-induced gene expression. This concentration was the dTP index in the ToxCast STM_hESC assay inferring a threshold for human developmental toxicity potential at 1.22 µM [58]. In the absence of a direct validation, these lines of evidence from different types of assays for GSID_48519 support each other and add confidence to the simulated outcomes. Chemicals may disrupt FGF signaling in multiple ways including FGF kinematics or effects on FGF-modulated genes, including Hox activation. Chimeric embryoids constructed from different proportions of normal and FGF-deficient cells showed a biologically relevant reduction of FGF signaling around 25% [53], consistent with the cybermorph predicted by ESABM. Computationally different degrees of chemical insult to the kinematics of FGF and/or ATRA systems can thus be modeled with ESABM for predictive toxicology. Although much of ESABM is based on biological information from the mouse, translation of developmental toxicity potential to a human epiblast model can be evaluated with in vitro ToxCast data from embryonic stem cell platforms for human and mouse [58].

5. Conclusions

Virtual in silico tissue models are a novel approach to visualize cellular trajectories, map toxicodynamics, and predict adverse phenotype (cybermorphs) in ways difficult to accomplish in vivo. The temporality of ESABM at present covers some of the morphogenetic series of events from pre-gastrulation (E6.25) through mid-gastrulation (E7.25). We consider the present state of ESABM as quasi-normal because the primitive streak does not elongate, and the Hox-pattern is over-simplified. Proof of concept shows that the model enables integration of in vitro chemical bioactivity data for specific molecular targets with known embryology for mechanistic prediction of altered mesoderm development, shown in a case example with FGF and ATRA pathways. An important limitation with the current ESABM prototype is that it does not simulate anterior extension of the neuraxis, heart field nor caudal growth zones. Those subsequent morphogenetic stages will be left for future model development. However, it’s capacity to translate NAMs data into potential mesodermal phenotypes in silico may be validated by in vitro studies using engineered microsystems such as gastruloids that recapitulate various morphogenetic processes.