An Algorithm for Cellular Reprogramming

TL;DR

This paper presents a data-driven dynamical model of cellular gene expression to optimize transcription factor use in cellular reprogramming, demonstrating the potential of mathematical and control theory approaches in biology.

Contribution

It introduces a novel dynamical modeling approach based on gene clustering and optimal control to identify effective transcription factors for reprogramming.

Findings

Validated several transcription factors for reprogramming

Demonstrated the effectiveness of dynamical models in biology

Highlighted potential for improved control of cellular processes

Abstract

The day we understand the time evolution of subcellular elements at a level of detail comparable to physical systems governed by Newton's laws of motion seems far away. Even so, quantitative approaches to cellular dynamics add to our understanding of cell biology, providing data-guided frameworks that allow us to develop better predictions about and methods for control over specific biological processes and system-wide cell behavior. In this paper we describe an approach to optimizing the use of transcription factors in the context of cellular reprogramming. We construct an approximate model for the natural evolution of a synchronized population of fibroblasts, based on data obtained by sampling the expression of some 22,083 genes at several times along the cell cycle. (These data are based on a colony of cells that have been cell cycle synchronized) In order to arrive at a model of moderate complexity, we cluster gene expression based on the division of the genome into topologically associating domains (TADs) and then model the dynamics of the expression levels of the TADs. Based on this dynamical model and known bioinformatics, we develop a methodology for identifying the transcription factors that are the most likely to be effective toward a specific cellular reprogramming task. The approach used is based on a device commonly used in optimal control. From this data-guided methodology, we identify a number of validated transcription factors used in reprogramming and/or natural differentiation. Our findings highlight the immense potential of dynamical models models, mathematics, and data guided methodologies for improving methods for control over biological processes.

Figures

Equations

$$x_{k+1}=A_k{x}_k+B{u}_k.$$

$$x_{j,k}:=\sum _{\begin{array}{c}i\text{ s.t.}\\ tad(i)=j\end{array}}{g}_{i,k}.$$

$$A_k$$

$$w_{i,m}:=a_m{s}_i{c}_{i,m},$$

$$b_{j,m}:=\sum _{\begin{array}{c}i\text{ s.t.}\\ tad(i)=j\end{array}}{w}_{i,m}$$

$$x_2$$

$$x_3$$

$$x_4$$

$$⋮$$

$$d=\Vert x_T-x_6(u_k)\Vert ,$$

$$\displaystyle\begin{array}[]{ll}\displaystyle\operatorname*{\text{minimize}}_{u_{k}}&\|x_{T}-x_{6}(u_{k})\|\\ \operatorname*{\text{subject to}}&\begin{cases}u_{m,k}\geq 0,&\text{if }m\in\hat{p}\\ u_{m,k}=0,&\text{if }m\not\in\hat{p}\end{cases}\end{array}$$

$$\displaystyle\begin{array}[]{ll}\displaystyle\operatorname*{\text{minimize}}_{\mathbf{\alpha}}&\|\mathbf{g}_{FIB,UM}-\mathbf{\alpha}\mathbf{g}_{FIB,GTEX}\|\end{array}$$

$$x_{k+1}=A_k{x}_k+B{u}_k.$$

$$x_4=A_3{A}_2{A}_1{x}_1+C\bar{u},$$

$$C:=A_3{A}_{2}B+A_{3}B+B$$

$$\underset{\bar{u}}{\min }\Vert x_T-A_3{A}_2{A}_1{x}_1-C\bar{u}\Vert .$$

$$x_T-A_3{A}_2{A}_1{x}_1\in \text{span}(C),$$

$$A_3{A}_2{A}_1{x}_1+C\bar{u}$$

$$C\bar{u}$$

$$\bar{u}$$

$$d=\bigg{\|}(I_{N}-CC^{\dagger})(x_{T}-A_{3}A_{2}A_{1}x_{1})\bigg{\|}$$

$$\mu :=d_0-d_{\ast }$$

$$e(u)=\sum _{k=0}^{K-1}{u}_k^T\cdot u_k={\mu }_2.$$

Full text

Significance Statement

Reprogramming the human genome toward any desirable state is within reach; application of select transcription factors drives cell types toward different lineages in many settings. We introduce the concept of data-guided control in building a universal algorithm for directly reprogramming any human cell type into any other type. Our algorithm is based on time series genome transcription and architecture data and known regulatory activities of transcription factors, with natural dimension reduction using genome architectural features. Our algorithm predicts known reprogramming factors, top candidates for new settings, and ideal timing for application of transcription factors. This framework can be used to develop strategies for tissue regeneration, cancer cell reprogramming, and control of dynamical systems beyond cell biology.

Introduction

In 1989, pioneering work by Weintraub et al. successfully reprogrammed human fibroblast cells to muscle cells via over-expression of transcription factor (TF) MYOD1, becoming the first to demonstrate that the natural course of cell development could be altered [1]. In 2007, Yamanaka et al. changed the paradigm further by successfully reprogramming human fibroblast cells to an embryonic-stem-cell-like state (induced pluripotent stem cells; iPSCs) using four TFs: POU5F1, SOX2, KLF4, MYC. This showed that a differentiated cell state could be reverted to a more pluripotent state [2].

These remarkable findings demonstrated that the genome is a system capable of being controlled via an external input of TFs. In this context, determining how to push the cell from one state to another is, at least conceptually, a classical problem of control theory [3]. The difficulty arises in the fact that the dynamics – and even proper representations of the cell state and inputs – are not well-defined in the context of cellular reprogramming. Nevertheless, it seems natural to treat reprogramming as a problem in control theory, with the final state being the desired reprogrammed cell. In this paper, we provide a control theoretic framework based on empirical data and demonstrate the potential of this framework to provide novel insights into cellular reprogramming [4, 5].

Our goal is to mathematically identify TFs that can directly reprogram human fibroblasts to a desired target cell type. As part of our methodology, we create a model for the natural dynamics of proliferating human fibroblasts. We couple data from bioinformatics with methods of mathematical control theory–a framework which we dub data-guided control (DGC). Using time series data and a natural dimension reduction through topologically associating domains (TADs), we capture the natural dynamics of the cell, including the cell cycle.

We use this model to determine a principled way to identify the best TFs for efficient reprogramming of a given cell type toward a desired target cell type. Previously, selection of TFs for reprogramming has been based largely on trial and error, typically relying on TF differential expression between cell types for initial predictions. Recently, Rackham et al. devised a predictive method based on differential expression as well as gene and protein network data [6]. Our approach is fundamentally different in that we take a dynamical systems point of view, opening avenues for investigating efficiency (probability of conversion), timing (when to introduce TFs), and optimality (minimizing the number of TFs and amount of input).

Using genomic transcription and architecture data, our method identifies TFs previously found to reprogram human fibroblasts into embryonic stem cell-like cells and reprogram fibroblasts into muscle cells. Our method also predicts TFs for conversion between human fibroblasts and many additional target cell types. In addition, we show the efficacy of using TADs for genome dimension reduction. Our analysis predicts the points in the cell cycle at which the insertion of TFs can most efficiently affect a desired change of cell state. Implicit in this approach is the notion of distance between cell types, which is measured in terms of the amount of transcriptional change required to transform one cell type into another. In this way, we are able to provide a comprehensive quantitative view of human cell types based on the respective distances between cell types.

Our framework separates into three parts:

Define the state. Use structure and function observations of the initial and target cell types’ genomes to define a comprehensive state representation. 2. 2.

Model the dynamics. Apply model identification methods to approximate the natural dynamics of the genome, from time series data. 3. 3.

Define and evaluate the inputs. Infer from bioinformatics (TF binding location and function) where TFs can influence the genome, then quantify controllability properties with respect to the target cell type.

The actual dynamics of the genome are undoubtedly very complicated, but as is often done in mathematical modeling studies, we use measurements to identify a linear approximation. This will take the form of a difference equation that is widely studied in the control systems literature, [7]

[TABLE]

In this case, the three items listed above correspond respectively to the value of the state $x_{k}$ at time $k$, the time dependent state transition matrix $A_{k}$, and the input matrix $B$ (along with the input function $u_{k}$).

Methods

Genome State Representation and Dimension Reduction: $\mathbf{x}_{k}$

The state representation $x$ in Eq. 1 is the foundation for any control system and is critical for controllability analysis. To fully represent the state of a cell, a high number of measurements would need to be taken, including gene expression, protein level, chromatin conformation, and epigenetic measurements. As an initial simplification, we assume that the gene expression profile is a sufficient representation of the cell state.

Gene expression for a given cell is dependent on a number of factors, including (but not limited to): cell type, cell cycle stage, circadian rhythm stage, and growth conditions. In order to best capture the natural fibroblast dynamics from population-level data, time series RNA-seq was performed on cells that were cell cycle and circadian rhythm synchronized in normal growth medium conditions (See SI). Prior to data collection, all cells were temporarily held in the first stage of the cell cycle, G0/G1, via serum starvation. Upon release into the cell cycle, the population was observed every $\Delta t=8$ hours (h) for 56 h, yielding 8 time points (at 0, 8, 16, $\dots$, 56 h). Let $g_{i,k}$ be the measured activity of gene $i=1,\dots,N$ at measurement time $k=1,\dots,8$, where $N$ is the total number of human genes observed (22,083). Analysis of cell-cycle marker genes indicated that the synchronized fibroblasts took between 32-40 h to complete one cell cycle post growth medium introduction, after which cells became largely unsynchronized (Fig. S1). Because of this, we define $K=5$ to be the total number of time points used for this model.

An obstacle to using $g$ to represent $x$ in a dynamical systems approach is the computational feasibility of studying a system with over 20,000 variables, necessitating a dimension reduction. Naïve dimension reductions such as partitioning the genome into 1 mega-base pair (Mb) bins ignores inherent structural organization of the genome and obscures important intricacies of finer resolutions. A comprehensive genome state representation should include aspects of both structure and function, and simultaneously have low enough dimension to be computationally reasonable. Along these lines, we propose a biologically inspired dimension reduction based on topologically associated domains (TADs).

The advent of genome-wide chromosome conformation capture (Hi-C) allowed for the studying of higher order chromatin structure and the subsequent discovery of TADs [8]. TADs are inherent structural units of chromosomes: contiguous segments of the 1-D genome for which empirical physical interactions can be observed [9]. Moreover, genes within a TAD tend to exhibit similar activity, and TAD boundaries have been found to be largely cell-type invariant [9, 10, 11]. TADs group structurally and functionally similar genes, serving as a natural dimension reduction that preserves important genomic properties. Fig. 1 depicts an overview of this concept. We computed TAD boundaries from Hi-C data via an algorithm that uses Fielder vector partitioning, described in Chen et al. (See SI) [12].

Let $tad(i):=j$ if gene $i$ is contained within TAD $j$. We define each state variable $x_{j,k}$ to be the expression level of TAD $j=1,\dots,\tilde{N}$ at time $k$, where $\tilde{N}=2,245$ is the total number of TADs that contain genes. Specifically, $x_{j,k}$ is defined as the sum of the expression levels of all genes within the TAD, measured in reads per kilobase of transcript per million (RPKM), i.e.

[TABLE]

The vector of all TAD activities at measurement $k$ is denoted with a single subscript $x_{k}\in\mathbb{R}^{\tilde{N}\times 1},k=1,\dots,K$.

State Transition Matrix: $\mathbf{A}_{k}$

Given the data we have, perhaps the most direct way to model the evolution of TAD activity level would be to assume a model of the form $x_{k+1}=x_{k}+y$, where $x_{k}$ and $x_{k+1}$ come from data, and $y$ is the vector difference of $x_{k+1}$ and $x_{k}$. However, the data could also be viewed in a different way. Taken over a full cycle, the average value of the expression level of the 2,245 TADs is known, within experimental error. Assuming that there is a function $f$ which maps $x_{k}$ to $x_{k+1}$, we can subtract the steady state average, $\bar{x}$, and focus on measuring the deviation from average as the cycle evolves. With this in mind, we have $f(x)=\bar{x}+A(x-\bar{x})$ where $A$ is allowed to depend on $x$’s location in the cell cycle. That is, we build a model for the variation of the cell cycle about $\bar{x}$. For the model to match data and capture variability over the cell cycle, we will need to have a different $A$ for each time step. Using the principle that $A$ should differ as little from the identity as possible, we let $A_{k}$ be the identity plus a rank one matrix chosen to match the data, for each time step $k$. In this case we have $x_{k+1}-\bar{x}=A_{k}(x_{k}-\bar{x})$.

We define a time dependent state transition matrix $A_{k}$.

[TABLE]

where $I_{\tilde{N}}$ is the $\tilde{N}\times\tilde{N}$ identity matrix.

Input Matrix and Input Signal: $\mathbf{B},\mathbf{u}_{k}$

With the natural TAD-level dynamics established in the context of our control Eq. 1, we turn our attention to quantifying methods for control.

A TF is a protein that can regulate a gene positively or negatively by binding to a specific DNA sequence near a gene and encouraging or discouraging transcription. This is accomplished, for example, by altering local chromatin conformation or by recruiting RNA polymerase II and other transcriptional machinery [13]. The degree to which a TF activates or represses a gene depends on the specific TF-gene interaction and most likely on a variety of nuclear subtleties and intricacies that are difficult to quantify. Let $w_{i,m}$ be the theoretical regulation weight of TF $m$ on gene $i$, where $w_{i,m}>0$ ($w_{i,m}<0$) if TF $m$ activates (represses) gene $i$, and $m=1,\dots,M$, where $M$ is the total number of well-characterized TFs. Weights that are bigger in absolute value, $|w_{i,m}|\gg 0$, indicate stronger transcriptional influence, and weights equal to zero, $w_{i,m}=0$, indicate no influence.

Extensive TF perturbation experiments would be needed to determine $w_{i,m}$ for each TF $m$ on each gene $i$. Instead, we propose an alternative (simplified) method to approximate $w_{i,m}$ from existing, publicly available data for TF binding sites, gene accessibility, and average activator/repressor activity. To determine the number of possible binding sites a TF $m$ recognizes near gene $i$, we scanned the reference genome for the locations of potential TF binding sites (TFBSs) (See SI). Position frequency matrices (PFMs), which give information on TF-DNA binding probability, were obtained for 547 TFs from a number of publicly available sources ($\therefore M=547$). Let $c_{i,m}$ be the number of TF $m$ TFBSs found within $\pm 5$kb of the transcriptional start site (TSS) of gene $i$ (Fig. S2). In our model, the magnitude of $w_{i,m}$ is proportional to $c_{i,m}$. False negatives would include distal TFBSs outside of the $\pm 5$kb window, while false positives would be erroneous TFBS matches.

Although many TFs can do both in the right circumstances, most TFs have tendency toward either activator or repressor activity [13]. That is, if TF $m$ is known to activate (repress) most genes, we can say with some confidence that TF $m$ is an activator (repressor), so $w_{i,m}\geq 0$ ($w_{i,m}\leq 0$) for all $i$. To determine a TF’s function, we performed a literature search for all 547 TFs and labeled 299 as activators and 124 as repressors (See SI). The remaining TFs were labeled unknown for lack of conclusive evidence for activator or repressor function. In the case of inconclusive evidence, the TF was evaluated as both an activator and a repressor in separate calculations. Here, we define $a_{m}$ as the activity of TF $m$, with 1 and -1 denoting activator and repressor, respectively, and the sign of $w_{i,m}$ will be determined by $a_{m}$.

TFBSs are cell-type invariant since they are based strictly on the linear genome. However, it is known that for a given cell type, certain areas of the genome may be opened or closed depending on epigenetic aspects. To capture cell type specific regulatory information, we obtained gene accessibility data through DNase-seq. DNase-seq extracts cell type specific chromatin accessibility information genome-wide by testing the genome’s sensitivity to the endonuclease DNase I, and sequencing the non-digested genome fragments. This data is used for our initial cell type to determine which genes are available to be controlled by TFs [14]. Here, we define $s_{i}$ be the DNase I sensitivity information (accessibility; open/close) of gene $i$ in the initial state, with 1 and 0 denoting accessible and inaccessible, respectively (See SI).

We approximate $w_{i,m}$ as

[TABLE]

so that the magnitude of influence is equal to the number of observed consensus motifs $c_{i,m}$, except when the gene is inaccessible ($s_{i}$ = 0) in which case $w_{i,m}=0$.

Since we are working off a TAD-dimensional model, our input matrix $B$ must match this dimension. Let $b_{m}$ be a 2,245-dimensional vector, where the $j^{th}$ component is

[TABLE]

and define a matrix B=\left[\begin{array}[]{cccc}b_{1}\;b_{2}\;\cdots b_{M}\end{array}\right].

The amount of control input is captured in $u_{k}$, which is a $\mathbb{R}^{M\times 1}$ vector representing the quantity of the external TFs we are inputting to the system (cell) at time $k$. This can be controlled by the researcher experimentally through manipulation of the TF concentration [15]. In this light, we restrict our analysis to $u_{k}\geq 0$ for all $k$, as TFs cannot be subtracted from the cell. $u_{m,k}$ is defined as the amount of TF $m$ to be added at time point $k$.

With all variables of our control Eq. 1 defined, we can now attempt to predict which TFs will most efficiently achieve cellular reprogramming from some $x_{I}$ (initial state; fibroblast in our setting) to $x_{T}$ (target state; any human cell type for which compatible RNA-seq data is available) through manipulation of $u_{k}$. An overview of our DGC framework is given in Fig. 2.

Selection of TFs

We consider different scenarios for the type of input regime. The first is assuming the input signal is constant $u_{1}=u_{k}=\bar{u}$, intended to mimic empirical regimes where TFs are given at a single time point. We will show that this method is theoretically inferior to inputting different TFs at different points in the cell cycle in a later section.

Eq. 1 has an explicit solution that is easily computed.

[TABLE]

Notice the expression for $x_{4}$ depends the input matrix $B$ and the input signal $u_{k}$.

If $x_{T}$ is a target condition, then the Euclidean distance $\|\cdot\|$ can be used to measure how close a state is to the target state, i.e.:

[TABLE]

where the notation $x_{6}(u_{k})$ is used to emphasize the dependence of $x_{6}$ on $u_{k}$. Considering all possible input signals, one can compute the optimal control that finds the minimum distance for a given initial and target cell type. Let $u_{*k}$ denote the optimal $u_{k}$ used to minimize $d$, and $d_{*}$ denotes this minimum distance value.

We note here that when determining which TFs can be used to reach $x_{T}$, it is often desirable and more experimentally feasible to minimize the number of distinct TFs given to the cells. Transfection of cells with multiple different TFs can lead to cell stress and death, and a lower efficiency of transfection overall. Moreover, many experimentally confirmed direct reprogramming regimes use $\leq$4 TFs to achieve reprogramming. For these reasons, we set all indices of $u_{k}$ equal to zero, except for indices corresponding to TFs that we choose.

We define $\hat{p}$ to be a set of positive integers that refer to the indices of $u_{k}$ (read: TFs) that are allowed to be non-zero (e.g. $\hat{p}=\{1,4,7\}$ refers to TFs 1, 4, and 7). Let $p$ be the number of elements in $\hat{p}$. Given a set of TFs ($\hat{p}$), we determine the quantity and timing of TF input ($u_{*k}$) that minimizes the difference ($d_{*}$) between the initial ($x_{I}$) and target ($x_{T}$) cell state. Mathematically, this can be written as

[TABLE]

We use MATLAB’s lsqnonneg function to solve Eq. 9, which gives $u_{*k}$ and $d_{*}$.

Let $d_{0}:=\|x_{T}-x_{I}\|$, be the distance between initial and target states with no external influence. We define the score $\mu:=d_{0}-d_{*}$, which can be interpreted as the distance progressed towards target. $\mu$ can be calculated for each $\hat{p}$ and sorted (high to low) to determine which TF or TF combination is the best candidate for direct reprogramming between $x_{I}$ and $x_{T}$.

Remark: Subsets of TFs were chosen for each calculation based on the following criteria: $\geq$10-fold expression increase in target state compared to initial state, and $\geq$4 RPKM in target state. These criteria are used to select differentially expressed TFs and TFs that are sufficiently active in the target state.

Results

Quantitative Measure Between Cell Types

In order to best utilize our algorithm to predict TFs for reprogramming, compatible data on target cell types must be collected. For this, we explore a number of publicly available databases where RNA-seq has been collected, along with RNA-seq data collected in our lab. The ENCODE Consortium has provided data on myotubes and embryonic stem cells (ESCs) (See SI) [16]. The GTEx portal provides RNA-seq data on a large variety of different human tissue types [17]. Although each GTEx experiment is performed on tissue samples, thus containing multiple different cell types, we use these data as more general cell state targets.

To give a numerical structure to cell type differences, conceptually similar to Waddington’s epigenetic landscape, we calculate $d_{0}$ between all cell types collected. Fig. 3A shows $d_{0}$ values for 32 tissue samples collected from the GTEx portal, along with ESC, myotube, and our fibroblast data (additional cell type $d_{0}$ values shown in SI). Warmer colors (red) denote further distances between cell types. GTEx RNA-seq data is scaled to keep total RPKM difference between time series fibroblast and GTEx fibroblast RNA-seq minimal (See SI).

TF Scores

To assess our method’s predictive power, a subset of target cell types are presented here that have either validated TF reprogramming methods or TFs highly associated with the target cell type. Additional predicted TFs for reprogramming are included in SI. We note that though experimentally validated TFs provide the best current standard for comparison, we believe experimental validation with our predicted TFs may provide more efficient reprogramming results. For all reprogramming regimes presented in this section, fibroblast is used as the initial cell type due to the availability of synchronized time series data, and all TFs are introduced at $k=1$ [9].

For conversion of fibroblast to myotubes, the top predicted single input TFs are MYOG and MYOD1, both of which are known to be crucial for myogenesis. While MYOD1 is the classic master regulator reprogramming TF for myotube conversion, activation of downstream factor MYOG is necessary for full conversion [18]. For fibroblast to ESC conversion, a number of TFs known to be necessary for pluripotency are predicted, including MYCN, ZFP42, NANOG, and SOX2 [19]. With the knowledge that no single TF has been shown to fully reprogram a fibroblast to an embryonic state, combinations of TFs are more informative for this analysis. The top scoring combination of 3 TFs is MYCN, NANOG, and POU5F1–three well-known markers for pluripotency [19]. Interestingly, POU5F1 scores poorly when input individually, but is within the top set of 3 TFs when used in combination with MYCN and NANOG. Left ventricle reprogramming includes TFs that are known to be necessary for natural differentiation in the top score for all 1-3 combinations. These include GATA4 (a known TF in fibroblast to cardiomyocyte reprogramming), HEY2, and IRX4 [20, 21, 22].

Time-dependent TF Addition

Fibroblast to ESC conversion was of particular interest in our analysis as this is a well-studied regime with a number of validated TFs (with a variety of reported efficiencies), and this conversion is promising for its regenerative medicine application. High scoring TFs yield many that are known markers for pluripotency, but the top combination of 3, MYCN, NANOG, and POU5F1, has not been used specifically together, to our knowledge.

Since our method incorporates dynamical RNA-seq data, analysis can be extended to determine the best time to input control for a given set of TFs. In our model, there are five possible input times: 0, 8, $\dots$, 32 h. We assume a TF continues to influence the system at a constant value once input until the final time (40 h). We restrict our analysis here to combinations of 3 TFs. This gives $5^{3}=125$ possible Time-dependent regimes to input the TFs; e.g. TF1, TF2, TF3 are input, respectively, at times 0,0,0, or 0,0,8, or 0,0,16, or $\dots$ or 32,32,24 or 32,32,32. Inputting a TF at time $k^{*}$ can be viewed mathematically as requiring $u_{m,k}=0$ for all $k<k^{*}$.

Time-dependent analysis of the top scoring ESC TFs reveals that scores vary widely depending on the time of input. MYCN and NANOG show a strong preference for input at the beginning of the cell cycle, while POU5F1 shows a slight preference for input towards the end of the cell cycle, with the highest score achieved when MYCN and NANOG are input at 0 h and POU5F1 is input at 32 h. Analysis on how the time of input control affects $\mu$ is shown in Fig. 3C. Time-dependent analysis was also conducted for the top combination of 3 TFs for fibroblast to left ventricle. This set includes GATA4, HEY2, and IRX4, all factors highly associated with the cardiac phenotype [20, 21, 22]. This analysis predicted that the best reprogramming results would occur if GATA4 is given immediately (0 h), with IRX4 and HEY2 given later (24 and 32 h, respectively).

Discussion

The results from this algorithm show promise in their prediction of known reprogramming TFs, and demonstrate the importance of including time series data for gene network dynamics. Time of input control has shown to have an impact on the end cell state, in line with what has been shown in natural differentiation [23].

While we believe that this is the best model currently available for predicting TFs for reprogramming, we are aware of its limitations and assumptions. TAD-based dimension reduction is based on the observation that genes within them correlate in expression over time, though we lack definitive proof of regulation by shared transcriptional machinery [9]. This assumption was deemed necessary for dimension reduction in the context of deriving transition matrix $A_{k}$. With finer time steps in RNA-seq data, the assumption may not be necessary for TF prediction, at the cost of increased computation time. Additionally, a 5kb window flanking the TSS of each gene was used to ensure that all potential regulators are found, at the cost of potential inclusion of false positive motifs.

GTEx data proved to be an invaluable resource for testing our algorithm, providing many target states for TF prediction. It is important to note that these data are collected from cadaver tissue samples; therefore the RNA-seq data is coming from a heterogeneous cell-type population and may be enriched for stress factors known to be elevated after death (e.g. HSF4). Ideally, RNA-seq data for target cells would be derived from a homogeneous population, with minimal experiment collection variables. Future work includes the extension of this DGC approach to other target cell types.

Although this program can score TFs relative to other TFs in a given reprogramming regime, it is difficult to predict a $\mu$ threshold that would guarantee conversion. Additionally, rigorous experimental testing will be required to validate these findings and determine how our $u$ vector translates to TF concentration. This is a product of the large number of assumptions that must be made to develop the initial framework for a reprogramming algorithm. With finer resolution in the time series gene expression, more subtle aspects of the genomic network may be observed, allowing for better prediction.

Our proposed data-guided control framework successfully identified known TFs for fibroblast to ESC and fibroblast to muscle cell reprogramming regimes. The framework rates individual TFs as well as sets of TFs. We employ a biologically-inspired dimension reduction via TADs, a natural partitioning of the genome. This comprehensive state representation was the foundation of our framework, and the success of our methods motivates further investigation of the importance of TADs as functional units to control the genome.

A dynamical systems view of the genome allows for analysis of timing, efficiency, and optimality in the context of reprogramming. Our framework is the first step toward this view. The successful implementation of time-varying reprogramming regimes would open new avenues for direct reprogramming. Experimental verification of predicted regimes and development of methods to identify optimal sets of TFs are planned for the near future. This template can be used to develop regimes for changing any cell into any other cell, for applications that include reprogramming cancer cells and controlling the immune system. Our DGC framework is well equipped for designing personalized cellular reprogramming regimes. Finally, this framework can serve as a general technique for investigating the controllability of networks strictly from data.

Methods and Materials

Hi-C and RNA-seq data were collected from cell cycle and circadian rhythm-synchronized proliferating human fibroblasts of normal karyotype. Data were collected every 8 h, spanning 56 h. Publicly available data was used for target cell types. Detailed materials and methods are provided in Chen et al. and in SI [9].

Acknowledgement

We thank Robert Oakes, Emily Crossette, and Sijia Liu for their critical reading of the manuscript and helpful discussions. We extend special thanks to James Gimlett and Srikanta Kumar at Defense Advanced Research Projects Agency (DARPA) for support and encouragement. This work is supported, in part, by the DARPA Biochronicity Program and the DARPA Deep-Purple and FunCC Program.

Supporting Information (SI)

Data

The fibroblast (FIB) data (Hi-C and RNA-seq) used for this application was originally collected and published in a paper by Chen et al. [9]. We refer the reader to this paper for a full description of technical protocols. Embryonic Stem Cell (ESC) and myotube (MT) data was downloaded from NCBI-GEO (GSE23316 ENCODE Caltech RNA-seq and GSE52529) [16]. 53 different tissue RNA-seq samples were downloaded from GTEx portal [17]. 51 different immune cell type RNA-seq samples were obtained from the BLUEPRINT Epigenome project [24].

Hi-C and Construction of TADs

We computed TAD boundaries from genome-wide chromosome conformation capture (Hi-C) data using an algorithm described in Chen et al. [12]. The algorithm was applied to averaged time series Hi-C data from proliferating human fibroblast (FIB) at 100 kilo-base pair (kb) resolution, which identified 2,562 TADs across all autosomal chromosomes (i.e. excluding Chromsomes X and Y). Of the 2,562 TADs, 317 contained no genes and were excluded from our analysis, leaving 2,245 TADs. These TADs ranged in size from a few hundred kb to several Mb, and contained on average 9-10 genes (standard deviation of 18 genes); one gene at minimum, and 249 genes maximum.

Construction of $\mathbf{B}$ Matrix

TF binding site position frequency matrix (PFM) information was obtained from Neph et al. and MotifDB, which is a collection of publicly available PFM databases such as, JASPAR, Jolma et al.  cispb_1.02, stamlab, hPDI, UniPROBE [25, 26]. TRANSFAC PFM information was included as well. Motif scanning of the human reference genome (hg19) was performed using FIMO of the MEME suite, in line with methods established by Neph et al. [25]. DNase-seq information for human fibroblasts was derived from ENCODE for fibroblast (GSM1014531). If a narrow peak is found within the $\pm$5kb of a gene TSS, the region is classified as open. TF function information was determined through an extensive literature search.

Scaling of RNA-seq

Due to differences in data collection procedures, the RNA-seq RPKM values obtained from the GTEx portal were of lower value, on average, compared to our fibroblast dataset, thus favoring repressor TFs for $\mu$ scoring. In order to account for this in our model, we scaled all GTEx RNA-seq data by a factor that solves the equation

[TABLE]

where $g_{FIB,UM}$ is the gene-level RNA-seq vector average of our fibroblast data, $g_{FIB,GTEX}$ is the gene-level RNA-seq vector of “Cells - Transformed fibroblasts” from the GTEx portal, and $\alpha$ is a scalar that solves this equation. For this data, $\alpha=2.6113$ and all GTEx data used as a target state was scaled by this factor.

Removal of MicroRNA

MicroRNA were removed from this analysis due to there high variance in RPKM values and unpredictable function.

TF Scores - Additional GTEx Data

For fibroblast to Adipose-Subcutaneous, the highest scoring factor is EBF1, a known maintainer of brown adipocyte identity, and a known promoter of adipogenesis in fibroblasts [27]. The 2nd highest scoring marker, PPARG, has been shown to be involved in adipose differentiation, and can be used individually to achieve reprogramming from fibroblasts [28]. Curiously, ATF3 is implicated here as being useful for adipocyte differentiation although its function has been shown to repress PPARG and stymie cell proliferation [29]. Upon further research using time dependent addition, ATF3 addition scores best when added towards the end of reprogramming process.

2 Brain tissue samples, Cerebellum and Hippocampus, both predict TFs necessary for natural differentiation. Interestingly, our algorithm select different TFs for each conversion, with factors linked specifically to each tissue. For Cerebellum, NEUROD1, has been shown to be required for granule cell differentiation, while ZIC1 and ZIC4 are both known to promote cellebular-specific neuronal function [30, 31]. The top scoring combination of 3 TFs are all similarly known to be important in neurogenesis (NEUROD1, ZBTB18, UNCX) [32, 33]. Hippocampus TF scoring includes FOXG1 as the top predicted factor, a factor specifically needed in hippocampus development. OLIG2, FOXG1, and GPD1 are the top scoring set for hippocampus reprogramming, all of which have been shown to been necessary for hippocampus function.

Colon TF scoring finds known differentiation factor in natural colon secretory linage development, ATOH1, as the highest scoring individual factor [34]. The top scoring combination of 2 TFs includes ATOH1 along with CDX2, another known factor necessary for full differentiation of colon cells, specifically small intestine maturation [35]. Liver cell reprogramming similarly finds known factors for differentiation in the top score of all 3 combinations: HNF4A, CUX2, PROX1 [36, 37, 38]. All factor play a role in correct development of hepatic progenitor cell-types and hepatic stem cells, the cell types just above in lineage differentiation.

TF Scores - BLUEPRINT Project Data

A number of immune cell types extracted from the BLUEPRINT Project revealed promising predicted TF results when fibroblast is used as the starting point [24]. $d_{0}$ values between cell types are shown in Fig. S4.

For fibroblast to macrophage direct reprogramming, a number of factors scoring highly in our algorithm are known to play a role in macrophage reprogramming or the differentiation. SPI1 (along with CEBPA) has been verified experimentally to reprogram fibroblasts into macrophage-like cells [39]. IKZF1 has been shown to be crucial for macrophage polarization via the IRF4/IRF5 pathway [40]. MYB has been shown to be crucial for the upstream cell type HSC [41].

For fibroblast to HSC direct reprogramming, the top scoring individual factor is highly associated with the target phenotype and has been shown to support HSC growth and regeneration [42]. ERG (in combination with GATA2, LMO2, RUNX1c, and SCL) is a confirmed reprogramming factor for fibroblast to HSC in mice [43].

For fibroblast to erythroblast reprogramming, ERG is a promising factor as it is required for the maintenance of the upstream cell type HSC [43]. NFIA is shown to promote the erythroid lineage from HSC differentiation [44].

Alternative Computation of $\mathbf{u}$

Below is an example of how $u$ can be computed without the constraint that $u_{k,m}\geq 0$. Assume $u_{k}:=\bar{u}$ is constant for all $t$. Then

[TABLE]

can be written as

[TABLE]

where

[TABLE]

We seek the control $\bar{u}$ that minimizes the distance between $x(3)$ and the target $x_{T}$:

[TABLE]

We can see that an exact solution exists if

[TABLE]

and is given by

[TABLE]

where $C^{\dagger}$ denotes the Moore-Penrose pseudoinverse of the matrix $C$, computed using the singular value decomposition of $C$. Even when Eqn (15) is not satisfied, it is well established that the control (18) solves (14) .

Define

[TABLE]

[TABLE]

$\mu$ can be used to compare between potential TFs for a defined initial state ($x_{I}$), target state ($x_{T}$), and TF number ($p$). The larger the value of $\mu$, the higher the relative score for its corresponding TF set.

We note that accurate TF predictions for some desired target cell type may not depend on minimizing distance alone, but also the amount of “energy" required for the system to reach $d_{*}$. We denote energy here with $\mu_{2}$ and define it as:

[TABLE]

$\mu_{2}$ is analogous to the amount of a TF that needs to be added to the system to achieve $d_{*}$. In the case where two different TFs achieve the same $\mu$ score, $\mu_{2}$ can be computed to decide the better candidate (i.e. lower $\mu_{2}$ is better TF).

SI Figures