Epigenetic Landscapes Explain Partially Reprogrammed Cells and Identify Key Reprogramming Genes

A common metaphor for describing development is a rugged “epigenetic landscape” where cell fates are represented as attracting valleys resulting from a complex regulatory network. Here, we introduce a framework for explicitly constructing epigenetic landscapes that combines genomic data with techniques from spin-glass physics. Each cell fate is a dynamic attractor, yet cells can change fate in response to external signals. Our model suggests that partially reprogrammed cells are a natural consequence of high-dimensional landscapes, and predicts that partially reprogrammed cells should be hybrids that co-express genes from multiple cell fates. We verify this prediction by reanalyzing existing datasets. Our model reproduces known reprogramming protocols and identifies candidate transcription factors for reprogramming to novel cell fates, suggesting epigenetic landscapes are a powerful paradigm for understanding cellular identity.


I. ATTRACTOR NEURAL NETWORKS: ADDITIONAL DETAILS
This supplementary text gives a brief introduction to Hopfield neural networks 1,2 and how they can be adapted to study epigenetic landscapes. We begin by reviewing the basic principles underlying the original Hopfield neural network. We then show how to generalize this to continuous spins 3 as well as discrete spins with correlated cell fates 4 (projection method). For an in-depth introduction to neural networks, please see the beautiful book by Amit 5 .
A. Discrete, Standard Hopfield There are N genes and each gene i is either on or off, with the output denoted by S i = ±1. Alternatively, we could use the variables S = 1 2 (S + 1) = 1, 0 with the corresponding substitutions in all equations below. The input to a given gene i is denoted by the local field where J ij is the interaction between gene i and gene j and B i is the external (i.e interaction independent) bias of gene i. Both J ij and B i are assumed to be independent of S i . The landscape H is given by where in equation 3 we have introduce the order parameter for the overlap (dot product or "magnetization") of a spin configuration with a given cell fate µ as m µ and also introduced the cell fate bias b µ . The overlap is defined in terms of the cell fate vectors ξ µ i as: To prove that H is a Lypanov function (i.e. has stable equilibrium states and follows the standard definition of an "energy"), it is necessary to show that H is a decreasing function and bounded below. To do so, consider flipping a single S i . The resulting change in H is When we have symmetric interactions, J ij = J ji , this simplifies to To determine the sign of ∆H we need the relation between h i and ∆S i . For deterministic (stochastic) dynamics, as long as ∆S i and h i are always (usually) the same sign, we always (usually) have ∆H < 0. Therefore, any set of dynamics that stochastically matches the sign of ∆S i and h i will lead to H being a Lypanov function. This implies that any choice of dynamics leads to the same stable fixed points, but may give rise to different trajectories, limit cycles, and sizes of basins of attraction for fixed points, see Amit 5 section 2.2 and 3.5 for a detailed analysis. Therefore, in this paper we focus on predictions that are independent of the exact dynamics. This is equivalent to thinking about the stationary properties of the model.
We will follow the standard convention for neural networks and physics and implement Glauber dynamics which is an asynchronous, stochastic update rule. In this update scheme, at each time step, one gene is selected at random and probabilistically updated according to its local field with h i defined above (or equivalently h i = − ∂H ∂Si ) and t time measured in discrete updates. Also, β = 1/T is the inverse temperature and characterizes the slope of the sigmoid function. When β → ∞, the sigmoid approaches a deterministic step function, while when β → 0 each state is equally likely. Now we need to specify the gene interaction J ij and establish the global minima of the system. There are p cell fates and the state of gene i in cell fate µ is given by ξ µ i . The gene interaction is a correlation based interaction and in the standard Hopfield neural network it is defined as In the standard Hopfield network, the cell fates have two assumptions. First, each cell fate is assumed to on average be unbiased (i.e. equal number of positive and negative spins) and second every pair of cell fates is approximately orthogonal These two assumptions can be relaxed in extensions of the standard Hopfield neural network, see later sections for one example (the projection method) that can incorporate correlated cell fates. Now we can prove that each cell fate is a global minima of the landscape. For no external fields, the landscape can be written as: This can be rewritten in terms of the overlap as: Then as long as N is large compared to p, whenever we are in a given cell fate the energy is H = −N/2 and this is the lowest bound since m 2 ≤ 1. We have shown that for p N , H is a decreasing, bounded function and hence is a Lypanov function. When p and N are both large, a full replica calculation shows that H remains a Lypanov function 6 .
While we have established that the landscape is a Lypanov function, we also need to examine the dynamical stability of the cell fates and the existence of spurious attractors. In the absence of stochastic update noise (β → ∞), we can examine the signal-to-noise ratio of the cell fates. If a state is dynamically stable, one needs S i h i > 0. When the state is in a given cell fate (without loss of generality assume cell fate 1), we have that which can be broken into a signal term (first term) and noise term (second term) as follows: For large N , the signal term approaches 1. We can evaluate the noise term by recognizing that it is an unbiased sum of (N − 1)(p − 1) ≈ N p random steps, and therefore has mean 0 and standard deviation √ pN , giving us Therefore as long as N is much larger than p, every cell fate is a fixed point. This rough signal-to-noise argument can be made more rigorous by a spin-glass replica calculation 6 which finds that cell fates are stable (in the case β → ∞) as long as the ratio of p/N is less than 0.138.
Here is an intuitive argument of why the landscape must be rugged, which implies the scaling of stable states with N . From looking at small systems, a naive guess would be that the number of stable states should scale with the size of the state space 2 N . This scaling could be achieved if each minima occurred when a single TF state is turned on while all the other TFs are off. However, this implies that each minima is only marginally stable; any spin flip will move the state out of the minima. In order to have a basin of attraction, more TFs are needed to determine the minima. A simple error correction or redundancy could be implemented by using r redundant TFs, but this would require exponentially more states r N . Instead, stable states could be determined by overlapping sets of TFs, as in the Hopfield neural network. This form of error-correction leads to frustration and Gaussian noise between the stable states, hence the scaling of stable states with N and not 2 N .
An unavoidable consequence of the non-linearity (ruggedness) of the Hopfield network is that in addition to the desired attractors (the input cell fates), there are additional spurious, metastable, attractors. There are a variety of spurious attractors, but the most common are symmetric mixtures of odd states 2 , for example without loss of generality we can make a spurious state with the first three cell types, S spur = major (ξ 1 + ξ 2 + ξ 3 ), where major stands for majority vote (equivalently the sign function) at each spin. The most common spurious attractor are symmetric mixtures of 3 states (as in the example above). A signal-to-noise analysis can also be done to establish that these spurious attractors are stable attractors, but with a smaller basin of attraction than the input cell fates (see Amit 4.3 for details 5 ).

B. Continuous, Standard Hopfield
The previous section describes the basic ideas of Hopfield neural networks. Here, we show how discrete Hopfield neural networks can be considered a limiting case of continuous differential equations of gene expression. We start by defining continuous spins, σ i , that can take on real number between −1 and 1. For continuous dynamics, we must modify the dynamics of the corresponding local field. In particular, if the local field decays in time with a time constant τ i we have where the J ij are the same as in the discrete case and the spin σ i is related to the local field by some monotonic Now the landscape is given by where the first two terms are the same as in the discrete case while the third is the new term for continuous only. Taking derivatives with respect to time gives us we can relate the derivative of h i to the derivative σ i . Then using the fact that g i is monotonically increasing we can show that the change in H is always negative: The decrease in H along with the fact that H is bounded below, shows that we have a Lypanuv function. It is easy to see that every discrete stable point is also a stable point in the continuous model; however, the continuous Hopfield neural networks can have additional stable points.

C. Continuous Gene Expression
A popular approach to model gene interactions is based on the genetic toggle switch 7 and represents gene interactions by a Hill function. For now, we will use the general variable σ ∈ [σ min , σ max ].
In the most general case, we have that where the input h i is in the range [−∞, ∞] and the output σ i is in the range If we rescale every gene by its dynamic range and center the Hill function at zero, we get that σ = σ ∈ [−1, 1] and Using the function above for σ i = g i [h i ] allows one to relate continuous Hopfield neural networks to gene expression using Hill coefficients.

D. Discrete as Limit of Continuous
How can we relate the continuous model of gene expression to the previous discrete model? There are two limits. First, if we take the discrete time limit with the update time much greater than the input memory, we get Second, in the genetic toggle switch language, when the cooperativity is large n 1, then S i → ±1. This gives us a deterministic, discrete model of gene expression. If we introduce stochasticity through Glauber dynamics, we completely recover the discrete Ising model of gene expression.

E. Discrete, Projection Method
The standard Hopfield attractor neural network assumes that the "memories" (cell fates) have nearly no correlations amongst themselves. However, cell fates are highly correlated (see Figure S1). Therefore, instead of the standard Hopfield attractor neural networks, we will implement the projection method neural networks 4 .
The correlation between cell fate µ and ν is given by Now the inferred correlation-based, TF interaction matrix is Then the landscape can be rewritten as where in equation 26 we have introduced the projection order parameter a µ which is the orthogonal projection of a spin vector onto the subspace spanned by the stable cell fates A simple geometric picture illustrates that H makes each cell fate a global minimum of the landscape. An arbitrary vector can be rewritten in terms of its projection in the cell fate subspace and its orthogonal component δS i , Then, the distance of an arbitrary vector S to the cell fate subspace is given by ∆, which can be rewritten as This allows us to rewrite the stabilizing term of the landscape as This provides a very clear interpretation of the landscape as the global distance of an arbitrary vector S to the natural cell fate subspace 4 .
Again, let's examine the signal-to-noise of cell fates in the absence of stochastic update noise. If a state is dynamically stable, one needs S i h i > 0. When the state is a given cell fate (without loss of generality assume cell fate 1), we have that Therefore, the stability of cell fate 1 has no noise interference from the other cell fates, and we have that cell fates are stable up to p/N = 1.