Providence St. Joseph Health Providence St. Joseph Health
Providence St. Joseph Health Digital Commons Providence St. Joseph Health Digital Commons
Articles, Abstracts, and Reports
2-28-2017
Cell population structure prior to bifurcation predicts e6ciency of Cell population structure prior to bifurcation predicts e6ciency of
directed differentiation in human induced pluripotent cells. directed differentiation in human induced pluripotent cells.
Rhishikesh Bargaje
Kalliopi Trachana
Martin N Shelton
Christopher S McGinnis
Joseph X Zhou
See next page for additional authors
Follow this and additional works at: https://digitalcommons.psjhealth.org/publications
Authors Authors
Rhishikesh Bargaje, Kalliopi Trachana, Martin N Shelton, Christopher S McGinnis, Joseph X Zhou, Cora
Chadick, Savannah Cook, Christopher Cavanaugh, Sui Huang, and Leroy Hood
Cell population structure prior to bifurcation predicts
efficiency of directed differentiation in human induced
pluripotent cells
Rhishikesh Bargaje
a,1
, Kalliopi Trachana
a,1
, Martin N. Shelton
a
, Christopher S. McGinnis
a
, Joseph X. Zhou
a
,
Cora Chadick
a
, Savannah Cook
b
, Christopher Cavanaugh
b
, Sui Huang
a,2
, and Leroy Hood
a,2
a
Institute for Systems Biology, Seattle, WA 98109; and
b
Institute for Stem Cell & Regenerative Medicine, University of Washington Medicine Research,
Seattle, WA 98109
Contributed by Leroy Hood, December 29, 2016 (sent for review December 7, 2016; reviewed by Alfonso E. Martinez Arias and Irving L. Weissman)
Steering the differentiation of induced pluripotent stem cells
(iPSCs) toward specific cell types is crucial for patient-specific
disease modeling and drug testing. This effort requires the capac-
ity to predict and control when and how multipotent progenitor
cells commit to the desired cell fate. Cell fate commitment repre-
sents a critical state transition or “tipping point” at which complex
systems undergo a sudden qualitative shift. To characterize such
transitions during iPSC to cardiomyocyte differentiation, we ana-
lyzed the gene expression patterns of 96 developmental genes at
single-cell resolution. We identified a bifurcation event early in the
trajectory when a primitive streak-like cell population segregated
into the mesodermal and endodermal lineages. Before this branch-
ing point, we could detect the signature of an imminent critical
transition: increase in cell heterogeneity and coordination of gene
expression. Correlation analysis of gene expression profiles at the
tipping point indicates transcription factors that drive the state
transition toward each alternative cell fate and their relationships
with specific phenotypic readouts. The latter helps us to facilitate
small molecule screening for differentiation efficiency. To this end,
we set up an analysis of cell population structure at the tipping
point after systematic variation of the protocol to bias the differ-
entiation toward mesodermal or endodermal cell lineage. We
were able to predict the proportion of cardiomyocytes many days
before cells manifest the differentiated phenotype. The analysis of
cell populations undergoing a critical state transition thus affords
a tool to forecast cell fate outcomes and can be used to optimize
differentiation protocols to obtain desired cell populations.
single-cell analysis
|
critical state transitions
|
iPSC to cardiomyocyte
differentiation
|
differentiation efficiency
|
prediction
T
he availability of human induced pluripotent stem cells
(iPSCs) with their potential to differentiate into virtually any
cell type creates unprecedented possibilities, not only to study
human development and disease but also to generate patient-
specific cells to determine personalized drug response (1, 2).
However, steering iPSCs efficiently into pure populations of a
specific cell type, such as cardiomyocytes, remains a challenge,
because the binary nature of cell fate decisions often causes the
“leakage” of cells into undesired lineages at each such decision
point. Additionally, optimizing established differentiation pro-
tocols for a specific genetic background (i.e., patient-specific
iPSC lines) to maximize differentiation efficiency is time-con-
suming because of the long time period (up to weeks) until cells
display the differentiated phenotype that informs about the
success of a differentiation protocol. Thus, it is paramount to
develop tools that not only reveal the critical regulators that
govern lineage-specific decision-making but at the same time,
also facilitate and shorten the optimization procedures for iPSC
differentiation protocols.
Toward this aim, longitudinal single-cell gene expression
analysis provides a new avenue to understand lineage commitment
in mouse and human pluripotent cells (3–5). Reconstructing
lineage trajectories at single-cell resolution captures cell fate
transitions in a large statistical ensemble of identical systems as
each individual progenitor cell in a differentiating population,
allowing us to dissect molecular and cellular patterns driving
lineage commitment. For instance, single-cell resolution analyses
have shown that cell types form discrete clusters when gene ex-
pression patterns are visualized in a low-dimensional space using,
for example, principle component analysis or t-distributed sto-
chastic neighbor-embedding plots (3, 6) (Fig. 1 A and B). This
pattern is consistent with the concept of attractors [i.e., stable
cell states of the gene regulatory network (GRN)], which corre-
spond to the valleys in Waddington’s “epigenetic landscape” (7).
A cell fate transition then corresponds to a switching between
distinct attractors via transient unstable states and can be analyzed
as coordinated shift of gene expression in a low-dimensional cell-
state space (8) (Fig. 1 A and B). This formalism enables us to study
universal patterns that underlie major transitions of GRN states
(hence transitions of cell states) independent of specific molecular
mechanisms, such as the specific structure of the regulatory net-
work that drives the transition. Such phenomenological analysis of
major state shifts in complex systems has been successfully used
for ecosystems and social systems (9, 10).
Specifically, we postulate that exit from pluripotency is not
simply a jump between attractors but instead, is initiated by the
gradual destabilization of the pluripotent stem cell attractor
triggered by exogenous signals (i.e., growth factors and modulators
Significance
Induced pluripotent stem cells (iPSCs) open new possibilities
for generating personalized disease models and drug testing.
However, iPSC differentiation to a specific cell type can take
weeks to complete, delaying the optimization process (maxi-
mize yield of desired cell types) for each patient’s iPSC. This
task can be accelerated if the destination cell type can be de-
termined early during cell lineage trajectory before cells man-
ifest the desired phenotype. Our results indicate such a
possibility: by quantifying the cell population structure during
a critical state transition, we identified key regulators of line-
ages commitment and predicted the percentage of desired cell
types for several protocol variations 2 wk in advance.
Author contributions: R.B., K.T., S.H., and L.H. designed research; R.B. and K.T. performed
research; J.X.Z., C. Chadick, S.C., and C. Cavanaugh contributed new reagents/analytic
tools; R.B., K.T., and C.S.M. analyzed data; and R.B., K.T., M.N.S., S.H., and L.H. wrote
the paper.
Reviewers: A.E.M.A., University of Cambridge; and I.L.W., Stanford University.
The authors declare no conflict of interest.
Freely available online through the PNAS open access option.
1
R.B. and K.T. contributed equally to this work.
2
To whom correspondence may be addressed. Email: sui.huang@systemsbiology.org or
Leroy.Hood@systemsbiology.org.
This article contains supporting information online at
www.pnas.org/lookup/suppl/doi:10.
1073/pnas.1621412114/-/DCSupplemental
.
www.pnas.org/cgi/doi/10.1073/pnas.1621412114 PNAS
|
February 28, 2017
|
vol. 114
|
no. 9
|
2271–2276
DEVELOPMENTAL
BIOLOGY
of signaling pathways). This response is akin to flattening of the
valley in the landscape, which facilitates exit from the attractor
state until, at a critical point, the pluripotency attractor suddenly
vanishes, providing access to two alternative cell fate attractors
(Fig. 1 A and B). A destabilization of an attractor until it vanishes
formally represents a bifurcation event (11, 12). The associated
sudden qualitative system changes are known to produce the
phenomenological signatures of a critical state transition (“tipping
point”) (9, 10). As the cell population approaches a tipping point,
we expect to observe two changes that can only be revealed by
analyzing gene expression patterns at single-cell resolution: (i)
increased cell population diversity because of the destabilization
of the attractor and diminished “attracting force” and simulta-
neously, (ii) a higher coordination of gene expression across the
cells as they move on a trajectory along which the attractor tran-
sition takes place (11, 12).
This framework of attractor destabilization can help us for-
malize how exogenous signals relate to the differentiation effi-
ciency that is usually measured as the percentage of the desired
cell type in the differentiated cell population. We hypothesized
that the signals conveyed by the treatment to cause iPSC dif-
ferentiation not only destabilize the attractor but also, bias the
destabilization toward a specific lineage. This bias would be
manifest before fate commitment. Thus, we examined by single-
cell analysis the population structure after the treatment but
before cell lineage commitment to determine if it can inform
about the future course of the differentiation trajectory. To
validate our hypotheses, we systematically varied the levels of
differentiation cues for cardiomyocyte differentiation and in-
vestigated how a range of signals affected differentiation effi-
ciency. We show that our analysis of the cell population structure
at the tipping point can help us forecast the preference by dif-
ferentiating iPSC for cardiac over other fate options (hence, to
predict the efficiency of the desired differentiation).
Results
We first monitored changes in transcript levels at single-cell
resolution during the first 6 days as cells exit pluripotency and
move toward the cardiomyocyte cell fate (
Fig. S1). Extensive
prior knowledge guided us to identify (i) intermediate cell states
at branch points of development, (ii) key transcriptional regu-
lators that control cell fate decisions, and (iii) instructive signals
that guide the differentiation process (2, 13–17). We used this
knowledge to select 96 gene markers for our study (
Dataset S1,
Table S1
). A standard method for induced pluripotent stem cell-
derived cardiomyocytes (iCMs) differentiation (Fig. 1C) con-
sisting of the sequential treatment of iPSCs with cytokines and
other molecules that induce cardiac mesoderm in vivo was used:
activin A (day 0), BMP4 (bone morphogenetic protein 4) com-
bined with a Wnt pathway activator (day 1), and a Wnt antag-
onist (day 3) mimicking, at least partially, the differentiation
signals that epiblast (E) cells are exposed to during heart de-
velopment in vivo (13, 17, 18). This widely used protocol yields
∼70% cardiomyocytes within 2 wk in culture (Materials and
Methods), although there is considerable variability depending
on the initial conditions (i.e., iPSC density plating) as well as
genetic background of the ES cell/iPSC line used (16).
To reconstruct the iPSC to iCM differentiation trajectory and
identify lineage branch points, we measured transcript expres-
sion of the selected genes in ∼1,900 individual differentiating
cells obtained during the first 6 days of differentiation (Fig. 1D,
Fig. S2, and Dataset S2). A major lineage branching took place
at day 3 when individual cells transitioned from a multipotent,
primitive streak (PS)-like progenitor state to either a more dif-
ferentiated mesodermal (M) state or an endodermal (En) state
as indicated by lineage-specific transcripts (Fig. 1E). This abrupt
disappearance of the progenitor state and its split into two gene
expression programs suggest a bifurcation in the dynamics of the
underlying GRN (12)—a critical state transition (9, 10).
The signature of a critical state transition that can be identi-
fied by single-cell resolution analysis of cell populations consists
of (i) a decrease of overall cell to cell correlation with respect to
gene expression and (ii) a concomitant increase in overall gene
to gene correlation across the cells (11). The first one manifests
an increase in cell diversity as the attractor destabilizing and
allowing access to new GRN states (lineage priming). However,
Fig. 1. Directed differentiation at single-cell resolution. (A) Theoretical and
(B) experimental framework to study cell differentiation as a transition be-
tween attractors. Before transition (time t0): Cells in state A (local minima in
a quasipotential landscape) are defined by a distinct GRN state (expressed and
nonexpressed genes are colored and gray, respectively). The state A attractor
manifests as either a dense cloud of points in a high-dimensional cell-state space
(as measured using single-cell qPCR) or a tight, uniform distribution of a single
gene/dimensi on (as measured by flow cytometry) as shown in B. The tipping point
(time t1): The attractor destabilizes (via changes in the quasipotential landscape),
and cells become primed toward a future attractor state(s). Cells in the poised state
A’ exhibit increased cell diversity, which manifests as a shift in the high-dimensional
state space or a wider distribution in a single dimension. Posttransition (time t2):
Stable states B and C emerge through the stabilization of mutually exclusive GRN
states that can be observed as two clouds occupying distinct positions in the high-
dimensional cell-state space or bimodal distribution of the marker gene as shown
in B.(C) Snapshot of the iPSC to iCM differentiation protocol used in this study.
Asterisks mark the perturbation time points (days 0, 1, and 3) that also correspond
to cell culture media exchanges. (D) Diffusion map (DM) of the iPSC to iCM dif-
ferentiation based on 1,934 single-cell gene expression vectors of 96 genes. Color
of each dot represents the day of collection during differentiation. Arrows indicate
the direction of cell-state trajectories. The dashed arrow points toward the un-
desired cell state. (E ) Dynamics of state-specific transcription factors. The violin
plots show the variability of gene expression (log2Ex) across each cell population
for five transcription factors: NANOG (stem cell marker), GSC (PS marker), MESP1
(posterior PS/cardiac mesoderm marker), GATA4 (mesoderm and endoderm
marker), and TBX2 (cardiac marker).
2272
|
www.pnas.org/cgi/doi/10.1073/pnas.1621412114 Bargaje et al.
the counterintuitive increase in overall gene to gene correlation
reveals a tight coordination of gene expression before the tran-
sition (Fig. 2 and
Fig. S3) (11). These changes can be summa-
rized by the critical transition index I
C
(t) computed for each
measured time point t (Fig. 2A and
Fig. S3), which is defined as
the ratio of the average of all pairs of gene to gene correlation
coefficients to the average of all pairs of cell to cell correlation
coefficients. Computing the I
C
(t) values (from day 0 to day 3)
revealed a significant increase as the differentiating cell pop-
ulation approached the M–En branch point, indicating a bi-
furcation (Fig. 2A and
Fig. S3). To show that the observed trend
was independent of the quality and quantity of genes or cells, we
calculated the Ic for randomly selected subsets of our dataset
(
Fig. S3). Cells were indistinguishable at day 0 (E state), and
there was no apparent correlation between pluripotency and
lineage-specific transcripts (Fig. 2B and
Dataset S1, Table S2)—
consistent with the theory that, in an attractor state, cell pop-
ulation diversity is mainly caused by symmetric fluctuations
around the “set point” caused by gene expression noise (7, 19).
Specifically, on destabilization of the E state triggered by the first
differentiation signal (activin A), cells diversified, and gene to
gene correlation between NANOG,anE state-specific marker,
and PS state-specific markers increased (Fig. 2C and
Dataset S1,
Table S3). Our data confirm previously reported interactions
between NANOG and the transcription factor EOMES (20, 21)
as a major regulatory interaction that drives exit from pluri-
potency toward PS state (
Fig. S4). At day 2, when the PS-like
cells are still uniform with respect to lineage-specific markers, we
observed a temporary decrease, still significant, in the critical
transition index (Fig. 2A), consistent with the PS state being a
distinct and observable, although transient, stabilized state. By
day 2.5, the value of I
C
increased again driven by the emergence
of correlations and anticorrelations in the expression of lineage-
specific transcription factors (Fig. 2D and
Dataset S1, Table S4).
After cells were committed to a specific lineage, cell-state vari-
ability (within each new subpopulation) decreased, thus lowering
I
c
for each individual day 3 cell subpopulation (Fig. 2A).
Combining the above findings with consensus clustering and
correlation analysis allowed us to build a comprehensive model
of early iPSC to iCM differentiation (Fig. 2E). Our data support
two distinct cell (sub)states after day 2 (
Fig. S5), which were
evident in the mutually exclusive expression of the fate-determining
transcription factors indicative of binary lineage branching (22).
The i dentified heterogeneity at this stage can be correlated with
distinct in vivo states during the anterior–posterior patterning
of the PS (
Figs. S6 and S7). In particular, SOX17 (23) and
HAND1 (24) appeared to display the familiar toggle switch-like
binary behavior that segregates the PS-like cells into two dis-
tinctly primed popu lations: i f HAND1 >> SOX17, cells were
primed toward M fates (posterior PS ); if HAND1 << SOX17,
they were primed to the En fate ( anterior PS). Similar obser-
vations have been reported by other single-cell studies for
mesoderm differentiation (3, 4). Howeve r, our ana lysis addi-
tionally revealed that the expression of the cell surface marker,
cKIT, correlated with this anterior vs. posterior PS specification
(
Fig. S6). Thus, we decided to investigate the distribution of the
cKIT protein expression phenotype and its ass ociation with
mesoderm–endoderm branching.
We found that, at the tipping point, cKIT protein expression
varied among individual cells, displaying the widest spread,
consistent with maximal cell to cell variability (Fig. 3A). Addi-
tionally, around the tipping point, the heterogeneous cell pop-
ulation transiently exhibited an M–En continuum, in which
individual cells expressed the molecular signature, indicating
priming toward either the desired cardiac (cKIT
−
and HAND1
+
/
SOX17
−
cells) or the undesired noncardiac (cKIT
+
and
HAND1
−
/SOX17
+
cells) fate. Singl e-cell gene expressio n pro-
filing of the extreme tails of the cKIT distribution (outliers),
cKIT
low
and cKIT
high
cells, mapped them to cell states primed
for the M and En lineages, respectively. Thus, information on
prospective fate is hidden in the bulk population distributions
Fig. 2. A critical transition signature for differentiation branch points.
(A) Time point-specific boxplots represent the distribution of I
C
(t) values
from 1,000 permutations of 25 randomly selected genes. After bifurcation,
we used cells that cluster as M lineage for day 3. The mean value corresponds
to the I
C
(t) value [X
(t)
= 96 genes × M cells]. **P value < 2e-10 for comparison
between the time points (Kolmogorov–Smirnov test and Wilcoxon rank sum
test). (B) Gene to gene (GxG) correlation plots for six lineage-specific transcription
factorsatday0(“in attractor”). The shade corresponds to the Pearson’s corre-
lation across all of the cells for each pairwise comparison, whereas the shape of
the data cloud shows the distribution of Pearson’s correlation across all of the
cells for each gene pair. (C and D) GxG correlation plots for six lineage-specific
transcription factors during two state transitions. We can observe distinct pat-
terns for individual genes, such as EOMES (important during E → PS but not PS →
M transition), or small regulatory circuits (i.e., day 2.5 shows anticorrelated
networks that are related with lineage segregation). (E) Early iPSC to iCM dif-
ferentiation model. Each cell state (stable or transitional) can be marked by
specific transcription factors.
Bargaje et al. PNAS
|
February 28, 2017
|
vol. 114
|
no. 9
|
2273
DEVELOPMENTAL
BIOLOGY