Theoretical modelling discriminates the stochastic and deterministic hypothesis of cell reprogramming

Embryonic stem cell, which is able to divide indefinitely while maintaining pluripotency, is expected to be useful in numerous subfields of life sciences, including toxicology, disease mechanism researchI1 and regeneration medicineI2 I3 . In 2006, Yamanaka and his colleagues achieved artificially induced pluripotent stem (iPS) cells from mouse fibroblastsI4 , and later from adult human fibroblastsI41 .

Induction of iPS cell is in a very low efficiency ($\sim 0.05\%$) at firstI4 I41 , and there are several models to explain the low efficiencyI42 . One of these models, the stochastic hypothesis, suggests that most or all the cells are competent for reprogramming but limited of randomness, only a part of cells are able to be reprogrammed in one experiment. On the other hand, elite model presupposes that only a few cells in a culture are competent for reprogramming and the process is deterministic. Generally the process of cell reprogramming was regarded as a stochastic processI5 I6 at first. There are numerous experiments indicated that the process of cell reprogramming consists of an early stochastic phase and a late hierarchic phase. The low efficiency also implies there is stochasticity. From the point of control theory, this model may be elucidated as a core control part, which is multistable and transits from one steady state to another one stochastically, and a downstream signal transduction cascadeI6 . However, these two models are soon challenged by other research which shows that cell reprogramming is a deterministic processI7 and especially almost all the cells are able to be reprogrammed into iPS cells, which conflicts with the elite model. The debate of whether the nature of cell reprogramming process is stochastic or deterministic hasn’t ended yet, since the gene regulation networks of cell differentiation are still poorly understoodMacArthur_2009 , thus it’s hard to predict how to induce iPS cells efficiently. Although the biological experiments have yet reveal the nature of cell reprogramming, there are other approaches to investigate this question. Because the elite model supposes only a fraction of cells can be reprogrammed, and this paper only focus on single cell dynamics, we use the term ’deterministic hypothesis’ instead of ’elite model’ to represent the hypothesis that cell reprogramming is a deterministic process.

Our aim is to discriminate these two hypotheses, the stochastic hypothesis and deterministic hypothesis, especially in the way that the experiments can apply to. These two hypotheses, apparently, are different since one model suppose that not all the cells can be reprogrammed statistically while the other one does. However, this difference isn’t able to be tested so far by experiments, because the efficiency of reprogramming may not only determined by the nature of reprogramming but also by the other factors including the environment and the method itself. Therefore, we turned to the dynamics of reprogramming. From a perspective of nonlinear dynamics, different cell fates are able to be regarded as a steady state of a dynamical system, which is the core gene regulation network for cell fate decision. This concept was noted by C. Waddington in the 1940s, which by now is called Waddington’s epigenetic landscape Waddington . In the Waddington’s metaphor, the cell fate is like a ball and the landscape epitomizes the constraints of dynamical system, thus the cell fate can transit by stochastically ’jump’ from one valley to another and the probability of transition depends on the magnitude of noise and the height of the hill between two valleys (the energy barrier)MacArthur_2009 Enver .

Under this scenario, the stochastic hypothesis represents that the cell fate decision network is a traditional Waddington’s landscape, consisting numerous valleys and hills. The ball (cell fate), is able to stochastically transit from one valley to another. Here each valley represents a cell fate and reprogramming process is to trigger cell from one fate to another. On the other hand, in the deterministic hypothesis stochasticity plays a trivial role, there is no transition between two steady state and the cell fate transition is determined by the continuous alternation of the landscape (fig.1).

Therefore, a one key difference between these two models is the reprogramming time in stochastic hypothesis depends not only on the ’energy barrier’ between two steady state but also on the magnitude of noise while the reprogramming time in deterministic hypothesis depends only on the dynamics shifting the landscape. This idea is intuitively obvious. However, the quantitative analysis of the timescale required for reprogramming still lacks, especially the contribution of noise to the reprogramming time. If the noise works trivially for the total time for reprogramming, alternating noise magnitude may not lead to any significant change of reprogramming time.

To calculate the time required for reprogramming of two hypotheses, we want to use one unified model to compare the two hypotheses. According to previous studies on lineage specification, mutual inhibition paradigm is shown to be the most significant concept.Zon2008 Huang2011 For example, two opposing transcription factors, GATA1 and PU.1 plays dominant role in erythroid/megakaryocyte and myelomonocytic lineage determination respectively. GATA1 and PU.1 mutual inhibited each other, and also control their own expression, forming an self-activation loop (fig.2a) Graf2009 . We will show that this simple model fits both stochastic hypothesis and deterministic hypothesis thus it provides a unified framework for us to compare the two distinct hypotheses.

As previous shown, the core network which in responsible for the cell fate decision is able to be generally regarded as mutual inhibition Enver M11 : two opposing transcription factors $x$ and $y$ inhibited each other while activated themselves respectively. This topology is able to be described by the following ordinary differential equations:

These equations are based on Hill functions which are generally used in describing gene expression regulation.Uri2006 The first item represents the self-activation of each transcription factor; the second item represents the mutual inhibition; and the third item describes the unregulated degradation. One of the benefit to use this simple model is that it exists two types of bifurcations (fig.2). In the type I bifurcation, as the parameter grows the steady state divides into two and the system becomes bistable (supercritical bifurcation). In the type II bifurcation, as the parameter grows, one steady state disappears and the other two change continuously (subcritical bifurcation).

In the deterministic hypothesis, the stochastic transition contributes trivially, which means the steady state has a deterministic trajectory during bifurcation. In the stochastic hypothesis, both the stem cell and the differentiated cell steady states are intrinsically existed. The cell is able to ’jump’ from one steady state to another stochastically. To calculate the reprogramming time, it’s necessary to calculate the height of the ’energy barrier’ between two steady states. For example, consider a nonlinear dynamic system which has more than one steady states

which the driving force F exhibits in a nonlinear manner. If the dynamic system is of one dimension, there is always a potential function $U$ satisfied $U_{AB}=\int_{x_{A}}^{x_{B}}F(x)dx$, which means $U$ is the primitive of $F(x)$. In this situation, the transition probability between two different steady states $x_{A}$ and $x_{B}$ is:

where $U_{AS}$ is the $\Delta U$ between $x_{A}$ and the highest energy point between $x_{A}$ and $x_{B}$ (the energy barrier). $\varepsilon$ is the magnitude of noise.

However, when the dimension of the system is greater than 1, the upper theory may not be able to be applied since the high-dimensional, non-equilibrium systems generally are not gradient system, which means the following equation may not be held:

Because only the gradient part of the driving force $\textbf{F}(\textbf{x})$ determines the transition rate between two steady states, it’s ideal to decompose the driving force into two parts: a gradient of potential $-\nabla U$ and the remanent as following:

where D is the diffusion coefficient tensor and here accounting for the magnitude of noise.M3 There are numerous methods to decompose the driving forceM2 , here we use the decomposition based on the flux of the probability which the probability can be described as a diffusion equation:

The diffusion equation describes a conservation of probability (local change is due to net flux). Flux vector J is defined as $\textbf{J}=\textbf{F}P-\textbf{D}\cdot\frac{\partial P}{\partial\textbf{x}}$.M4 representing the speed of the flow of the probability in concentration space x.

In the steady state $\frac{\partial P}{\partial t}=0$, then $\nabla\cdot\textbf{J}(\textbf{x},t)=0$. If the system is an equilibrium system, $\nabla\cdot\textbf{J}(\textbf{x},t)=0$ always leads to $\textbf{J}=0$. However, for a nonequilibrium system, this condition does not mean that J have to vanish because the detailed balance condition is satisfied. Therefore, deviating from the definition of J, the flux in the steady state is defined as: $\textbf{J}_{ss}=\textbf{F}P_{ss}-\textbf{D}\cdot\frac{\partial P_{ss}}{\partial\textbf{x}}$, so:

where $U=-lnP_{ss}$. So $U$ reflects the potential and $\textbf{F}_{r}=\textbf{J}_{ss}/P_{ss}$.

The quasi-potential $U$ is able to be calculated by $U=-lnP_{ss}$ which $P_{ss}$ is calculated through a Fokker-Planck equation:

This equation was solved numerically by finite difference method with the second boundary condition.(fig.4)

As previous discussed, there are two types of bifurcation generally in the model (see Methods). The first type is illustrated in fig.2b, which the bifurcation is controlled by the inhibition coefficients $\beta$. This bifurcation could be elucidated from a biological perspective as a deterministic transition from one state (e.g. ESC state) to another state (e.g. stomatic cell state). The second type of bifurcation is illustrated in fig.2c, which the degradation rate $k$ controls the systems with 3 steady states to 2 steady states. This type may be elucidated as the ESC state goes to the differentiated state like stomatic cell through bifurcation or stochastic transition from one steady state to another and reprogramming required stochastic state transition but not only bifurcation.

Our aim is to discriminate the two hypotheses by their characteristic time for reprogramming. To simulate the experiments approach of reprogramming, we used a simple exponential equation to represent the effect of transformation of genes into somatic cell to reprogrammeI4 or add chemical small molecules to reprogrammeR1 , which trigger the control parameter’s change.

Here $\theta$ is the control parameter which in the type I bifurcation is $\beta$, and in the type II bifurcation, $k$. $\theta_{b}$ is the basal level of $\theta$, which means the final $\theta$ value; $\theta_{0}$ is the original level of $\theta$, representing the $\theta$ value of differentiated cell state; and $r$ represents how fast $\theta$ declines. When the system with $\theta_{0}$ reaches its steady state, $\theta$ begins to change as in Equation 10, then the time for reaching the new steady state is calculated.