Induced pluripotent stem cells (iPS cells) come from the reprogramming of somatic cells via overexpression of four exogenous transcription factors Oct4, Sox2, Klf4, and c-Myc (OSKM) in somatic cells.^[1] A large number of experiments have demonstrated that the reprogramming process of somatic cells includes early and late phases.^[2–8] In the early phase, a sequence molecular events takes place; for instance, the transcriptional shift from differentiated status, the mesenchymal to epithelial transition, the downregulation of lineage-specific markers, and the upregulation of early pluripotency genes such as alkaline phosphatase (AP) and embryonic surface marker SSEA-1. In the late phase, the reprogramming process includes several genetic-related changes; for instance, the transgene silencing, telomere elongation, X chromosome reactivation, erasure of epigenetic memory, and the activation of endogenous pluripotency genes Nanog, Oct4 and Sox2 (NOS). It was also found that the early phase is a stochastic (or probabilistic) phase of genes activation after the somatic cell is induced by exogenous OSKM factors, and the late phase is a deterministic (or hierarchical) phase of gene activation.^[9–11]

To understand the behavior of induced pluripotency in the reprogramming of somatic cells, Muraro et al.^[12] proposed a kinetic model which includes the feedforward loop, cascade, and AND gate gene regulatory network motifs. MacArthur and Lemischka^[13] argued that the pluripotent state can be amenable to analysis using the tools of statistical mechanics and information theory, and the pluripotent state is a statistical property of stem cell populations. Morris et al.^[14] introduced mathematical approaches to help map the landscape between cell states during reprogramming, it was found that the exogenous OSKM expression can change potential barriers between cell states and facilitate the transitions between cell states. By using single cell transcript profiling and mathematical modeling, Chung et al.^[15] found that reprogrammed cells infected with exogenous OSKM follow two trajectories, and the stochastic phase is an ordered probabilistic process with independent gene specific dynamics.

The reprogramming process of somatic cells has two important properties: the first is the delay kinetic behavior, which is a successively delayed activation of downstream factors in the reprogramming process; and the second is the irreversible switch behavior, that is an irreversible switch from the transgene-dependent phase to the transgene-independent phase. Meanwhile, fluctuations are ubiquitous in various biological systems; for instance, the biochemical interaction processes,^[16–19] the gene regulatory networks,^[20–22] and the neural signal systems,^[23–30] etc. Although the stochastic models were proposed to understand the induced pluripotent stem cell generation,^[10] most previous investigations of somatic cell reprogramming process focused onto the dynamics of induced pluripotency and its behaviors captured, and the intrinsic fluctuation of key genes expression level has not been investigated in the somatic cell reprogramming process.

In this paper, based on the gene network motifs according to the timeline of molecular events taking place during induced pluripotency,^[12] a kinetic model of coherent feedforward transcription regulation loops is proposed to study the intrinsic fluctuations and the sensitivity of response to variation in parameters in somatic cell reprogramming process. First, the formulas of Fano factor for key gene expression level, a measure of the relative size of the intrinsic fluctuations of gene expression level, are analytically derived from the kinetic model of somatic cell reprogramming process around the steady-state. By employing our Fano factor formulas, the effects of parameters (including the activated ratio by upstream transcription factor, the self-activation ratio, and the self-degradation ratio) on intrinsic fluctuations are discussed. It can be seen that our theoretical results obtained by Langevin theory are coincident with those obtained by the Gillespie algorithm. Second, susceptibility is used to measure the sensitivity of response to variation in systemic parameters. Finally, we draw our conclusions.

Based on the feedforward transcription regulation loops in somatic cell reprogramming process,^[12] we hypothesized that the mechanisms of transcription regulation are governed by Hill function form. In the deterministic description, the time evolution of expression levels of three key genes in the coherent feedforward transcription regulation loops (see Fig. 1) can be given by following equations in the dimensionless

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Scheme of coherent feedforward transcription regulation loops in somatic cell reprogramming process. S denotes the exogenous OSKM (Oct4, Sox2, Klf4, and c-Myc) factors, X represents the earliest pluripotency genes (such as the AP), Y represents the pluripotency genes (such as the SSEA-1), and Z denotes the endogenous pluripotency genes (Nanog, Oct4, and Sox2).

In this kinetic model, parameters k_i ( , S₂, X, Y) represent the activated ratio by the upstream transcription factor, (i = X, Y, Z) are the self-activation ratio, (i = X, Y, Z) are the self-degradation ratio, m_i ( , S₂, X₁, X₂, Y₁, Y₂, Z) are the critical value of Hill function. For the sake of simplicity, the other parameter values are set as , , , , , , , , , , , , , , , and . Under those parameter values, the kinetic model equations (1)–(3) involved in comprising several feedforward loops can explain the observed properties (the delay kinetic and irreversible switch behaviors) of reprogramming induced pluripotency by exogenous OSKM factors.^[12,31]

To study the intrinsic fluctuation and susceptibility in the somatic cell reprogramming process, the Hill coefficient n is set as 1, then the steady-state ( , , ) of kinetic model Eqs. (1)–(3) is obtained accordingly

In the stochastic description, the kinetics of feedforward transcription regulation loop is given by following Langevin equations

The linearization of Langevin equations (7)–(9) around the steady-state ( , , ) leads to

The Fokker–Planck equation corresponding to Langevin equations (14)–(16) can be written as

To study the relationship between steady-state and parameters, a susceptibility, which measures the sensitivity of a response to variation in parameter ** , is defined as

Based on the timeline of the molecular events taking place during induced pluripotency, the stochastic kinetic model of coherent feedforward transcription regulation loops is proposed to study the intrinsic fluctuations in the somatic cell reprogramming. The Fano factor formulas of three key genes expression level were derived by using of Langevin theory in last section. It is well-known that the Fano factor is independent of volume, and the open systems at thermodynamic equilibrium obey Poisson statistics, and for those the Fano factor is always 1. The effects of different parameter values on the intrinsic fluctuations are discussed below.

3.1. Effects of the activated ratio by upstream transcription factor on intrinsic noises

In the kinetic model of somatic cell reprogramming process, there are four activated ratios by upstream transcription factors; i.e., the parameters , , k_X, and k_Y. The effects of those activation ratios on Fano factors are shown in Fig. 2.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Intrinsic noises of gene expression levels of key genes X, Y, and Z versus the activated ratios ( , , X, Y) by upstream transcription factors. The lines are the Fano factor from theoretical formulas Eqs. (18)–(20), and the dots are the results of Fano factor from the Gillespie algorithm with size parameter .

Figure 2(a) shows that the internal fluctuations of the expression level of key genes (X, Y, and Z) are decreased with the increasing of activated ratio by exogenous OSKM factors (S). However, the internal fluctuation of expression levels of three key genes are not observably changed with the increasing of activated ratios ( and k_Y) as shown in Fig. 2(b) and 2(d). It is interesting that there is a minimum for the internal fluctuation of expression levels of gene Y when the activated ratio k_X by upstream transcription factor X is around 0.25 as shown in Fig. 2(c).

3.2. Effects of the self-activation ratio and self-degradation ratio on intrinsic noises

There are three self-activation ratios, that is, for transcription factor X, for transcription factor Y, and for transcription factor Z, the effects of those self-activation ratios on Fano factor of X, Y, and Z are shown in Fig. 3.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Intrinsic noises of the expression levels of key genes X, Y, and Z versus the self-activation ratios (i = X, Y, Z). The lines are the Fano factor from theoretical formulas Eqs. (18)–(20), and the dots are the results of Fano factor from Gillespie algorithm with size parameter .

As the self-activation ratio (i = X, Y, Z) increases, Figure 3 shows that the internal fluctuations of gene expression levels ( and Z) mainly depend on itself activation ratio. It is found that the Fano factor approaches one for low values of self-activation ratio, reaches a maximum when the self-activation ratio is increased, and then decrease back to one for high values of self-activation ratio. However, with the increase of the self-activation ratio of X, the Fano factors of Y and Z decrease, and saturate to a plateau value; as shown in Fig. 3(a). These results show that when the self-activation rate is very high or very low, the intrinsic noise of key genes X, Y, and Z approach the bare Poissonian limit (i.e., that is, the somatic cell reprogramming process induced by exogenous OSKM factors is generally effective in attenuating intrinsic noises of three key genes.

In the kinetic model of somatic cell reprogramming process, there are three self-degradation ratio; i.e., the parameters , and . The effects of those self-degradation ratios on Fano factors are shown in Fig. 4.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Intrinsic noises of expression levels of key genes X, Y, and Z versus the self-degradation ratios (i = X, Y, Z). The lines are the Fano factor from theoretical formulas Eqs. (18)–(20), the dots are the results of Fano factor from Gillespie algorithm with size parameter .

With the increasing of self-degradation ratio (i = X, Y, Z), Figure 4 shows that the internal fluctuations of expression level of three key genes (X, Y, and Z) mainly depend on itself self-degradation ratio. It can be found that the Fano factor approaches one for low values of self-degradation ratio, reaches a maximum when self-degradation ratio is increased, and then decreases for high values of self-degradation ratio. However, with the increasing of self-degradation ratio of X, the Fano factor of Y increases, and saturates to a plateau value as shown in Fig. 4(a). When the self-degradation rate is very high or very low, these results show that the intrinsic noise of key genes X, Y, and Z approach the bare Poissonian limit (i.e., ); that is, the somatic cell reprogramming process is generally effective in attenuating intrinsic noises of three key genes.

In our kinetic model, S denotes the exogenous OSKM (Oct4, Sox2, Klf4, and c-Myc) factors, X represents the earliest pluripotency genes (such as the AP), Y represents the pluripotency genes (such as the SSEA-1), and Z denotes the endogenous pluripotency genes (Nanog, Oct4, and Sox2). The delayed reprogramming dynamics and low efficiency properties can be explained by the genes regulation network structure which is described by a cascade of genes regulation motifs.^[12] It has been demonstrated that the reprogramming process of somatic cells includes early and late phases,^[2–8] and a stochastic model can predicts that most or all cells are competent for the reprogramming process.^[10]

Although the early phase of reprogramming process is a stochastic phase when the somatic cell is induced by exogenous OSKM factors, Figure 2 shows that, with the increasing of activated ratio (e.g., , , k_X, and k_Y), the intrinsic fluctuation ( , the Fano factor) of downstream transcription factor expression level is not observably affected by the activation of upstream transcription factor. Figure 3 and 4 show that the intrinsic fluctuations (i.e., the Fano factor) of three key genes (i.e., the earliest pluripotency genes such as the AP, the pluripotency genes such as the SSEA-1, and the endogenous pluripotency genes (Nanog, Oct4, and Sox2)) expression level are mainly determined by itself-activation ratio (i = X, Y, Z) or itself-degradation ratio (i = X, Y, Z), and there is a critical itself-activation ratio or itself-degradation ratio at which the intrinsic fluctuation of expression level of transcription factor is maximum. These results imply that the intrinsic noise of upstream transcription factor can not be transmitted to the expression of downstream transcription factor in the cascade of genes regulation motifs, which might be the reason why the reprogramming process is slow (or time delay dynamics) and the reprogramming is inefficient during induced pluripotency by exogenous OSKM factors when compared with cell fusion nuclear reprogramming.

3.3. Sensitivity of response to variation in parameters

Susceptibility is used to measure the sensitivity of response to variation in systemic parameters. Under different activated ratio by upstream transcription factor, the variation of self-activation or self-degradation ratios to the sensitivity of steady-state is considered in following, respectively.

Under different values of activated ratio by exogenous OSKM factors, the sensitivities of steady-state response to the variation in the self-activation and self-degradation ratios are shown in Fig. 5. The susceptibility of increases first, reaches a maximum, and then decreases for large values of self-activation ratio and self-degradation ratio . It can be found that the magnitude of maximum is increased with the increasing of activated ratio by exogenous OSKM factors.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The susceptibility of steady value to variation in self-activation and self-degradation ratios with different values of activated ratio ( ) by upstream transcription factor.

Under different values of activated ratio k_X by upstream transcription factor X, the sensitivities of steady-state response to the variation in the self-activation and self-degradation ratios are shown in Fig. 6, the susceptibility of increases first, reaches a maximum, and then decreases for large values of self-activation ratio and self-degradation ratio . It can be found that the magnitude of maximum is increased as the activated ratio k_X is increased by upstream transcription factor X; however, the maximum of sensitivity will disappear when the activated ratio k_X is too large (e.g., ).

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. The susceptibility of steady value to variation in self-activation and self-degradation ratios with different values of activated ratio (kX) by upstream transcription factor.

Under different values of the activated ratio k_Y by the upstream transcription factor Y and different activated ratios by exogenous OSKM factors, the sensitivities of steady-state response to the variation in the self-activation and self-degradation ratios are shown in Fig. 7, respectively. The susceptibility of increases first, reaches a maximum, and then decreases for large values of the self-activation ratio and self-degradation ratio . It can be found that the magnitude of maximum is increased with as the activated ratio k_Y is increased by upstream transcription factor Y (see Figs. 7(a)–7(b)) and with the activated ratio by exogenous OSKM factors (see Figs. 7(c)–7(d)).

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. The susceptibility of steady value to variation in parameters the self-activation and self-degradation ratios with different values of activated ratios (kX and ) by upstream transcription factor.

Figure 5–7 show that the susceptibility of steady-state response to the variation in systemic parameters has a nonmonotonic behavior (seem as a resonant phenomenon). The nonmonotonic behavior of susceptibility of steady-state response to the variation in parameters implies that the reprogramming process of somatic cells can be triggered through forced expression of a set of transcription factors (i.e., the earliest pluripotency genes such as the AP, the pluripotency genes such as the SSEA-1, and the endogenous pluripotency genes (Nanog, Oct4, and Sox2)) under certain value of self-activation (or the self-degradation ratios) of three key genes, at which the susceptibility is maximum.

In this paper, a kinetic model of coherent feedforward transcription regulation loops is proposed to study the intrinsic fluctuations and the sensitivity of response to variation in parameters in somatic cell reprogramming process. The formulas of Fano factor for three key gene expression level are analytically derived from the kinetic model of the somatic cell reprogramming process around steady-state. Susceptibility is used to measure the sensitivity of response to variation in systemic parameters.

It is found that the internal fluctuation of the expression levels of three key genes are not observably changed with the increasing of activated ratios ( and k_Y), and there is a minimum for the internal fluctuation of the expression level of gene Y when the activated ratio k_X by the upstream transcription factor X is around 0.25. As the self-activation ratio (or self-degradation ratio) increases, the internal fluctuations of the expression level of three key genes mainly depend on the activation ratio and degradation ratio. The Fano factor approaches one for low values of the self-activation ratio (or self-degradation ratio), reaches a maximum when the self-activation ratio (or self-degradation ratio) is increased, and then decreases for high values of self-activation ratio (or self-degradation ratio). The theoretical results obtained by Langevin theory are coincident with those obtained by the Gillespie algorithm. Under a fixed activated ratio by the upstream transcription factor, with as the self-activation ratio (or self-degradation ratio) increases, the susceptibility of steady-state increases at first, reaches a maximum, and then decreases. The magnitude of maximum is increased with the increasing of the activated ratio by the upstream transcription factor.

The theoretical results of this paper imply that: (i) the dynamics of coherent feedforward transcription regulation loops^[12,31] can be used to explain the observed delay kinetics and irreversible switch behavior of reprogramming induced pluripotency by exogenous OSKM factors; (ii) the intrinsic noise of upstream transcription factor can not be transmitted to the expression of downstream transcription factor in the cascade of genes regulation motifs; (iii) the susceptibility of steady-state response to the variation in systemic parameters has a nonmonotonic behavior, and the reprogramming process of somatic cells might be triggered through forced expression of a set of transcription factors. In fact, under certain circumstances, it was found that the amplification of low-level fluctuations in transcriptional status may be sufficient to trigger reactivation of the core pluripotency switch and reprogramming to a pluripotent state.^[34] Consequently, our results can provide new insights into the roles of internal fluctuation and susceptibility in the dynamics of induced pluripotency and the behaviors captured of somatic cell reprogramming.