Neural differentiation dynamics controlled by multiple feedback loop motifs in a comprehensive molecular interaction network

doi:10.21203/rs.2.9709/v2

Download PDF

Research article

Neural differentiation dynamics controlled by multiple feedback loop motifs in a comprehensive molecular interaction network

https://doi.org/10.21203/rs.2.9709/v2

This work is licensed under a CC BY 4.0 License

Version 2

posted

You are reading this latest preprint version

Background: Computational simulation using mathematical models is a useful method for understanding the complex behavior of a living system. The majority of studies using mathematical models to reveal biological mechanisms use one of the two main approaches: the bottom-up or the top-down approach. When we aim to analyze a large-scale network, such as a comprehensive knowledge-integrated model of a target phenomenon, for example a whole-cell model, the variation of analyses is limited to particular kind of analysis because of the size and complexity of the model. Results: To analyze a large-scale network of neural differentiation, we developed a hybrid method that combines both approaches. To construct a mathematical model, we extracted network motifs, subgraph structures that recur more often in a metabolic network or gene regulatory network than in a random network, from a large-scale network, detected regulatory motifs among them, and combined these motifs. We confirmed that the model reproduced the known dynamics of HES1 and ASCL1 before and after differentiation, including oscillation and equilibrium of their concentrations. The model also reproduced the effects of overexpression and knockdown of the Id2 gene. Our model suggests that the characteristic change in HES1 and ASCL1 expression in the large-scale network is controlled by a combination of four feedback loops, including a novel large loop discovered in this study. Conclusion: The model extracted by our hybrid method has the potential to reveal the critical mechanisms of neural differentiation. The hybrid method is applicable to other biological events.

Bioinformatics

neural differentiation

regulatory motif

feedback regulation

signaling pathway

mathematical models

Computational simulation using mathematical models is useful for understanding complex systems. Mathematical models of biological mechanisms use either the bottom-up or top-down approach [1]. Both approaches have their own challenges. The bottom-up approach has problems in acquiring regulatory relationships, whereas the top-down approach has problems in parameter determination. In the bottom-up approach, a few molecules of interest may be the pivot of an extended molecular network underlying the target biological event. For example, a mathematical model of the budding yeast cell cycle has been constructed using the bottom-up approach by combining known biochemical reactions [2]. This approach could yield a widely applicable model to simulate each specific analysis subject, but it is difficult to apply to the analysis of comprehensive biochemical networks because of limited information on the regulatory interactions of molecules and the parameters of mathematical models. The top-down approach is used to determine the intrinsic control mechanisms of target biological events. A simple model that simulates characteristic local dynamics is constructed using pathway databases or in a data-driven manner based on omics data by excluding non-essential factors from a comprehensive network. In the absence of a comprehensive pathway map, a statistical model is one of the first choices to construct a comprehensive network [3–5]. Yet, analysis of a large-scale network is still challenging, given that many parameters need to be determined in advance. Estimation of all the parameters of a whole-cell model has not yet been appropriately resolved [6]. To analyze the dynamics of a large-scale network, it is often divided into feasible-size modules, which are defined as small networks of functional units amenable to simulation and analysis [7]. Overall, the bottom-up approach can be used to analyze part of a biological event around a source molecule, but it is not applicable to a large-scale network; the top-down approach cannot be applied to analyze the dynamics because of the problem of parameter estimation. A whole-cell mathematical model of the bacterium Mycoplasma genitalium, which contains 525 genes, was built on the basis of enormous experimental data [8]; however, the use of a life-cycle model is not a simple way to analyze the mechanisms of dynamic control and requires a lot of information about the particular species.

To comprehensively analyze the dynamics of regulatory mechanisms, we developed a hybrid method that combines the bottom-up and top-down approaches. We aimed to decrease the number of modeled elements without losing the characteristic dynamics; therefore, we focused on network motifs that are important for the dynamics. Network motifs are subgraph structures that recur more often in a metabolic network or gene regulation network than in a random network. They are important for determining intrinsic regulation mechanisms derived from network characteristics [9]. In the current method, we extracted motifs from a large-scale network and used them to reconstruct a simple network, which reflects the original dynamics of the entire network because it contains these important motifs.

We then modeled the differentiation of neural stem cells (NSCs) and analyzed the characteristic changes of dynamics during differentiation. NSCs replicate and differentiate into neurons, astrocytes, or oligodendrocytes [10]. Some models simulate early differentiation or functional neurons [11, 12], but no model enables us to simulate and analyze the dynamics of the large-scale network of neuronal differentiation. The basic helix-loop-helix type transcription factors HES1, ASCL1, and OLIG2 show characteristic differences in their dynamics before and after differentiation [13]. They also maintain oscillatory dynamics during cell replication. If the concentration of ASCL1 is higher than that of HES1 during the non-oscillatory state, the NSCs differentiate into neurons. If the concentration of HES1 or OLIG2 is higher than that of ASCL1, the NSCs differentiate into astrocytes or oligodendrocytes, respectively [13]. In the current study, we constructed a comprehensive molecular-interaction network of NSC differentiation into neurons using the available data; we then developed a mathematical model that maintained the original dynamics of the network by integrating network motifs. The model could simulate the characteristic dynamic changes before and after differentiation. The model also reproduced the effect of overexpression or knockdown of the Id2 gene, which encodes an inhibitor of HES1 dimerization [14]. We suggest that stabilization of oscillations and characteristic dynamic changes are regulated by the combination of multiple feedback loops.

Our method allows the analysis of a comprehensive network by collecting information exhaustively and scaling down the network according to the rationality of dynamics without arbitrariness. On the basis of the analysis of the dynamics of neuronal differentiation, we here propose that a combination of multiple motifs is important to define the major dynamics of an entire network.

Signaling network of neuronal differentiation

To construct a signaling network of NSC differentiation into a mature neuron, we collected publicly available information about the switch from NOTCH, a molecular marker of differentiation, to neuronal markers such as Tau. The network was constructed by using 54 molecules and contained five modules (Fig. 1A). The first module was the differentiation switch from NOTCH. The second module was the expression of transcription factors that are early neural markers [15–19]. The third module was the transition from an immature neuron to a mature neuron. The fourth module was regulation to gain mature neuron functions [20, 21]. Finally, beta-catenin, a molecule related to the function of mature neurons, controls a bHLH-type transcription factor to adjust differentiation [22]. Converting this entire signaling network into a mathematical model can be challenging because multiple parameters need to be taken into account. Because the network dynamics are controlled by the dynamics of individual network motifs [23], we focused on feedback loop motifs that are important to nonlinear dynamics. For example, a feedback loop may confer the ability of homeostasis, ultra-sensitivity, hysteresis, and bistability [24]. A positive feedback loop is defined as a loop structure containing zero or an even number of negative regulations, and a negative feedback loop is defined as a loop structure containing an odd number of negative regulations (Fig. 2A). Because it was difficult to extract feedback loop motifs from a large-scale network, the entire network was contracted to a loop-only structure (Fig. 2B, Additional file 1: Figure S1 and Table S1). After contraction, feedback loop structures were extracted to obtain the core network, which contained 14 molecules (Fig. 1B). This network included four feedback loop motifs: (1) a negative-feedback loop formed by HES1 dimerization, (2) a positive-feedback loop between PI3K and aPKC_PAR3_PAR6, (3) a negative-feedback loop between PTEN and GSK3B, and (4) a negative-feedback loop between beta-catenin and HES1. The first three loops have been previously investigated [16, 20] (the first one has been especially well analyzed [25, 26]), but the negative-feedback loop between beta-catenin and HES1, which was the largest in our model, has not been reported. The core network retained the feedback loop motifs of the original entire network, and therefore was expected to maintain the original dynamics. To analyze the dynamics of the core network, we constructed a mathematical toy model.

Mathematical model of the core network

Because the core network was still too complex to convert into a mathematical model, we contracted it to the minimum nodes but maintained the feedback loop structures, including the differentiation switch NOTCH and the indicators of characteristic dynamics before and after differentiation (Fig. 3, Additional file 3). The self-negative-feedback loop formed by HES1 dimerization was reconstructed as a three-molecule loop by adding mRNA of HES1 (mHES1) as a transcriptional process, because the self-negative feedback loop represented by transcriptional and translational processes was actually composed of the generation of mRNA and translation of the mRNA (see Methods for details). Although previous studies introduced a delay into the model [11, 12], we convolved the delay with the parameters of the HES1 translation and dimerization steps. The non-delay model of the HES1 self-loop generated oscillations (Additional file 1: Figures S2, S3, and Table S2). To minimize the network, most feedback loops were converted to three-molecule loops by cascade contraction. The contracted cascades were represented with the nodes in our toy model. The names of the nodes were determined based on the name of the first or the last node of the original cascade with the addition of the suffix _ca or _ci. The suffix _ca denoted the active form of a cascade, and _ci denoted the inactive form. The positive-feedback loop between PI3K and aPKC_PAR3_PAR6 was reconstructed as a three-molecule loop, and the contracted nodes were named PIP_ca (PIP_ci), aPKC_ca (aPKC_ci), and PI3K_ca (PI3K_ci) (Additional file 1: Figure S4). The negative-feedback loop between PTEN and GSK3B was reconstructed as a four-molecule loop, and the contracted nodes were named PTEN_ca (PTEN_ci) and GSK3B_ca (GSK3B_ci) (Additional file 1: Figure S5). The negative-feedback loop between beta-catenin and HES1 was reconstructed as a six-molecule loop (Additional file 1: Figure S6). This network translated into a deterministic mathematical toy model governed by the Hill equation and Michaelis–Menten kinetics (Table 1). The final toy model was constructed using 9 molecules (16 nodes) and 20 reactions, which reduced the network size by 84.3% in comparison with the entire signaling network. Thereafter, the dynamics of neuronal differentiation were analyzed by using the toy model.

Simulation of the oscillatory dynamics

To simulate the oscillatory state of HES1 and ASCL1 in undifferentiated NSCs as a basal condition (as reported by Imayoshi et al. [13]), the toy model was investigated using an oscillatory parameter set. The ranges of 40 parameters (Table 2) allowed the model to maintain oscillations. To check whether these parameter values were physiologically relevant, we compared them to the general parameter values of Michaelis–Menten kinetics based on the BioNumbers database [27]. The minimal values of the parameters in Hill equation that allowed to reconstruct the expected behaviors of the model were chosen. Almost all parameters of enzymes, except kM_re2_s24 and kM_re2_s9, exhibit K_Mvalues above 0.1 µM [28], and almost all of the K_M values in the model were above 0.1 µM in the oscillatory range. The maximum (37,596,000 h⁻¹) and minimum (264 h⁻¹) K_cat values were acquired from 27 records for mammals in BioNumbers (Additional File 2); all K_cat values in the model were within this range. These data indicated that the oscillatory state could be established under physiologically relevant conditions. The model could simulate the HES1 and ASCL1 oscillation within the 2.5-h period reported by Imayoshi et al. [13] in the undifferentiated state as the basal condition (Fig. 4A, Additional file 1: Figure S7).

Model validation

During neural differentiation, the characteristic dynamics of the transition from the oscillatory to the non-oscillatory state of HES1 and ASCL1 are controlled by the concentration of NOTCH [13]. In the non-oscillatory state after differentiation, the concentration of ASCL1 was higher than that of HES1 in a neuron, whereas the concentration of HES1 was higher than that of ASCL1 in an astrocyte. To simulate the physiological condition at the initiation of neuronal differentiation, we set a low concentration of NOTCH during simulation. As a result, ASCL1 became dominant in a non-oscillatory state. At the same time, the concentration of GSK3B_ca, which included GSK3B (a negative regulator [20]), decreased and the concentrations of aPKC_ca, which included aPKC_PAR3_PAR6, and PI3K_ca, which included PI3K (positive regulators [20]), increased (Fig. 4B). These results agreed with the previously reported experimental results during neuronal differentiation [13, 20]. Conversely, we set a high concentration of NOTCH and executed simulation. The oscillation of the concentrations of HES1 and ASCL1 disappeared and the equilibrium concentration of HES1 became higher than that of ASCL1 (Fig. 4C). This result agrees with the previously reported initiation of glial differentiation [13]. Both ASCL1 and HES1 maintained oscillations at physiologically relevant NOTCH concentrations in an NSC (Fig. 4D) [29]. This NOTCH-dependent dynamic transition was consistent with the experimental results of Imayoshi et al. [13].

We also simulated the overexpression and knockdown of the Id2 gene, which encodes an inhibitor of HES1 dimerization [14]. The knockdown of Id2 promotes neuronal differentiation by suppressing HES1 and enhancing ASCL1 expression, and the overexpression of Id2 inhibits neuronal differentiation by enhancing HES1 and suppressing ASCL1 expression [14]. To simulate the inhibition of HES1 dimerization by Id2, an inhibition parameter, kSm_Id (an indicator of the amount of Id2), was introduced into equations 1 and 3 (Table 1) to obtain equations 1′ and 3′ (Table 3), respectively. When Id2 knockdown was simulated by setting kSm_Id to 0.1, oscillations disappeared, and a neuronal differentiation state with ASCL1 domination was observed at each NOTCH concentration (Fig. 5A). When Id2 overexpression was simulated by setting kSm_Id to 10, oscillations also disappeared, and a non-neuronal differentiation state with HES1 domination was observed at NOTCH concentrations over 0.7 µM (Fig. 5B). Thus, the toy model maintained the original dynamics of the entire network.

Network motif analysis

We analyzed the contribution of each motif to the results of simulation of NOTCH responsiveness by bifurcation analysis. To simulate a loss-of-function mutation that inhibits the self-feedback regulation of HES1, we collapsed the self-feedback loop of HES1 in our model and executed simulation. When the loop was collapsed by changing equation 2 to equation 2’ (Table 4), the oscillatory state, which is important for maintaining the undifferentiated state, disappeared (Fig. 6A). This result was consistent with the reported importance of the HES1 self-feedback loop in neuronal differentiation [25, 26]. To simulate the specific inhibition of the effect of GSK3B on PTEN, we collapsed the negative feedback loop between GSK3B and PTEN by changing equations 5 and 6 to equations 5’ and 6’, respectively (Table 4). The oscillatory state shifted to higher NOTCH concentrations (Fig. 6B), which means that this loop was involved in the sensitivity of oscillations to NOTCH concentration. GSK3B may promote or inhibit NOTCH signaling under different conditions depending on the loops involved [30, 31, 32, 33]. The dichotomic characteristics of GSK3B may enable it to arbitrate the response to NOTCH, given that our simulation suggested that the GSK3B−PTEN loop regulated the sensitivity of oscillations to NOTCH. To simulate the specific inhibition of the effect of PI3K to PIP2, we collapsed the positive feedback loop between aPKC_PAR3_PAR6 and PI3K by changing equations 7 and 8 to equations 7’ and 8’, respectively (Table 4). The relative changes in HES1 and ASCL1 concentrations with NOTCH concentration were maintained, but the oscillatory state disappeared (Fig. 6C); therefore, the aPKC_PAR3_PAR6−PI3K loop is required to maintain the oscillatory state. We estimated the effect of the collapse of this loop by inhibiting PI3K. A PI3K inhibitor, LY294002, inhibits proliferation of neural progenitor cells [22, 33]; our result agreed with the published data. To simulate the specific inhibition of the effect of beta-catenin on HES1, we collapsed the negative feedback loop between beta-catenin and HES1 by changing equation 2 to equation 2’' (Table 4). The oscillations disappeared, and the ASCL1-dominant state became narrower (Fig. 6D), indicating that this loop is required for maintaining oscillations and upregulating ASCL1 at low concentrations of NOTCH. The collapse of this loop equals beta-catenin inhibition, which represses proliferation of neural progenitor cells and accelerates glial differentiation [34, 35]. Glial differentiation is induced without oscillation under ASCL1 repression [13]. Our result that the collapse of this loop leads to ASCL1 repression is consistent with the reported acceleration of glial differentiation by beta-catenin inhibition. These results suggest that multiple feedback loops are essential for the characteristic dynamics of neural differentiation.

We generated an NSC differentiation network containing four feedback loops on the basis of publicly available data. Our toy model was constructed through cascade contraction of the comprehensive network with preservation of feedback loops. Three types of HES1 and ASCL1 states regulated by NOTCH concentration were consistent with NOTCH-dependent neural differentiation suggested by Imayoshi et al. [13]. Although experimental data show complex waveforms of HES1 and ASCL1, we simulated the main-frequency waves, which have a period of 2 to 3 h [13], and analyzed the dynamics of the transition from the oscillatory to the non-oscillatory state qualitatively using a toy model. The results of the simulation of GSK3B, aPKC_PAR3_PAR6, and PI3K, represented as GSK3B_ca, aPKC_ca, and PI3K_ca in the toy model respectively, are consistent with previous experimental results [20] (Fig. 4B, C). The results of the simulation of Id2 knockdown or overexpression are also consistent with experimental results [14]. Therefore, our toy model could adequately simulate the dynamics not only of HES1 and ASCL1 but also of other molecules. Our model suggests that three loops (HES1 negative self-feedback, positive feedback between aPKC_PAR3_PAR6 and PI3K, and negative feedback between GSK3B and HES1) are important to maintain undifferentiated state oscillations. We suggest that the newly discovered negative-feedback loop between beta-catenin and HES1 in the comprehensive network is most important because of its greatest contribution to the characteristic dynamics (Fig. 6D). A relation between beta-catenin and HES1 plays a role in tumorigenesis [36]. As HES1 controls cancer stem cells [37], the negative feedback loop that we newly discovered may be related to proliferation and differentiation of cancer stem cells. These findings could only be made by using the analysis based on a large-scale network, thus highlighting the effectiveness of our approach.

We demonstrated that focusing on feedback loop motifs instead of the whole network when constructing a model was sufficient for agreement with experimental results. Our approach could be applied to analysis of various biochemical networks by simulation. By streamlining large-scale network construction, our approach could help to analyze various biological phenomena, such as cell differentiation, cell division, or pathogenesis. However, the large-scale network will be probably insufficient and heterogeneous when it is constructed using the available data alone. To overcome this false-negative problem (the relations that exist but cannot be detected), many data-driven network reconstruction methods have been developed. These statistical approaches are mainly classified under two categories, expression based [38] and sequence based [39]. Although the methods of both categories could reveal undiscovered relations that could not be inferred manually, the reconstructed network includes many false-positive regulations. A nonlinear model intrinsically causes a complex behavior. With an increase in the number of false-positive regulations, an increase in the number of nonlinearities becomes avoidable. Based on this mathematical background, a model with a high number of false-positive regulations seemingly generates the real behavior, but is different from the real system; therefore, it is desirable to build a mathematical model only from reliable elements. Addition of false-positive regulation to the model could have a considerable effect and complicate the conversion of a data-driven network into a mathematical model. Currently, the manual methods are better than the data-driven methods for construction of a mathematical model.

Our large-scale network of neuronal differentiation may lack some components and thus may not completely represent neuronal differentiation. In our network, ASCL1 is directly affected by HES1; ASCL1 oscillation controls proliferation and differentiation [13] and affects NOTCH receptors of adjacent cells via activation of DLL [40]. Neural differentiation is also affected by adjacent cells [41]. We excluded this information because we focused on the dynamics of a single cell. To reveal the entire mechanism of neural differentiation, adding a path to adjacent cells, for instance via DLL, might be required. The analysis of multiple cells might provide a model that can simulate dynamics other than state transitions. The concentrations of both HES1 and ASCL1 decrease in a non-oscillatory state when an NSC differentiates into an oligodendrocyte [13]. To simulate this transition, we need to add a signaling pathway focusing on the oligodendrocyte marker OLIG2, which oscillates with a period of 400 min in NSCs [13]. This period is much longer than that of HES1 or ASCL1, and the dispersion of the oscillation is very high; therefore, OLIG2 regulation might involve a delay mechanism to elongate the period and a mechanism to amplify dispersion. Recently, a similar method was used to analyze oligodendrocyte differentiation [42]. Similar to our study, the authors used a manual method to construct a network; however, they also introduced publicly available interaction data from omics databases. In comparison with our method, this approach may reduce the number of false-negative interactions. On the other hand, the study [42] focused only on two- to four-node feedback loops. Our contraction method may detect a larger regulatory system of oligodendrocyte differentiation. The complete mechanism of neural differentiation may be simulated by integrating this information and methods.

Our model can simulate the dynamics of the NSC-to-neuron transition and exemplify the reverse transition by increasing the concentration of NOTCH, but differentiation is mostly irreversible. Therefore, it is difficult to validate the results of reversing from a neuron to an NSC. Some hypotheses suggest the core factors of differentiation that also inhibit reprogramming or control the mechanisms generated by neurogenic niches. Specific network structure such as positive feedbacks or micro-environmental factors may be important for hysteresis in differentiation, and a more detailed, larger network needs to be analyzed. Although we focused on NOTCH signaling in this study, our network also includes FGF as another input signal. Analysis of the behavior of the network stimulated with FGF may show variate responses and as a result may reveal other mechanisms of neural differentiation. Our method can be used to analyze the dynamics of a new large-scale network when new information becomes available.

The construction of a large-scale network of neuronal differentiation based on publicly available data led to identification of a new feedback loop motif, which we expect to regulate differentiation. The large-scale network was modeled mathematically after network contraction, and the extracted toy model simulated the characteristic dynamics of HES1 and ASCL1, which were suggested to be regulated by multiple loops. More information about neural differentiation and further motif analyses will deepen our understanding of the mechanism of neural differentiation. Our approach is applicable to other biological models for which detailed mechanisms are unknown.

Construction of a neuronal differentiation network

To construct the complete signaling network, relations among molecules involved in neuronal differentiation were collected from the literature and pathway databases [15–22], mainly the NCI Pathway Interaction Database [45] and WikiPathways [46]. Molecular interactions related to glial differentiation were excluded to focus on the differentiation of NSCs into neurons. The signaling networks were constructed using CellDesigner v4.3 [47] by merging binary relations from each data source and were saved in Systems Biology Markup Language (SBML) format. The complete signaling network was constructed by manually merging all the SBML files. Then, every molecule was color-coded according to its type (receptor, enzyme, transcription factor, or other). Finally, the enzymes and other molecules were divided into active and inactive forms. The graphical representation conformed to the proposed set of symbols in CellDesigner [48].

Contraction of the network

To simplify the network without losing the essential loop structures that are important to the original dynamics, the sequences of unidirectional edges (cascade motifs) were converted to a single edge between two molecules. The rate of a cascade reaction depends on the rate-limiting reaction; therefore, each cascade was represented as a single rate-limiting reaction. The contraction was continued until a network contained hub nodes only (or nearly so). After cascade contraction, feedback loop structures were extracted as a core network. At this time, a transcription–translation self-feedback loop, which was initially constructed from one or two molecules, was reconstructed to include three molecules by adding a transcription or translation event, because such events require much more time than signal transduction events, and this investment of time is related to nonlinear dynamics. Finally, to minimize the network size, feedback loops in the core network were contracted again until at least three nodes. After this contraction, the nodes of contracted cascade were named on the basis of the node name of the original cascade.

Mathematical model construction

The contracted network was converted into a mathematical model (toy model), in which the dynamics of the original network were maintained. The model was constructed using CellDesigner. The contracted network was drawn and kinetics with parameters were assigned to every edge using the SBMLsqueezer [49] plug-in in CellDesigner. In SBMLsqueezer, the ‘Reversibility’ option was set to ‘Use information from SBML’. Every type of enzyme kinetics was set to the Michaelis–Menten equation, which is one of the best-known models of enzyme kinetics. That of transcriptional or translational reactions was set to the Hill equation, which is the simplest way to describe sigmoidal responses. The parameters were estimated in the range of physiologically relevant values which could simulate the experimental results reported by Imayoshi et al. [13]. The range was determined from records in BioNumbers [27, 28]. The initial concentration of each protein was set to 1.0 µM because the typical concentration range of a signaling protein is 10 nM–1 µM [29]. For enzymes and other molecules, the concentration of the active and inactive form was set to 0.5 µM each.

Simulation and analysis

The dynamics of the toy model were validated by steady-state simulation and parameter analysis. Simulation was performed using CellDesigner, and the analysis was conducted using COPASI 4.14 (Build 89) [50]. The simulation was calculated using SOSlib [51] with the error tolerance set to −6. To simulate the change in dynamics due to the induction of differentiation, an event that perturbs the concentration of NOTCH (differentiation control factor) at an arbitrary time point was configured. The parameter search and bifurcation analysis were conducted using the Parameter Scan function in COPASI, which records the time course of an arbitrary molecule for 500-h while changing an arbitrary parameter or concentration within the defined range. The time course was calculated using LSODA [52] with the following parameters: Integrate Reduced Model, 0; Relative Tolerance, 1e–06; Absolute Tolerance, 1e–12; Max Internal Steps, 10,000. The scope of a parameter scan was set to 0.001–1000, and the range of NOTCH concentrations as a differentiation control factor was set to 0–2.8µM.

aPKC: atypical protein kinase C

aPKC_PAR3_PAR6: aPKC in complex with PAR3 and PAR6

ARP2/3: actin-related protein 2/3
ASCL1: achaete-scute family bHLH transcription factor 1

bHLH: basic-helix-loop-helix
CDC42: cell division control protein 42 homolog
DLL: delta-like canonical Notch ligand

GSK3B: glycogen synthase kinase 3 beta
HES1: hairy and enhancer of split 1 

Id2: inhibitor of DNA-binding 2 
NSCs: neural stem cells

OLIG2: oligodendrocyte transcription factor 2

PAR3: partitioning defect 3

PAR6: partitioning defect 6
PI3K: phosphatidylinositol-3 kinase
PIP2: phosphatidylinositol 4,5-bisphosphate

PIP3: phosphatidylinositol 4,5-trisphosphate
PTEN: phosphatase and tensin homolog

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and material

Simulation and analysis codes of CellDesigner and COPASI are available upon request.

Competing interests

The authors declare that they have no conflict of interest.

Funding

Sumitomo Dainippon Pharma Co., Ltd. provided support in the form of salaries for T. Iwasaki and K. Yamazaki and as part of the research budget to N. Hiroi. The company did not play any role in the study design, data collection, or analysis. The specific roles of these authors are described in the Authors’ Contributions section. This work was supported by the Japan Society for the Promotion of Science (JSPS), a Grant-in-Aid for Scientific Research (Kiban C) 15K06925 awarded to N. Hiroi.

Authors’ contributions

Conceived and designed the experiments: TI, TY, NH, AF, KY. Performed the analyses: TI, RT, TY, TH. Wrote the paper: TI, TY, NH, AF. Critical advice on the manuscript: TH, KY. All authors have read and approved the manuscript.

Acknowledgements

We thank Dr. Tetsuya J. Kobayashi (Institute of Industrial Science, the University of Tokyo) for his constructive comments and encouragement. We thank ELSS (Tsukuba, Japan) for their critical advice on English usage and proofreading of the manuscript.

Bruggeman FJ, Westerhoff HV. The nature of systems biology. Trends Microbiol. 2007; 15(1):45–50.
Chen KC, Calzone L, Csikasz-Nagy A, Cross FR, Novak B, Tyson JJ. Integrative analysis of cell cycle control in budding yeast. Mol Biol Cell. 2004;15(8):3841–62.
Vignes M, Vandel J, Allouche D, Ramadan-Alban N, Cierco-Ayrolles C, Schiex T, et al. Gene regulatory network reconstruction using Bayesian networks, the Dantzig Selector, the Lasso and their meta-analysis. PLoS One. 2011;6(12):e29165.
Chai LE, Loh SK, Low ST, Mohamad MS, Deris S, Zakaria Z. A review on the computational approaches for gene regulatory network construction. Comput Biol Med. 2014;48:55–65.
Park Y, Kellis M. Deep learning for regulatory genomics. Nat Biotechnol. 2015;33(8):825–6.
Karr JR, Williams AH, Zucker JD, Raue A, Steiert B, Timmer J, et al. Summary of the DREAM8 Parameter Estimation Challenge: Toward Parameter Identification for Whole-Cell Models. PLoS Comput Biol. 2015;11(5):e1004096.
Yugi K, Kubota H, Toyoshima Y, Noguchi R, Kawata K, Komori Y et al. Reconstruction of insulin signal flow from phosphoproteome and metabolome data. Cell Rep. 2014;8(4):1171–83.
Karr JR, Sanghvi JC, Macklin DN, Gutschow MV, Jacobs JM, Bolival B Jr, et al. A whole-cell computational model predicts phenotype from genotype. Cell. 2012;150(2):389–401.
Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U. Network motifs: simple building blocks of complex networks. Science. 2002;298(5594):824–7.
Louvi A, Artavanis-Tsakonas S. Notch signalling in vertebrate neural development. Nat Rev Neurosci. 2006;7(2):93–102.
Monk NA. Oscillatory expression of Hes1, p53, and NF-kappaB driven by transcriptional time delays. Curr Biol. 2003;13(16):1409–13.
Zeiser S1, Müller J, Liebscher V. Modeling the Hes1 oscillator. J Comput Biol. 2007;14(7):984–1000.
Imayoshi I, Isomura A, Harima Y, Kawaguchi K, Kori H, Miyachi H, et al. Oscillatory control of factors determining multipotency and fate in mouse neural progenitors. Science. 2013;342(6163):1203–8.
Bai G, Sheng N, Xie Z, Bian W, Yokota Y, Benezra R, Kageyama R, Guillemot F, Jing N. Id sustains Hes1 expression to inhibit precocious neurogenesis by releasing negative autoregulation of Hes1. Dev Cell. 2007;13(2):283–97.
Kageyama R, Ohtsuka T, Hatakeyama J, Ohsawa R. Roles of bHLH genes in neural stem cell differentiation. Exp Cell Res. 2005;306(2):343–8.
Kageyama R, Ohtsuka T, Kobayashi T. Roles of Hes genes in neural development. Dev Growth Differ. 2008;50 Suppl 1:S97–103.
Seki, T., Sawamoto, K., Parent, J. M., Alvarez-Buylla, A, editors. Neurogenesis in the Adult Brain, vol. 1. Tokyo: Springer; 2011.
Roybon L, Mastracci TL, Ribeiro D, Sussel L, Brundin P, Li JY. GABAergic differentiation induced by Mash1 is compromised by the bHLH proteins Neurogenin2, NeuroD1, and NeuroD2. Cereb Cortex. 2010;20(5):1234–44.
Bhat KM, Maddodi N, Shashikant C, Setaluri V. Transcriptional regulation of human MAP2 gene in melanoma: role of neuronal bHLH factors and Notch1 signaling. Nucleic Acids Res. 2006;34(13):3819–32.
Arimura N, Kaibuchi K. Neuronal polarity: from extracellular signals to intracellular mechanisms. Nat Rev Neurosci. 2007;8(3):194–205.
Hand R, Bortone D, Mattar P, Nguyen L, Heng JI, Guerrier S, et al. Phosphorylation of Neurogenin2 specifies the migration properties and the dendritic morphology of pyramidal neurons in the neocortex. Neuron. 2005;48(1):45–62.
Shimizu T, Kagawa T, Inoue T, Nonaka A, Takada S, Aburatani H, et al. Stabilized beta-catenin functions through TCF/LEF proteins and the Notch/RBP-Jkappa complex to promote proliferation and suppress differentiation of neural precursor cells. Mol Cell Biol. 2008;28(24):7427–41.
Riccione KA, Smith RP, Lee AJ, You L. A synthetic biology approach to understanding cellular information processing. ACS Synth Biol. 2012;1(9):389–402.
Eberhardt M, Lai X, Tomar N, Gupta S, Schmeck B, Steinkasserer A, Schuler G, Vera J. Third-kind encounters in biomedicine: immunology meets mathematics and informatics to become quantitative and predictive. Methods Mol Biol. 2016;1386:135–79.
Hirata H, Yoshiura S, Ohtsuka T, Bessho Y, Harada T, Yoshikawa K, Kageyama R. Oscillatory expression of the bHLH factor Hes1 regulated by a negative feedback loop. Science. 2002;298(5594):840–3.
Kageyama R, Ohtsuka T, Kobayashi T. The Hes gene family: repressors and oscillators that orchestrate embryogenesis. Development. 2007;134(7):1243–51.
Milo R, Jorgensen P, Moran U, Weber G, Springer M. BioNumbers—the database of key numbers in molecular and cell biology. Nucleic Acids Res. 2010;38:D750–3.
Bar-Even A, Noor E, Savir Y, Liebermeister W, Davidi D, Tawfik DS, Milo R. The moderately efficient enzyme: evolutionary and physicochemical trends shaping enzyme parameters. Biochemistry. 2011;50(21):4402–10.
Legewie S, Herzel H, Westerhoff HV, Blüthgen N. Recurrent design patterns in the feedback regulation of the mammalian signalling network. Mol Syst Biol. 2008;4:190.
Foltz DR, Santiago MC, Berechid BE, Nye JS. Glycogen synthase kinase–3beta modulates notch signaling and stability. Curr Biol. 2002;12(12):1006–11.
Guha S, Cullen JP, Morrow D, Colombo A, Lally C, Walls D, Redmond EM, Cahill PA. Glycogen synthase kinase 3 beta positively regulates Notch signaling in vascular smooth muscle cells: role in cell proliferation and survival. Basic Res Cardiol. 2011;106(5):773–85.
Jin YH, Kim H, Oh M, Ki H, Kim K. Regulation of Notch1/NICD and Hes1 expressions by GSK–3alpha/beta. Mol Cells. 2009;27(1):15–9.
Kim WY, Wang X, Wu Y, Doble BW, Patel S, Woodgett JR, Snider WD. GSK–3 is a master regulator of neural progenitor homeostasis. Nat Neurosci. 2009;12(11):1390–7.
Ye F, Chen Y, Hoang T, Montgomery RL, Zhao XH, Bu H, Hu T, Taketo MM, van Es JH, Clevers H, Hsieh J, Bassel-Duby R, Olson EN, Lu QR. HDAC1 and HDAC2 regulate oligodendrocyte differentiation by disrupting the beta-catenin-TCF interaction. Nat Neurosci. 2009;12(7):829–38.
Zhang C, Zhang Z, Shu H, Liu S, Song Y, Qiu K, Yang H. The modulatory effects of bHLH transcription factors with the Wnt/beta-catenin pathway on differentiation of neural progenitor cells derived from neonatal mouse anterior subventricular zone. Brain Res. 2010;1315:1–10.
Peignon G, Durand A, Cacheux W, Ayrault O, Terris B, Laurent-Puig P, Shroyer NF, Van Seuningen I, Honjo T, Perret C, Romagnolo B. Complex interplay between -catenin signalling and Notch effectors in intestinal tumorigenesis. Gut. 2011;60(2):166–76.
Liu ZH, Dai XM, Du B. Hes1: a key role in stemness, metastasis and multidrug resistance. Cancer Biol Ther. 2015;16(3):353–9.
Chasman D, Fotuhi Siahpirani A, Roy S. Network-based approaches for analysis of complex biological systems. Curr Opin Biotechnol. 2016; 39:157–166.
McLeay RC, Bailey TL. Motif Enrichment Analysis: a unified framework and an evaluation on ChIP data. BMC Bioinformatics. 2010; 11:165.
Morimoto M, Nishinakamura R, Saga Y, Kopan R. Different assemblies of Notch receptors coordinate the distribution of the major bronchial Clara, ciliated and neuroendocrine cells. Development. 2012;139(23):4365–73.
Ramos C, Rocha S, Gaspar C, Henrique D. Two Notch ligands, Dll1 and Jag1, are differently restricted in their range of action to control neurogenesis in the mammalian spinal cord. PLoS One. 2010;5(11):e15515.
Cantone M, Küspert M, Reiprich S, Lai X, Eberhardt M, Göttle P, Beyer F, Azim K4, Küry P, Wegner M, Vera J. A gene regulatory architecture that controls region-independent dynamics of oligodendrocyte differentiation. Glia. 2019;67(5):825–843.
Hikichi T, Matoba R, Ikeda T, Watanabe A, Yamamoto T, Yoshitake S, et al. Transcription factors interfering with dedifferentiation induce cell type-specific transcriptional profiles. Proc Natl Acad Sci U S A. 2013;110(16):6412–7.
Real C, Glavieux-Pardanaud C, Le Douarin NM, Dupin E. Clonally cultured differentiated pigment cells can dedifferentiate and generate multipotent progenitors with self-renewing potential. Dev Biol. 2006;300(2):656–69.
Schaefer CF, Anthony K, Krupa S, Buchoff J, Day M, Hannay T, et al. PID: the Pathway Interaction Database. Nucleic Acids Res. 2009;37:D674–9.
Kelder T, van Iersel MP, Hanspers K, Kutmon M, Conklin BR, Evelo CT, et al. WikiPathways: building research communities on biological pathways. Nucleic Acids Res. 2012; 40:D1301–7.
Funahashi A, Matsuoka Y. Jouraku A, Morohashi M, Kikuchi N, Kitano H. CellDesigner 3.5: A Versatile Modeling Tool for Biochemical Networks. Proceedings of the IEEE. 2008;96(8): 1254–1265
Kitano H, Funahashi A, Matsuoka Y, Oda K. Using process diagrams for the graphical representation of biological networks. Nat Biotechnol. 2005;23(8):961–6.
Dräger A, Hassis N, Supper J, Schröder A, Zell A. SBMLsqueezer: a CellDesigner plug-in to generate kinetic rate equations for biochemical networks. BMC Syst Biol. 2008;2:39.
Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, et al. COPASI—a COmplex PAthway SImulator. Bioinformatics. 2006;22(24):3067–74.
Machné R, Finney A, Müller S, Lu J, Widder S, Flamm C. The SBML ODE Solver Library: a native API for symbolic and fast numerical analysis of reaction networks. Bioinformatics. 2006;22(11):1406–7.
Petzold L. Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM J Sci Stat Comput. 1983; 4:136–48.

Due to technical limitations, all Tables are only available as a download in the supplemental files section.

Additional file 1—Additional_file1.pdf

Intermediate network before core network. Figure S1.

Nodes removed by contraction. Table S1.

Detailed information about the HES1 self-loop model. Figures S2, S3, and Table S2.

Network structures of feedback loops. Figures S4, S5, and S6.

Simulation result for all molecules. Figure S7.

Additional file 2—Additional_file2.xlsx

K_cat values acquired from 65 records in BioNumbers.

Additional file 3—Additional_file3.xml

Mathematical model of the core network.

Download PDF

Version 2

posted

You are reading this latest preprint version

Neural differentiation dynamics controlled by multiple feedback loop motifs in a comprehensive molecular interaction network

Status:

Version 2

Abstract

Figures

Background

Results

Signaling network of neuronal differentiation

Mathematical model of the core network

Simulation of the oscillatory dynamics

Model validation

Network motif analysis

Discussion

Conclusion

Methods

Construction of a neuronal differentiation network

Contraction of the network

Mathematical model construction

Simulation and analysis

Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and material

Competing interests

Funding

Authors’ contributions

Acknowledgements

References

Tables

Additional File Legend

Supplementary Files

Status:

Version 2