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.