High-Dimensional Design-Of-Experiments Extracts Small-Molecule-Only Induction Conditions for Dorsal Pancreatic Endoderm from Pluripotency.

The derivation of endoderm and descendant organs, such as pancreas, liver, and intestine, impacts disease modeling and regenerative medicine. Use of TGF-β signaling agonism is a common method for induction of definitive endoderm from pluripotency. By using a data-driven, High-Dimensional Design of Experiments (HD-DoE)-based methodology to address multifactorial problems in directed differentiation, we found instead that optimal conditions demanded BMP antagonism and retinoid input leading to induction of dorsal foregut endoderm (DFE). We demonstrate that pancreatic identity can be rapidly, and robustly, induced from DFE and that such cells are of dorsal pancreatic identity. The DFE population was highly competent to differentiate into both stomach organoids and pancreatic tissue types and able to generate fetal-type β cells through two subsequent differentiation steps using only small molecules. This alternative, rapid, and low-cost basis for generating pancreatic insulin-producing cells may have impact for the development of cell-based therapies for diabetes.

Introduction

Endoderm is the germ layer that creates the majority of cells within most of the internal organ systems, such as lung, stomach, pancreas, liver, and gut. The ability to robustly generate endodermal descendant tissues will impact the studies and therapy modalities of multiple human diseases. Until present, almost all efforts on inducing endoderm from pluripotent cells have relied on using a TGF-β pathway agonist, most commonly Activin A (AA), to push pluripotent cells through an in vitro gastrulation event (D'amour et al., 2005; Gadue et al., 2006). This results in an endodermal population that can be successfully used for generating multiple descendant fates including intestinal (Spence et al., 2011), pancreatic (Kroon et al., 2008; Rezania et al., 2014; Pagliuca et al., 2014), and liver (Sampaziotis et al., 2015). Generation of more anterior endodermal fates, such as lung, has been achieved by providing patterning inputs at a subsequent stage (Green et al., 2011). However, recent studies argue that initial patterning of definitive endoderm may occur during its generation (Matsuno et al., 2016; Loh et al., 2014).

The pancreas is of particular interest for cell-based therapy in diabetes, which is characterized by defects in, or loss of, insulin-producing cells. The pancreas is formed from two spatially distinct primordia arising on the dorsal and ventral sides of the primitive gut tube, which subsequently fuse. Although both pancreatic buds are capable of generating all lineages of the adult pancreas (Matsuura et al., 2009), distinct transcriptional programs control the initial induction of the pancreatic domains on opposing sides of the gut tube. In mice, the ventral pancreatic bud forms first at approximately embryonic day 8.5 (E8.5) from a region of endoderm possessing bipotential competence for pancreas and liver (Angelo et al., 2012; Deutsch et al., 2001; Tremblay and Zaret, 2005). This early ventral endoderm field consists of a progenitor population that co-expresses Pdx1/Sox17 transiently, which by E9.5 splits to form the ventral pancreas and the extra-hepatobiliary system, respectively (Spence et al., 2009). Specification of the ventral pancreas relies on HHex expression. Gene ablation models have demonstrated complete ventral agenesis without affecting dorsal pancreatic bud formation (Bort et al., 2004). In contrast, the dorsal pancreatic bud in mice emerges at approximately embryonic day 9.0 and forms from an outgrowth caudal to the antral stomach region. Studies in mice have also identified factors involved in dorsal pancreatic specification with no effect on ventral organogenesis. Mnx1 (Hlxb9) knockout models have shown dorsal agenesis occurs without a ventral phenotype (Li et al., 1999). Mnx1 expression is observed in the ventral field but only following Pdx1 expression, whereas in the dorsal field, Mnx1 precedes Pdx1 expression. Raldh2 knockout models resulted in a dorsal-specific agenesis attributed to the loss of Pdx1 and Prox1 expression in the dorsal bud (Martin et al., 2005; Molotkov et al., 2005). Furthermore, studies in chick have shown that the initial budding of the dorsal pancreas is dependent on the selective inhibition of SHH within the dorsal midgut (Hebrok et al., 1998). Although it is unclear if the murine system is conserved between species, a recent study using laser capture followed by deep sequencing analysis described some fundamental differences between the ventral and dorsal pancreas during human development (Jennings et al., 2017).

Despite differential pathway utilization and distinct cell intrinsic factors the dorsal and ventral pancreatic programs have much in common. HNF1β (Tcf2) is required for pancreas specification in both pancreatic buds and is critical through pancreatic development. Tcf2 knockout mice fail to generate a ventral pancreas and have a greatly reduced dorsal bud incapable of differentiating or proliferating (Haumaitre et al., 2005). HNF1β is expressed in the pre-pancreatic foregut and functions at the apex of a sequential transcriptional cascade resulting in the activation of Hnf6 (Oc1) followed by Pdx1 (Poll et al., 2006). Conditional inactivation of HNF1β results in a loss of Glis3 and Ngn3 expression and results in a pancreas characterized with cystic ducts and a loss of the pro-endocrine field (de Vas et al., 2015). In human development, the importance of HNF1β is highlighted by the occurrence of “maturity-onset diabetes of the young type 5” (MODY5) syndrome, a condition attributed to mutations in the HNF1β gene. Although a heterozygous mutation in HNF1β does not display a phenotype in mouse studies, in humans heterozygous mutation of HNF1β have been shown to be associated with MODY5 or complete pancreatic agenesis suggesting a more important role for HNF1β in human pancreatic development than in mouse (Body-Bechou et al., 2014).

Most current laboratory efforts at directing the differentiation of pluripotent stem cells rely on emulating developmental signaling event(s) leading to the generation of desired cell type. This is generally accomplished by assaying a single pro-differentiation factor at time. The limitation of this methodology is that it relies on studying each pro-differentiation factor individually, hindering the detection of synergistic or systemic influences. A way to address this problem is using a systems biology approach (Kitano, 2002; Carinhas et al., 2012) capable of assaying multiple factors simultaneously in a manner capable of elucidating individual and synergistic effects. Because experimental design size increases exponentially as additional factors are incorporated, this greatly limits traditional methods from approaching a systems biology level of interrogation. To overcome this limitation, we have developed a novel approach focusing on key aspects of a manufacturing process (key concepts are defined within Box 1). Using Design-of-Experiments (DoE) mathematics (Chakrabarty et al., 2013), we are able to greatly increase the dimensionality of our differentiation experiments by relying on a compression of the design space (Gerin et al., 2014; Mendes et al., 2016). This systematic approach minimizes the number of experimental runs needed to interrogate multiple parameters simultaneously within a single experimental design (Rathore et al., 2014). Combined with a deep set of lineage-informative transcript level measurements a better understanding of the cell culture behavior is obtained within the design space (Mercier et al., 2013). Such an approach is integral to a Quality-by-Design (QbD) process (Juran, 1993; McConnell et al., 2010; Swain et al., 2018). It provides process understanding, allowing for consistent product manufacturing (Kumar et al., 2014; Lipsitz et al., 2016). In the present study, we use DoE-based optimization to produce a pancreatic directed differentiation protocol defining the critical process parameters relevant to the differentiation process and future manufacturing. We monitor gene expression throughout the differentiation process as a critical material attribute that defines the cellular phenotype. We then identify which pathway control elements are the critical process parameters that must be controlled to ensure proper differentiation.

Box 1

Definition of Key Terms Used within This Study.

Contrary to most methods of endoderm induction, we demonstrate that effective and regionalized patterned endoderm can be robustly differentiated directly from pluripotency without the use of TGF-β agonism. Exploring for optimal endodermal fate conversion conditions, we used HNF1β expression as an initial waypoint for pancreas and other endodermal derivatives. This led to a novel and highly robust protocol for inducing specialized human endoderm representative of the dorsal foregut region of the gut tube. We demonstrate that this population can be effectively converted into dorsal pancreatic progenitors that subsequently are able to adopt endocrine fates, including the generation of fetal-like beta cells. By inspecting the critical process parameters we created a three-stage protocol that converts PSCs into fetal-like beta cells using a series of small molecules.

Results

Application of Systems Biology to Understand Developmental Models

We combined a high-dimensional application of Design of Experiments (HD-DoE) methodology with deep response set measurements to generate predictive models for pluripotent culture forward differentiation. This method allowed interrogation of pathway interactions, providing a more comprehensive understanding of the biological inputs impacting endodermal differentiation. By challenging pluripotent cultures with a perturbation matrix composed of multiple morphogen pathway agonists and antagonists simultaneously, and measuring multiple fate determining genes, we extracted a systems-level effector/response dynamic. The application of DoE substantially compressed the number of experimental runs compared with a full factorial design. Yet it remained possible to statistically determine any first-order pathway/pathway interactions. It also revealed the system behavior covered by the tested dimensions within the concentration ranges used (i.e., “known space”).

The morphogens (effectors) used in initial experimental designs were AA, BMP4, FGF2, WNT3a, SHH, and RA. We also included respective small-molecule antagonists for each pathway except RA. These factors were chosen because previous publications within the stem-cell field have shown these inputs to elicit forward differentiation from pluripotency into various descendant fates. We generated D-optimal DoE designs to test all of these effectors in a single experiment. Designs specified the combinations of effectors to be included in each well of a 96-well plate. The designs were constrained to prevent unproductive combinations of agonist/antagonist of the same pathway being used together in the individual reactions. Robotic assembly of experimental media (i.e., the Perturbation Matrix) eliminated human errors while ensuring accuracy. Response measurements were custom-chosen early lineage-determining genes. Following experiment execution, a multivariable regression model was generated for each response gene as a function of effector contribution. All response data were mathematically fitted to maximize predictive power (Q2 maximization). As a result, we obtained an in silico representation of the behavior of each response gene as related to each effector input. These models allowed us to predict conditions that would achieve desirable induction, or suppression, of any of the genes monitored. We refer to these interrogations as in silico predictive analysis (ISPA), noting that such predictions were calculated on the basis of statistical models resting on the entire set of the DoE design.

For ISPA, various tools were needed to achieve specific desirable outcomes and identify critical process parameters. Coefficient Plots may be generated for each individual response gene. These plots display the coefficient for each effector term in the regression model for the respective gene. coefficient plots are scaled and centered and thus also provide graphical representations of model term significance. We used coefficient plots to inspect individual effectors' contributions to the activation of individual genes. The “optimizer” function used the regression models to identify which media compositions contributed to a desired differentiation event. The optimizer provides the relative Factor Contribution (FC) for each effector. FC is proportional to how important the individual effector is to the differentiation event. FC thereby helps identify criticality of a process parameter. For practical proposes, we considered FC < 10 as low relevance, 10–20 as relevant, and >20 as highly relevant. Particularly valuable, we used the “Dynamic Profiling” to visualize expression behavior for multiple genes simultaneously at any given input condition (set point).

Two Separate and Distinct Pathways Exist for Endoderm Induction

We initially set out to predict conditions needed to define a definitive endoderm population as suggested by the literature. This was accomplished by defining an anterior primitive streak (APS) population by modeling for the maximal expression of MESP1, EOMES, and BRACHYURY/T while minimizing EVX1 (posterior primitive streak marker) (Loh et al., 2014). Through ISPA, the conditions predicted to generate this differentiation event consisted of low tolerance to Wnt inhibition (FC = 18.65) and high levels of AA (FC = 18.85) (Figures 1A and 1B). These conditions agree with current protocols for generating definitive endoderm (DE) (D'amour et al., 2005). Indeed, when using this condition to differentiate pluripotent cultures, a FOXA2+/SOX17 + population was obtained within a 3-day period (Figure 1C). Through ISPA, a number of other genes were predicted to be highly expressed under these conditions including COL6A1, HHEX, MESP2, SOX17 (Figure 1D and data not shown). These genes are all known to be elevated in definitive endoderm. However, through ISPA, we noted that not all known early endoderm-expressed genes were uniquely maximized through the APS conditions. Inspecting the known space from the aforementioned experiment, a quite different solution set could be obtained focusing on HNF1β (TCF2) expression, also known to be expressed in definitive endoderm. Maximizing for HNF1β expression also led to expression of accompanying genes such as FOXA2, HNF4A, MNX1, CXCR4, MTF1, all known markers of endoderm, whereas expression of HHEX and SOX17 remained low (ISPA results not shown). We went on to characterize these distinct states and the requirements for their induction.

Figure 1

Effective Endoderm Induction in Absence of AA/WNT

(A) Schematic of endodermal generation through AA induction of Anterior Primitive Streak (APS) versus endoderm induction optimizing for HNF1β induction.

(B) Predicted conditions satisfying APS gene induction (red) versus HNF1β induction (green).

(C) Validation of endodermal marker induction using the two separate protocols.

(D) Coefficient plots (primary effectors only) for select endodermal genes. Error bars within the Coefficients Plots represent 95% confidence intervals. Significant terms are identified as terms with confidence intervals that do not overlap the y axis.

Retinoic Acid and Bmp Inhibition Synergistically Induce an Endodermal Program Mutually Exclusive to Activin-Induced Endoderm

Using ISPA, we inspected the fundamental logics governing endodermal gene induction. For the APS-derived DE, it was clear that many early endoderm genes were under the direct control of TGF-β signaling, displaying strong and positive coefficient terms from AA in their complex regulatory models. These genes included, but were not limited to, SOX17, CXCR4, LEFTY1, MIXL1, and HHEX (Figure 1D). This was not the case for multiple other known endoderm markers. A sub-group of endoderm genes did not respond to AA stimulation but were directly dependent on RA signaling and to a lesser extent required the inhibition of the BMP pathway. These genes included FOXA2, EPCAM, ONECUT1, CDX2, and MNX1 (Figure 1D) and were predicted to be directly controlled through the synergistic effects of retinoic acid and BMP inhibition with high factor contributions including FOXA2 (FC for RA = 24.8, BMPi = 25.3), HNF1β (FC for RA = 30.1, BMPi = 31.7), and MNX1 (FC for RA = 22.7, BMPi = 30.5) (Figure 1B and data not shown). Of note, AA was predicted to have no contribution to activating these genes; rather, inhibition of the TGF-β pathway was predicted to benefit the expression of these genes. Factor contributions for Alk5i were 7.59 for HNF1β, 8.88 for FOXA2, and 18.98 for MNX1 (Figure 1B and data not shown). To assess the robustness of the methodology, the same perturbation matrix design was used on a second pluripotent cell line, H9 (female) (Figure S1). Both the BMPi (Factor Contribution of 11.0) and the provision of RA (Factor Contribution 29.9) were again shown to be critical factors for the activation of HNF1β. Although the factor contribution for BMPi decreased, there was a corresponding increased importance in the absence of AA (Factor Contribution 20.1 that AA is not included) (Figures S1A’–S1B′). Altogether, this suggests that retinoic acid input when provided in the absence of TGF-β signaling is a critical process parameter for induction of HNF1β. Comparative analysis using dynamic profiling at the HNF1β set point revealed that effectors controlling HNF1β, MNX1 (both retinoic acid and TGF-β inhibition responsive), and HHEX (AA responsive) were very similar between these cell lines, which differ in sex and prior culturing conditions (the male H1 in Essential 8 and the female H9 in mTesR media) (Figures S1C and 1C′). The predicted conditions for HNF1β optimization (HNF1βOpt) were tested on differentiating pluripotent cells for validation. Based on protein expression, HNF1β, MNX1, and FOXA2 could all be activated as expected. The target gene for optimization, HNF1β, was present in 97.9 ± 0.7% of the cells (Figures S2A and 2C).

The two paths to endoderm activation were fundamentally distinct and rested on conflicting input logic suggesting that the pathways were mutually exclusive. We tested this by creating hybrid protocols assaying the effects of RA and BMPi in the presence of AA (Figure 2A). Inclusion of RA into the APS-based DE-generating protocols (D'amour et al., 2005) proved to only moderately increase gene expression for retinoic acid-responsive genes; only CDX2 (data not shown) and OSR1 (Figure 2B) were significantly up-regulated in this manner. Conversely, known TGF-β-responsive genes were shown to be significantly down-regulated when retinoic acid was included in APS-type DE-generating reactions including HHEX, SOX17, and GSC (Figure 2B and data not shown). This demonstrates that presence of either of the key protocol drivers (AA versus RA) will suppress the other. Furthermore, also as predicted by ISPA, the genes up-regulated in the presence of retinoic acid and LDN3189 were activated more efficiently when AA was excluded from these reactions. Importantly, the two key protocol inputs for the HNF1βOPT conditions, RA and LDN3189, sufficed to initiate differentiation comparable with the full-input HNF1βOpt conditions (Figure 2B) with 97.1 ± 1.8% of the cells within the culture expressing HNF1β (Figures S2A and S2C).

Figure 2

Non-APS-Derived Endoderm Is Critically Activated by Retinoic Acid and BMP Inhibition and Demonstrates a Dorsal Foregut Character

(A) Schematic of the experiments performed.

(B) Graphs showing the relative expression of several endodermal genes in response to retinoic acid and BMP inhibition. Experiments consist of quadruplicate biological replicates performed in parallel experiments within a single TC plate. Genes were normalized to the average expression of the endogenous levels of YWHAZ, GAPDH, and TBP. DE, definitive endoderm; RA, retinoic acid; LDN, LDN3189, a BMP inhibitor; WIN, Win 18446, an ALDH inhibitor. Error bars represent standard deviation/triplicate assays.

(C) Heatmap containing key endodermal genes from pluripotent cultures either subjected to classic definitive endoderm differentiation conditions or differentiated using conditions predicted in the HNF1β optimizer.

Retinoic Acid/TGF-β Inhibition-Induced Endoderm Is of a Dorsal Foregut Character

To gain a better understanding of the differing nature of the endodermal populations, we subjected cultures for RNA sequencing. Common endodermal genes were expressed at similar levels in both populations; these included CXCR4, FOXA2, EPCAM, GATA4, and GATA6 (Figure 2C). However, significant differences were observed for genes associated with patterning revealing that HNF1βOpt induced endoderm was enriched in genes characteristic of known dorsal (MNX1 and PAX6) and foregut endoderm (HOXA1, HOXA3, HNF4A, and HNF1β) (Figure 2C). In contrast, APS-type DE showed an enrichment for genes representative of ventral endoderm (NR5A2, HHEX, and SOX17) and more posterior endoderm (SOX17 and AFP), although not CDX1 (Figure 2C). Since the HNF1βOpt culture appeared to have a stronger dorsal foregut endoderm (DFE) phenotype, we challenged it for differentiation competence toward stomach, pancreas, and liver using conditions previously shown to induce these fates from APS-derived DE (Figure S3A). These tissues are derived from posterior foregut, and liver is a ventral derivative only. We also differentiated APS-derived DE for comparison. Differences in competence were observed between the DFE and APS-DE populations. In all cases, APS-DE cultures activated liver genes to higher levels (APOB, HHEX, and EVX1 (Figure S3B) than DFE, whereas stomach (OSR1) and pancreas (PDX1) genes were activated at higher levels in the DFE cultures. We also challenged the DFE culture to generate stomach organoids (Figures S3C and S3D [Mccracken et al., 2014]). OSR1 and PDX1 co-expression and SOX2 and PDX1 co-expression were both observed, suggesting that the stomach organoids preferentially converted into antral-type, posterior-most stomach (Figure S3E).

Organ-Field Specification Mechanisms from Dorsal Foregut Endoderm

Since DFE patterned cultures displayed a competence for PDX1 activation (Figure S3A), we next evaluated DFE pancreatic potential through a sequential DoE modeling experiment (as outlined Table S2). DoE designs included effectors previously shown to induce pancreas as well as effectors known to pattern along the mid-section of the developing gut tube. These included FGF2, FGF4, PD0325901, SHH, Sant1, BMP4, LDN3189, EGF, AA, A8301, and Wnt3a with agonist/antagonist constraints. Results from ISPA maximization of the expression of SOX2, OSR1, or PDX1 were compared (Figures 3A and 3B). Dynamic profiling analysis (Figure 3C) revealed that both OSR1 (gastric) and PDX1 (pancreatic and wider) were under similar mechanisms of control. Both genes were strongly responsive to RA with FC = 30.68 and FC = 31.38, respectively (Figures 3B and 3C). However, differential responsiveness to FGF signaling were predicted. FGF4 was predicted to be important for OSR1 (gastric) activation with an FC of 13.45, whereas MEK pathway inhibition with the small-molecule inhibitor PD0325901 was shown to strongly contribute to PDX1 activation with a substantial FC of 29.62 (Figures 3B and 3C). Thus, the bipotentiality of gastric/pancreatic fates is resolved by FGF/FGFi inputs, respectively. Of note, SOX2 expression was strongly decreased by RA and highly increased by SHH with respective FC of −28.59 and 29.58 (Figures 3B and 3C). Our data argue that retinoic acid secures a posterior antral field and active inhibition of SHH contributes to the down-regulation of SOX2 expression, hereby allowing for a switch from a gastric to a pancreatic field.

Figure 3

DFE Can Generate Pancreatic Endoderm with a Highly Dorsalized Nature

(A) Schematic showing modeling of the optimization for the stomach genes SOX2, OSR1 and the pancreatic gene PDX1. (B) The corresponding optimizers for the predicted maximal induction of SOX2, OSR1, and PDX1, respectively.

(C) Dynamic profiles for the effectors most responsible for the respective gene activation.

(D) Representative IHC of PDX1Opt.

(E) KeyGenes prediction for the respective DE- and DFE-derived pancreatic endoderm.

(F) Heatmap assessing the differential expression of several pancreatic and dorsal-specific pancreatic genes between the two protocols.

DFE-Derived Pancreas Is of a Dorsal Identity

Using the ISPA-defined PDX1Opt conditions (shown in Figure 3B), we next demonstrated that the PDX1-expressing DFE-derived cultures co-expressed several known pancreatic progenitor markers including FOXA2, NR5A2, GATA4, and SOX9 indicating that a true pancreatic endodermal (PE) state was rapidly induced (Figure 3D). This DFE-derived pancreatic induction was shown to be reproducible in the H9 female embryonic stem cell line, as well as in iPSC culture (Figure S4). An RNA-seq-based KeyGenes analysis (Roost et al., 2015) was used to verify that this DFE-derived PDX1-induction was truly a pancreatic fate (Figure 3E). We compared DFE-derived PE with previously published (Xie et al., 2013) APS-type DE-derived stage 4 PE (Figure 3F). Interestingly, both of the starting populations, DE and DFE cultures, initially displayed a similarity to brain, which was lost through the sequential stage in both protocols (Figure 3E). This loss of similarity to brain is likely attributable to emergent expression of neuronal fate suppressor genes (Jennings et al., 2017) activated during the sequential stage(s) in both populations (Figure 3F). KeyGenes analysis showed that the published APS-derived DE population displayed similarity to lung, whereas the DFE population did not (Figure 3E), presumably because the lung buds are exclusively derived from ventral endoderm. A hierarchical clustering comparison between DFE, DE, and their derived PE populations demonstrated that the DFE and DE populations were so similar that the DFE population clustered in-between the triplicate set of DE samples (Figure S5), attesting that both protocols attain a fundamental endodermal program. However, when comparing PE derived from either DFE or DE with genes highly enriched in the dorsal pancreas during human development (Jennings et al., 2017) we found that the DFE population already expressed several of these genes including DLL1, CNR1 FRZB, HOXA1, and ARMC3, whereas the APS-DE population displayed only low expression of CNR1 and FRZB. Also, expression of MNX1 (Hlxb9), a previously described dorsal marker, was expressed throughout the DFE culture, whereas detectable MNX1 expression within the DE-derived PE only began at the PFG stage (Figure 3F). Subsequently, DFE-derived PE continued to express the vast majority of the dorsal-specific genes (Figure 3F).

Figure 4

Biological Equivalence Testing of Small-Molecule-Based Pancreatic Endocrine Induction Protocols

(A) Schematic of the reactions performed in which the HNF1βOpt and PDX1Opt were replaced with only the CPP factors as identified as the effectors with the highest predicted factor contributions, retinoic acid and LDN3189 (R/L) or retinoic acid and PD0325901 (R/P), respectively.

(B) Transcript analysis of the four respective reaction conditions at either the PE stage (indicated in yellow) or the endocrine stage (indicated in red) and normalized to the averaged expression of GAPDH, TBP, and YWHAZ. Experiments consist of quadruplicate biological replicates performed in parallel experiments within a single TC plate.

Time Dependency of DFE-Derived Fetal-like Endocrine Cells

The importance of NOTCH pathway inhibition for terminal differentiation of endocrine cells is well known (Apelqvist et al., 1999; Jensen et al., 2000; Afelik et al., 2012). We examined the temporal effects of NOTCH inhibition toward endocrine commitment using three media inputs. All included a NOTCH pathway inhibitor (γ-secretase inhibitor XX), and we also evaluated single-SMAD (Alk5 inhibitor) and dual-SMAD inhibition (Alk5i/LDN) (as outlined in Figure S6A). Strikingly, prolonging the period of PE induction diminished endocrine competence (Figure S6B). Of note, increasing the duration at the PE stage led to an increase in the competence toward acinar and, finally, ductal differentiation (Figure S6C) as assayed through the expression of CPB1, and MIST1 (BHLHA15) and F3, HNF1β and PROM1, respectively. Inspecting the resulting endocrine cultures, we found PDX1+/NKX6.1+ co-expression throughout the culture with patches of NKX2.2/insulin C-peptide co-expression throughout the culture (Figure S7B). Furthermore, mono-hormonal (INS, GCG, and SST) as well as polyhormonal expressions of INS+/GCG + or INS+/SST+ (Figure S7B) were observed. Maximum C-peptide levels occurred after 10 days of exposure to Notch/Alk5 inhibition (Figure S7C), and cellular aggregates demonstrated dithizone retention (Figure S7D). Classical GSIS and microfluidic analysis assays (Adewola et al., 2010; Wang et al., 2012) were consistent and demonstrated an immature βcell physiology capable of synthesizing and storing insulin but with limited functional profile when responding to glucose fluctuations (Figures S7E and S7F). This DFE-derived endocrine cells also display a fetal-like βcell state as previously observed from APS-DE-derived β cells (Hrvatin et al., 2014).

Development of a Small-Molecule Method for Pancreatic Induction

Inspection of the PDX1Opt conditions identified that retinoic acid and PD0325901 had the highest factor contributions of 31.38 and 29.62, respectively (Figure 4A), whereas all other effectors tested had low factor contributions <10. Seeking to achieve a minimally complex method for pancreatic induction we combined the identified critical process parameters for HNF1βOpt (RA and LDN [R/L]) with the critical process parameters of the PDX1Opt (RA and PD0325901 [R/P]). PDX1 levels at the PE stage varied little between these different combinations of the induction methods. Inducing DFE using the full HNF1βOpt and then switching between the PDX1Opt and just the defined CPP components, RA and PD0325901, resulted in cultures that were either 85.6 ± 2.5% or 85.7 ± 2.5% positive for PDX1 expression, respectively (Figures S2B and S2C). However, when PDX1 induction was assayed on DFE induced using only RA and LDN3189, PDX1 levels varied between 81.6 ± 3.2% (for RA/LDN – PDX1Opt) or 86.7 ± 3.6% (for RA/LDN – RA/PD0325901) positive (Figures S2B and S2C). Endocrine competence was then confirmed with an endocrine fate conversion using the aforementioned combination of NOTCH and ALK5 inhibition. DFE induction using HNF1βopt followed by R/P significantly increased levels of the pancreatic markers assayed, and this increase was amplified when followed by an endocrine induction (Figure 4B). A further increase in pancreatic genes was observed when the HNF1βOpt was replaced with the identified CPP (R/L). Using the R/L step, followed by PDX1opt conditions, generated the highest levels of the pancreatic genes INS, GCG, NKX6.1, CHGA, GLP1R, and PDX1 (Figure 4B). Yet, comparative changes were observed when both the HNF1βOpt and PDX1Opt were replaced with R/L and R/P in conjunction. Slight decreases in INS, GCG, NKX6.1, CHGA, GLP1R, and PDX1 expression were observed under these conditions. However, moderate compensatory increases were observed in the expression of NEUROD1, NGN3, PAX4, MAFA, and SST. We conclude that the only critical process parameters examined within this study that need to be controlled for the directed differentiation of pluripotent cells to pancreatic endoderm are retinoic acid, BMP, and FGF pathways and that fetal β cells can be rapidly induced through three sequential stages using five small molecules: RA, LDN3189, PD0325901, γ-XX, and A8301.

Discussion

Others have argued that biological complexity is a barrier to fulfilling the potential of biotechnology (Sadowski et al., 2016). A solution to this barrier was proposed to require large numbers of complex experiments, combined with sophisticated software and hardware (Sadowski et al., 2016). Here we describe a generally applicable method for extracting the critical process parameters for media components needed for the manufacturing of a cellular product. Although our focus was on properly controlling pathways during a differentiation event, the HD-DoE process itself is amenable to identifying broader aspects, which may impact cellular identity such as basal media formulations, oxygen dependencies, seeding density, mechanical forces, and even timing of differentiation events as examples. We refer to this process as High-Dimensional Design of Experiments (HD-DoE).

The HD-DoE method circumvents many of the limitations imposed by the hypothesis-driven scientific process. By relying on computer-based experimental designs, such as D-optimal designs, it is possible to extract maximal information per cost unit from the system interrogated. From the experimenters' perspective, it is possible to do virtual experimentation (ISPA) in a highly predictive, fully data-driven manner. This approach allows for hypothesis testing to occur freely within the known space covered by the design geometry after the experiment has been performed. The HD-DoE method is compliant with the basic principles of Quality-by-Design (QbD). We argue that the HD-DoE approach is applicable for defining optimal differentiation conditions for almost any recalcitrant combinatorial problem related to mammalian cell culture. Industrial applications in need of process understanding include production of advanced biologics, specialized cells, and complex tissues for regenerative medicine. Drawbacks for process optimization using current methodologies, where individual components are tested serially, are that process understanding is lacking and that critical process parameters remain undefined, which could result in the process remaining unstable. Thus, for industrial manufacture of human cells, the current reductionist approach—iteratively moving from hypothesis to hypothesis—remains slow and fails to provide a statistical basis upon which manufacture can succeed. For example, Loh et al. used serial testing covering a very large number of experimental conditions (>3,200) to derive a logic for endoderm induction from pluripotency (Loh et al., 2014). Comparatively, each modeling experiment performed using HD-DoE covers an experimental condition space of >4,000 possible conditions, and each condition is effectively arranged within the 12-dimensional design geometry. The labor intensity, inherent bias, and lack of design interactions coverage using serial experimentation makes machine-based experimental designs superior. Combining such with robotics and deep measurements provides the necessary systems-level interrogation and hereby enables the ISPA process.

Limitations of Study and Methodology

Although interrogations of human biology are amplified by implementation of high-dimensional experimentation, limitations exist. Generally, the ability to make good effector/response choices is highly dependent on developmental knowledge of the system being studied. The approach is not as useful for screening unknown effector inputs as compared with solving optimization problematics. Through our applications, we have understood that modeling experiments gives superior results when they are performed on homogeneous input cultures and that culture heterogeneity significantly blunts strength of response gene modeling. This is to be expected considering that multiple effector competencies would exist in the different cells, and therefore, the HD-DoE method is not ideal, and probably unable, to resolve problems at the subsequent stages for cultures that are highly heterogeneous. Also, interpretation of results relies deeply on the known behavior of chosen responses. Considering the fact that many genes are not limited to a single expression domain, but rather can be expressed in multiple tissues at different developmental stages, ambiguity follows during interpretation. In these cases, findings can only be resolved through understanding of the developmental system being studied. As an example, although PDX1 induction is usually chosen for determination of pancreatic fate, PDX1 expression is wider than in pancreas, occurring in the stomach, duodenum, and gall bladder. Only though examining multiple genes simultaneously and considering alternative PDX1-expressing fates was it possible to determine the best conditions for pancreatic fate. The most important tool for such analysis during ISPA is the dynamic profiling tool set that allows system inspection from the vantage point of, e.g., maximal expression of a given factor.

Although resting on large designs, the effector/response modeling always remains limited to factors being tested within the design space. Currently applied designs in this study were two-level D-optimal interaction designs, and therefore, non-linear responses within the concentration ranges tested would not be appropriately modeled. Inspecting the distribution of standardized residuals provides a means to detect non-linearity. Should a response gene display a response profile exceeding the capability of the second-degree polynomial fit method (Partial Least Squares), transformation of data is possible, and recommended, to achieve a more optimal fit. That said, the most critical limitation of the approach is that models cannot predict conditions that are not inherently embodied within the design space. Conclusions are also limited to the cellular stage tested (as cells typically change competency for input effectors upon differentiating). Consequently, creating robust forward differentiation processes for any specialized human cell using the HD-DoE method involves serial conduction of modeling experiments.

Mechanisms of DFE Specification

In consideration of these aforementioned methodology limitations we have demonstrated a novel protocol capable of rapidly converting pluripotent cells into a regionalized endodermal population (DFE). This endodermal population is competent to form dorsal pancreatic progenitors and undergoing endocrine conversion through a three-stage protocol relying on the use of only five small molecules. Whether this DFE population has an equivalent state during human development or whether differentiation of pluripotent cells to the DFE state occurs in absence of a gastrulation-like event is unknown. Considering that the developing embryo begins retinoic acid production during gastrulation (Ulven et al., 2000), and that retinoic acid patterns endoderm toward dorsal fates (Davenport et al., 2016), it is possible that locally produced retinoids pattern gastrulating cells during development. Because provision of the ALK5 inhibitor positively affects DFE generation, it is suggested that Nodal signaling does not occur in this protocol. Alternatively, since it is known that not all gut endoderm derives from gastrulating cells (Kwon et al., 2008, Mcdonald and Rossant, 2014) and that Nodal signaling is not needed for generating dorsal fates in more primitive developmental models (Rottinger et al., 2015), it is possible that dorsal endoderm arises from a non-migratory population adjacent to the node (Ulven et al., 2000). Given that plasticity is known to exist in undifferentiated endoderm (Kumar et al., 2003) the hypothesis that cell fate can be directed through positional cues is also plausible. Interestingly, many of the genes shown to be directly activated by RA in this study are known to be induced during Wnt/AA-mediated generation of definitive endoderm (most notably FOXA2 and CXCR4) and therefore potentially responding to RA stronger than they do to the WNT/AA inputs. NODAL induction of DE differs from AA-induced DE-induction (Chen et al., 2013), where NODAL is the physiologically relevant molecule expressed within the node. It is possible that the impact of AA in most currently used protocols creates ventral patterning due to unknown differences in the functional activity of these different morphogens. Others have argued that the descendant population of an AA-induction step is heterogeneous (Green et al., 2011), further substantiated by the use of AA at lower concentrations to attain mesodermal induction. Whether RA contributes early in AA induction of DE is not fully established. Future studies are needed to clarify to what extent APS-derived endoderm might be supported by trace levels of in-culture produced RA or if retinoids provided through serum during APS-type DE formation contributes to the endoderm induction (D'amour et al., 2005).

Dorsal versus Ventral Pancreas Induction

Here we have demonstrated that it is possible to obtain patterned endoderm directly from pluripotency and have shown that it is competent for pancreatic induction. In comparison with previously published pancreatic protocols, we have demonstrated that three previously described differentiation stages can be reduced into a single stage. As a basis for our interpretation, RNA-seq data demonstrated that DFE is of a dorsal character using a set of 13 previously described genes defining a dorsal identity (Jennings et al., 2017). When compared with a pancreatic directed differentiation protocol (Xie et al., 2013) we concluded that DE-derived PE adopted a pancreatic phenotype during the published “PGT” stage as evident through the expression of HHEX, HNF1β, RFX6, HNF4A, PDX1, PROM1, and PROX1 but had no discernible dorsal identity until the cultures were exposed to retinoic acid at the “PFG” stage. This implies that the dorsal phenotype previously attributed to this protocol (Jennings et al., 2017) was a result of redirecting the population toward a dorsal fate at stage 3 of the protocol. For the DFE, the dorsal identity carries forward to pancreas and this state is permissive for induction of all pancreatic lineages. Previous studies have demonstrated that DE-derived pancreatic progenitors suffer from stray hepatic fates controlled through BMP signaling (Mfopou et al., 2010), a phenomenon that could be attributed to the ventral pancreas having bipotential competency for liver induction (Angelo et al., 2012; Deutsch et al., 2001; Tremblay and Zaret, 2005; Bort et al., 2004) and that early lateral plate-derived BMPs instruct this precursor toward hepatic fates (Chung et al., 2008). DFE-derived PE showed a bipotential competence for antral stomach induction, and we propose it is controlled through SHH signaling. This observation is supported by developmental studies in chick (Hebrok et al., 1998). Regardless of the dorsal/ventral origin of the endodermal population, both DFE and APS-DE readily give rise to pancreatic endoderm capable of generating endocrine cells. As observed for APS-DE derived cells, the DFE-derived endocrine cells are functionally more similar to fetal β cells than to fully mature glucose-responsive β cells. Whether DFE-derived insulin-producing cells can undergo in vivo maturation, as shown for APS-DE-derived cells (Kroon et al., 2008), needs to be determined. From a diabetes cell therapy development perspective, functional maturation toward glucose-dependent insulin release in vitro is a desirable goal that can be addressed using the HD-DoE method.

Resource Availability

Lead Contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact Jan Jensen (jensenj2@ccf.org or jjensen@trailbio.com, Phone: + 001 216 445 0990/+001 216 479 9754).

Materials Availability

Not applicable.

Data and Code Availability

RNA-seq data generated in this study are available at the NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra).

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.