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

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.

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][4][5]. Yet, analysis of a largescale 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 the 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][16][17][18][19]. The third module was the transition from a prematured 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, ultrasensitivity, 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 15 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 CDC42, (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 was represented by transcription and translation processes by a two-molecule loop (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, feedback loops were converted to three-molecule loops by cascade contraction. The positive-feedback loop between PI3K and CDC42 was reconstructed by selecting PI3K, PIP, and CDC42 from the core network (Additional file 1: Figure S4). The negativefeedback loop between PTEN and GSK3B was reconstructed by selecting PTEN and GSK3B to the loop of these three molecules (Additional file 1: Figure S5). The negative-feedback loop between betacatenin and HES1 was reconstructed by connecting HES1 and GSK3B with beta-catenin (Additional file 1: Figure S6). This network translated into a deterministic mathematical 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 basic 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)  conditions. The model could simulate the HES1 and ASCL1 oscillation within the 150-min period reported by Imayoshi et al. [13] in the undifferentiated state as the basic 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. In order to follow the physiological condition at the initiation of neural differentiation, we set the concentration of NOTCH decreased during simulation. As a result, ASCL1 became dominant in a non-oscillatory state. At the same time, the concentration of GSK3B (a negative regulator [20]) decreased, and the concentration of CDC42 and PI3K (positive regulators) increased (Fig. 4B). These results corresponded to the previously confirmed experimental results during a neural differentiation. Conversely, we set the high concentration of NOTCH and execute simulation. Then the oscillation of the concentration of HES1 and ASCL1 disappeared and the equilibrium concentration of HES1 became higher than that of ASCL1 (Fig. 4C). This result corresponds to the previously confirmed initiation of glial differentiation.
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, 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 10.0, 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 0.1, 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 the loss of functional mutation of amino acids that inhibit the selffeedback regulation of HES1, we collapsed the self-feedback loop of HES1 in our model and execute 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 to PTEN, we collapsed the negative feedback loop between GSK3B and PTEN in our model and execute simulation.
When the negative feedback loop was collapsed by changing equations 6 and 7 to equations 6' and 7', 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 positive feedback loop between CDC42 and PI3K in our model and execute simulation. When the positive feedback loop was collapsed by changing equations 8 and 9 to equations 8' and 9', 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 CDC42−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 to HES1, we collapsed the negative feedback loop between beta-catenin and HES1 in our model and execute simulation. When the negative feedback loop was collapsed 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.

Discussion
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, CDC42, and PI3K (Fig. 4B, C) are consistent with previous experimental result [20]. The results of the simulation of Id2 knockdown or overexpression are also consistent with experimental result [14]. Therefore, our toy model could adequately simulate the dynamics not only of HES1 and ASCL1 but also of other molecules. In relation to the characteristic dynamic responses to NOTCH concentration, our model suggests that three loops (HES1 negative selffeedback, positive feedback between CDC42 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 the most important motif because of its greatest contribution to the characteristic dynamics (Fig. 6D).
A relation between beta-catenin and HES1 is also known as tumorgenesis [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 discovered 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. As the increase of the number of false positive regulations, the increase of the number of nonlinearity cannot be avoidable. Based on this mathematical background, the model with high number of false positive regulations seemingly generates the real behavior, but is apart from the real system; therefore, it is desirable to build a mathematical model from purely 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 core factors of differentiation that also inhibit reprogramming [43] or control mechanisms generated by neurogenic niches [44]. Some positive feedback mechanisms or regulation of multiple cells might underlie this irreversibility, and a more detailed, larger network needs to be analyzed. Our method can be used to analyze the dynamics of a new large-scale network when new information becomes available.

Conclusion
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][16][17][18][19][20][21][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 colorcoded 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 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 to three-molecule loops. After this contraction, two molecules per edge were assigned arbitrarily from among the molecules 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 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 10,000 min 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-