{"id":478514,"date":"2024-01-05T19:00:00","date_gmt":"2024-01-06T00:00:00","guid":{"rendered":"https:\/\/platohealth.ai\/epigenetic-oct4-regulatory-network-stochastic-analysis-of-cellular-reprogramming-npj-systems-biology-and-applications\/"},"modified":"2024-01-07T02:17:24","modified_gmt":"2024-01-07T07:17:24","slug":"epigenetic-oct4-regulatory-network-stochastic-analysis-of-cellular-reprogramming-npj-systems-biology-and-applications","status":"publish","type":"post","link":"https:\/\/platohealth.ai\/epigenetic-oct4-regulatory-network-stochastic-analysis-of-cellular-reprogramming-npj-systems-biology-and-applications\/","title":{"rendered":"Epigenetic OCT4 regulatory network: stochastic analysis of cellular reprogramming – npj Systems Biology and Applications","gt_translate_keys":[{"key":"rendered","format":"text"}]},"content":{"rendered":"
In this section, we first introduce the model of the epigenetic OCT4 gene regulatory network developed in this paper. We then provide detailed explanations of the results obtained from our computational analysis. Specifically, we first studied how DNA methylation affects the dynamics of cellular differentiation, which we captured in our model by the progressive inactivation of the OCT4 gene starting from a high OCT4 level, corresponding to the pluripotent state49<\/a>,50<\/a><\/sup>. We then compared different reprogramming approaches by studying the impact of overexpressing OCT4 alone and in combination with TET1 and JMJD2 on the dynamics of OCT4 gene reactivation. Specifically, we examined process efficiency and latency variability for each approach, as these metrics have been commonly employed in published experimental studies to assess the success of cellular reprogramming21<\/a>,22<\/a>,23<\/a>,24<\/a>,31<\/a><\/sup>. To this end, we performed a computational study of the temporal trajectories of the system by simulating the reactions associated with the Epi OCT4 GRN (Fig. 2<\/a>) with Gillespie\u2019s Stochastic Simulation Algorithm (SSA)43<\/a><\/sup>.<\/p>\n Our model includes one TF gene, OCT4, and two chromatin modifier genes, TET1 and JMJD2. We consider OCT4 as the only TF in our model because it has been shown that overexpression of OCT4 alone leads to reprogramming15<\/a>,19<\/a>,35<\/a>,36<\/a>,37<\/a><\/sup>. The chromatin modifiers JMJD2 and TET1 are enzymes that catalyze the erasure of histone modification H3K9me3 and DNA methylation, respectively26<\/a>,38<\/a>,39<\/a>,40<\/a>,41<\/a>,42<\/a><\/sup>, two chromatin modifications associated with compacted chromatin state and gene silencing51<\/a><\/sup>. Specifically, while JMJD2 directly erases H3K9me326<\/a>,38<\/a>,39<\/a>,40<\/a><\/sup>, TET1 recognizes CpGme dinucleotides and converts methylated CpG to carboxylcytosine through multiple intermediate forms41<\/a>,42<\/a><\/sup>, none of which is recognized by DNMT1, the enzyme responsible for copying the CpGme pattern on the nascent DNA strand during DNA replication51<\/a><\/sup>. Transcription factor OCT4 recruits writers of H3K4me3 to its own gene52<\/a><\/sup> and to the JMJD253<\/a><\/sup> and TET154<\/a><\/sup> genes. This leads to a model in which OCT4 self-activates and also activates TET1 and JMJD2 by recruiting writers of activating chromatin modifications, while TET1 and JMJD2 self-activate and mutually activate each other and OCT4 by erasing repressive chromatin modifications. We call this system the epigenetic OCT4 gene regulatory network (Epi OCT4 GRN).<\/p>\n The chromatin modifiers TET1 and JMJD2 act on the chromatin modification circuit within each gene, which has been developed in ref. 44<\/a><\/sup>. This circuit includes H3K9 methylation (H3K9me3), DNA methylation (CpGme), H3K4 methylation\/acetylation (H3K4me3\/ac), and their known interactions. The first two modifications are associated with a repressed gene state51<\/a><\/sup>, while H3K4me3\/ac is associated with an active gene state26<\/a>,55<\/a><\/sup>.<\/p>\n We next describe the chromatin modification circuit and the model of gene expression (see ref. 44<\/a><\/sup> for details). In terms of species, in this model we have D (nucleosome with DNA wrapped around it), DA<\/sup> (nucleosome with H3K4me3\/ac), ({{{{rm{D}}}}}_{2}^{{{{rm{R}}}}})<\/span> (nucleosome with H3K9me3), ({{{{rm{D}}}}}_{1}^{{{{rm{R}}}}})<\/span> (nucleosome with CpGme), and ({{{{rm{D}}}}}_{{{{rm{12}}}}}^{{{{rm{R}}}}})<\/span> (nucleosome with both H3K9me3 and CpGme) (Fig. 1<\/a>a). Then, the expression rate of each gene will be determined by the number of nucleosomes with activating (DA<\/sup>) or repressive (({{{{rm{D}}}}}_{1}^{{{{rm{R}}}}},{{{{rm{D}}}}}_{2}^{{{{rm{R}}}}},{{{{rm{D}}}}}_{12}^{{{{rm{R}}}}})<\/span>) chromatin modifications. In terms of molecular interactions, both histone modifications and DNA methylation can be de novo established (process encapsulated in reactions \u24ea<\/span>, \u2460<\/span>, and \u2467<\/span>). Then, histone modifications can enhance the establishment of marks of the same kind to nearby nucleosomes via a read-write mechanism, generating auto-catalytic reactions (encapsulated in \u2461<\/span>, \u2462<\/span>). Analogously, repressive histone modifications enhance the establishment of DNA methylation, and vice versa, generating cross-catalytic reactions (encapsulated in \u246e, \u246f). Finally, each modification can be passively removed through dilution, due to DNA replication (reactions \u2463<\/span>, \u2464<\/span>, and \u2468<\/span>), or through the action of eraser enzymes (basal erasure) (reactions \u2465<\/span>, \u2466<\/span>, and \u2469). These erasers can be also recruited by the opposite modifications (recruited erasure), that is, repressive modifications recruit activating modification\u2019s erasers and vice versa (reactions \u246a, \u246b, \u246c, and \u246d). In this model, the rate of the processes described above for H3K9me3 (DNA methylation) is assumed not to change if the other repressive mark is present on the same nucleosome. All the reactions described above are collected in the list of Fig. 1<\/a>b, in which the reactions involving TET1 and JMJD2 are shaded in yellow and pink respectively. A diagram of the chromatin modification circuit corresponding to the reactions in Fig. 1<\/a>b is provided in Fig. 1<\/a>c. In this system, the transcriptional self-activation is modeled as a Hill function with cooperativity 1. Specifically, ({k}_{W}^{A})<\/span> (Fig. 1<\/a>b) is a monotonically increasing function of the abundance of X (X<\/i>), which can be written as ({k}_{W}^{A}={tilde{k}}_{W}^{A}(X\/{K}_{A})\/(1+(X\/{K}_{A})))<\/span>, in which ({tilde{k}}_{W}^{A})<\/span> is a coefficient that does not depend on X<\/i>, and K<\/i>A<\/i><\/sub> is the dissociation constant of the binding reaction between X and DNA44<\/a><\/sup>.<\/p>\n a<\/b> Nucleosome modifications considered in the model. b<\/b> Reactions associated with the chromatin modification circuit of each gene X, with X\u2009=\u2009O (OCT4), T (TET1), J (JMJD2). Here, reactions associated with activating marks, H3K9me3 and DNA methylation are enclosed in green boxes, pink boxes, and yellow boxes, respectively. Dark shades are associated with reactions describing the establishment of the modifications and light shades are associated with reactions describing the erasure of the modifications. Furthermore, shaded boxes enclose reactions involving J (pink) and T (yellow). Finally, each reaction rate constant is subscripted with W<\/i>, M<\/i>, or E<\/i> to indicate the association of its corresponding reaction with the writing (establishment), maintenance (auto\/cross-catalytic reactions), or erasure of a chromatin modification. c<\/b> Diagram representing the chromatin modification circuit for each gene X. d<\/b> Reactions associated with the gene expression model. Specifically, reactions associated with production (dark gray), dilution\/degradation (light gray), and artificial overexpression u<\/i>x<\/i><\/sub> (blue) of the gene product X. The reactions are described in the \u201cModel of the epigenetic OCT4 gene regulatory network\u201d subsection. e<\/b> Definitions and interpretations of r<\/i>, \u03b5<\/i>d<\/i><\/sub>, \u03b5<\/i>e<\/i><\/sub>, ({varepsilon }^{{prime} })<\/span>, (tilde{mu })<\/span>, and ({tilde{mu }}^{{prime} })<\/span>.<\/p>\n<\/div>\n<\/div>\n In our model, transcription is allowed only by DA<\/sup>, and transcription and translation are lumped together (reaction \u2470). The transcription by D is considered negligible because it is assumed that transcription of D by RNA polymerase II occurs concurrently with H3K4me3 deposition (i.e., conversion of D to DA<\/sup>), as observed in ref. 56<\/a><\/sup>. Furthermore, the gene product X is subject to dilution due to cell division and degradation (reaction \u2471). Finally, the production of X can be artificially increased through overexpression (reaction \u2472). A diagram representing these reactions is shown in Fig. 1<\/a>d. Then, based on the interactions among OCT4, TET1, and JMJD2, we can wire the three chromatin modification circuits to obtain the Epi OCT4 GRN circuit that we analyze in this paper (Fig. 2<\/a>). In the single gene\u2019s chromatin modification circuit proposed in ref. 44<\/a><\/sup> the concentrations of TET1 and JMJD2 impact the rate constants as they do here. However, these concentrations are constant parameters themselves. By contrast, in the Epi OCT4 GRN introduced here the concentrations of TET1 and JMJD2 are state variables and, as such, they vary dynamically under the effect of each other concentrations and the concentration of OCT4.<\/p>\n Diagram of the Epi OCT4 GRN, in which OCT4 self-activates and also activates TET1 and JMJD2 by recruiting writers of activating chromatin modifications on all genes (black arrows), while TET1 and JMJD2 self-activate and mutually activate each other and OCT4 by recruiting erasers for repressive chromatin modifications on all genes (pink and yellow arrows, respectively). For simplicity of illustration, in each gene\u2019s chromatin modification circuit we used gray for the solid arrows indicating the establishment and erasure and we did not represent the dashed arrows indicating recruitment and catalysis.<\/p>\n<\/div>\n<\/div>\n Let D<\/i>tot<\/i><\/sub>\u2009=\u2009Dtot<\/sub>\/\u03a9, in which Dtot<\/sub> represents the total number of modifiable nucleosomes within the gene of interest and \u03a9 represents the reaction volume, and the normalized time (tau =t{k}_{M}^{A}{D}_{tot})<\/span>. Let us define the dimensionless parameters (alpha ={k}_{M}\/{k}_{M}^{A})<\/span>, (bar{alpha }={bar{k}}_{M}\/{k}_{M}^{A})<\/span>, and ({alpha }^{{prime} }={k}_{M}^{{prime} }\/{k}_{M}^{A})<\/span>, in which \u03b1<\/i> represents the normalized rate constant of the auto-catalytic reaction for repressive histone modifications, and (bar{alpha })<\/span> and ({alpha }^{{prime} })<\/span> represent the normalized rate constants of the cross-catalytic reactions between repressive histone modifications and DNA methylation. Let us also introduce:<\/p>\n $$r=frac{{alpha }^{{prime} }}{alpha }=frac{{k}_{M}^{{prime} }}{{k}_{M}},$$<\/span><\/p>\n \n (1)\n <\/p>\n<\/div>\n that is, the ratio between the rate at which repressive histone modifications enhance the establishment of DNA methylation through cross-catalytic reactions (({k}_{M}^{{prime} })<\/span>) and the rate at which repressive histone modifications enhance their own establishment through auto-catalytic reactions (k<\/i>M<\/i><\/sub>). Furthermore, we use \u03b7<\/i> to represent the efficiency of the maintenance process of DNA methylation by DNMT126<\/a><\/sup>. The expression of \u03b7<\/i> was derived in ref. 44<\/a><\/sup> by accounting for the dilution of DNA methylation due to DNA replication and the maintenance process in which DNMT1 replicates CpG methylation on the newly synthesized DNA strand based on the pattern of the mother strand26<\/a>,57<\/a><\/sup>. The expression of \u03b7<\/i> is given by (eta ={delta }^{{prime} }\/delta)<\/span>, in which \u03b4<\/i> represents the rate constant of the passive erasure through dilution and ({delta }^{{prime} })<\/span> represents the effective passive erasure rate constant obtained from the balance between the dilution and the maintenance process. In particular, \u03b7<\/i>\u2009=\u20091 if DNMT1 is completely absent (no maintenance) and \u03b7<\/i>\u2009=\u20090 if the maintenance process is 100% efficient. Now, we define:<\/p>\n $$tilde{mu }=frac{{tilde{k}}_{E}^{R}{D}_{tot}}{{k}_{E}^{A}},,,{tilde{mu }}^{{prime} }=frac{{tilde{{k}^{{prime} * }}}_{T}{D}_{tot}}{{k}_{E}^{A}},,,{varepsilon }_{d}=frac{delta }{{k}_{M}^{A}{D}_{tot}},,,{varepsilon }_{e}=frac{{bar{k}}_{E}^{A}}{{k}_{M}^{A}{D}_{tot}},,,{varepsilon }^{{prime} }=frac{{k}_{E}^{A}}{{k}_{M}^{A}},$$<\/span><\/p>\n \n (2)\n <\/p>\n<\/div>\n with (tilde{beta }=O(1))<\/span> and (tilde{b}=O(1))<\/span> such that (({hat{k}}_{T}^{{prime} }{D}_{tot})\/{bar{k}}_{E}^{A}=tilde{beta }{tilde{mu }}^{{prime} })<\/span> and (({hat{k}}_{E}^{R}{D}_{tot})\/{bar{k}}_{E}^{A}=tilde{b}tilde{mu })<\/span>, respectively. Specifically, (tilde{mu })<\/span> is a dimensionless parameter quantifying the asymmetry between the erasure rates of repressive and activating histone modifications, while ({tilde{mu }}^{{prime} })<\/span> is a dimensionless parameter quantifying the asymmetry between the erasure rates of DNA methylation and activating histone modifications. Furthermore, ({varepsilon }_{d}=delta \/{k}_{M}^{A}{D}_{tot})<\/span> is the normalized rate constant associated with dilution due to DNA replication, and, since ({hat{k}}_{E}^{R}\/{k}_{M}^{A}=tilde{b}tilde{mu }{varepsilon }_{e})<\/span>, ({hat{k}}_{T}^{{prime} }\/{k}_{M}^{A}=tilde{beta }{tilde{mu }}^{{prime} }{varepsilon }_{e})<\/span>, ({tilde{k}}_{E}^{R}{D}_{tot}\/{k}_{M}^{A}=tilde{mu }{varepsilon }^{{prime} })<\/span> and ({tilde{k}}_{T}^{{prime} * }{D}_{tot}\/{k}_{M}^{A}={tilde{mu }}^{{prime} }{varepsilon }^{{prime} })<\/span>, the dimensionless parameter \u03b5<\/i>e<\/i><\/sub> (({varepsilon }^{{prime} })<\/span>) scales the ratio between the rate of the basal erasure (recruited erasure) and the one of the auto\/cross-catalysis of chromatin modifications. We collect the definitions and interpretations of these parameters in Fig. 1<\/a>e.<\/p>\n Finally, for modeling the gene expression process, we introduce ({bar{p}}_{O}={alpha }_{O}\/{gamma }_{O})<\/span>, ({bar{p}}_{T}={alpha }_{T}\/{gamma }_{T})<\/span>, ({bar{p}}_{J}={alpha }_{J}\/{gamma }_{J})<\/span>, ({bar{u}}_{O}=({u}_{O}Omega )\/{gamma }_{O})<\/span>, ({bar{u}}_{T}=({u}_{T}Omega )\/{gamma }_{T})<\/span>, ({bar{u}}_{J}=({u}_{J}Omega )\/{gamma }_{J})<\/span> and ({n}_{X}^{A})<\/span>, with X<\/i>\u2009=\u2009O<\/i>, T<\/i>, J<\/i>, that is, the total amount of nucleosomes modified with activating chromatin modifications for each gene X<\/i>.<\/p>\n We next describe the relationship between the abundance of MBD proteins (B) and the rate coefficients associated with the erasure of DNA methylation by the action of TET1. Specifically, these coefficients are ({tilde{k}}_{T}^{{prime} * })<\/span> for the reactions in which TET1 is recruited by DA<\/sup> (reactions \u246d) and ({hat{k}}_{T}^{{prime} })<\/span> for the reactions in which TET1 is not recruited by DA<\/sup> (reactions \u2469). As derived in44<\/a><\/sup>, ({hat{k}}_{T}^{{prime} })<\/span> and ({tilde{k}}_{T}^{{prime} * })<\/span> can be written as:<\/p>\n $${hat{k}}_{T}^{{prime} }=frac{{F}_{1}}{{F}_{2}frac{B}{{K}_{B}}+1},,,{tilde{k}}_{T}^{{prime} * }=frac{{F}_{3}}{{F}_{4}frac{B}{{K}_{B}}+1},$$<\/span><\/p>\n \n (3)\n <\/p>\n<\/div>\n in which K<\/i>B<\/i><\/sub> is the dissociation constant of the binding reaction between B and methylated DNA and F<\/i>1<\/sub>, F<\/i>2<\/sub>, F<\/i>3<\/sub> and F<\/i>4<\/sub> are parameters independent of B or K<\/i>B<\/i><\/sub>. Then, ({hat{k}}_{T}^{{prime} })<\/span> and ({tilde{k}}_{T}^{{prime} * })<\/span> increase when B<\/i>\/K<\/i>B<\/i><\/sub> decreases. From (2<\/a>), ({tilde{mu }}^{{prime} })<\/span> is therefore an increasing function of ({tilde{k}}_{T}^{{prime} * })<\/span> that can be written as:<\/p>\n $${tilde{mu }}^{{prime} }=frac{{tilde{k}}_{T}^{{prime} * }{D}_{tot}}{{k}_{E}^{A}}=frac{{D}_{tot}}{{k}_{E}^{A}}frac{{F}_{3}}{{F}_{4}frac{B}{{K}_{B}}+1}=frac{{F}_{5}}{{F}_{2}frac{B}{{K}_{B}}+1},$$<\/span><\/p>\n \n (4)\n <\/p>\n<\/div>\n in which ({F}_{5}={F}_{3}{D}_{tot}\/{k}_{E}^{A})<\/span> is independent of B or K<\/i>B<\/i><\/sub>. From (4<\/a>), we can conclude that knocking down MBD proteins (B<\/i>\u2009=\u20090) or locally preventing their binding to methylated DNA (K<\/i>B<\/i><\/sub>\u2009\u2192\u2009\u221e<\/i>) allows to increase ({tilde{mu }}^{{prime} })<\/span>.<\/p>\n It has been experimentally observed that, during cellular differentiation, the OCT4 gene undergoes progressive silencing, with concurrent establishment of DNA methylation50<\/a><\/sup>. Here, we investigate how the rate constant of the DNA methylation erasure process by TET1 ({tilde{mu }}^{{prime} })<\/span> affects the first time that, without any external stimulus, the OCT4 gene reaches the repressed state ({n}_{O}^{A}approx 0)<\/span> corresponding to a differentiated state36<\/a><\/sup>, starting from the active state ({n}_{O}^{A}\/{{{{rm{D}}}}}_{{{{rm{tot}}}}}approx 1)<\/span> (Fig. 3<\/a>). For values of ({tilde{mu }}^{{prime} })<\/span> sufficiently large, none of the trajectories reach the OCT4 repressed state in the time window of the simulation, indicating stability of the high OCT4 expression state. By reducing ({tilde{mu }}^{{prime} })<\/span>, most of the trajectories reach, although in a stochastic manner, the OCT4 repressed state. Finally, if we keep reducing ({tilde{mu }}^{{prime} })<\/span>, then the trajectories reach ({n}_{O}^{A}approx 0)<\/span> quickly and in a more synchronous fashion. This result suggests that, even if the initial state, where OCT4 gene is active, is devoid of DNA methylation, the stability of this active state is contingent on a sufficiently fast erasure of DNA methylation.<\/p>\nModel of the epigenetic OCT4 gene regulatory network<\/h3>\n
<\/a><\/div>\n
<\/a><\/div>\n
Effect of DNA methylation on the dynamics of OCT4 during differentiation<\/h3>\n