Unraveling the differential dynamics of developmental fate in central and peripheral nervous systems

Bone morphogenetic protein 2 (BMP2), differentially regulates the developmental lineage commitment of neural stem cells (NSC’s) in central and peripheral nervous systems. However, the precise mechanism beneath such observations still remains illusive. To decipher the intricacies of this mechanism, we propose a generic mathematical model of BMP2 driven differentiation regulation of NSC’s. The model efficiently captures the dynamics of the wild-type as well as various mutant and over-expression phenotypes for NSC’s in central nervous system. Our model predicts that the differential developmental dynamics of the NSC’s in peripheral nervous system can be reconciled by altering the relative positions of the two mutually interconnected bi-unstable switches inherently present in the steady state dynamics of the crucial developmental fate regulatory proteins as a function of BMP2 dose. This model thus provides a novel mechanistic insight and has the potential to deliver exciting therapeutic strategies for neuronal regeneration from NSC’s of different origin.

regulates Mash1 level and results in neurogenesis. (m) Overexpression of Id1 gene (G idt =5 X WT), (n) Overexpression of Hes5 gene (G 5t =5 X WT) and (o) Overexpression of Hes1 gene (G t =5 X WT) at high BMP2 (BMP2=20 s.u.) result in the down-regulation of Mash1 level and indicate gliogenic fate commitment in all the three cases. In all the figures k bmp2 =100 min -1 and k bmp22 =3e-02 min -1 are used. Other parameters are same as depicted in SI Text. Rest of the parameters. We were unable to track SN 1 in few cases (as we were unable to locate the position of SN 1 while performing the numerical analysis by using XPPAUT). Orange bar signifies movement of the saddle node (SN 1 ) towards higher BMP2 and dark green bar signifies movement of the saddle node (SN 1 ) towards lower BMP2 than the WT CNS case. Parameters are increased individually at an amount of 20% of the model parameters (SI Text) keeping all other parameters constant. Right panel shows bifurcation diagram of total Mash1 protein as a function of BMP2 at k bmp22 =5e-02 min -1 and k dmpid =4.5e-01 min -1 . (c) Left and middle panels show sensitivities of k bmp22 and k hsynp towards SN 1 and SN 4 in CNS. Right panel shows bifurcation diagram of total Mash1 protein as a function of BMP2 at k bmp22 =4e-02 min -1 and k hsynp =1.0 min -1 . Orange bar signifies movement of the saddle node towards higher BMP2 and dark green bar signifies movement of the saddle node towards lower BMP2 than the WT CNS case. Both the parameters are increased individually at an amount of 20% of the model parameters (SI Text) keeping all other parameters constant. In all the cases (Figs. S7a-c, right panels) BMP2 induces developmental cell fate change from gliogenic state (low Mash1 expression) to neurogenic state (high Mash1 expression).

Detailed description of Model construction
The model, schematized in Fig. 1 and Supplementary Fig. 1 This BMP2 driven self-inducing positive feedback loop of Mash1 promotes neuronal fate commitment in neural crest stem cells 1,3-5 . Also BMP2 stimulation is known to up regulate Id1 expression 1 . Activation of BMP2 leads to an increase in the expression level of Id1 1 . Id1 antagonizes Mash1 in two ways. Id1 competitively binds with E47 protein by sequestering E47 protein (corresponds the Id1 to E47 repression link in Module-1 of Fig. 1 (main text)) away from Mash1, which leads to the formation of transcriptionally inactive complexes and also it induces degradation (corresponds the Id1 to Mash1 repression link in Module-1 of Fig. 1 (main text)) of Mash1 monomer 2 . Mash1 protein stability is critically governed by E47/Id1 expression ratio 2 .
Keeping all the observations in mind module-1 (Fig. A1a) is constructed where we consider BMP2 stimulates both Mash1 and Id1 gene expressions at k bmp2 and k bmp21 rates respectively ( Supplementary Fig. 1). Mash1 mRNA (M m ) is transported from nucleus to cytoplasm. In cytoplasm Mash1 mRNA is designated as M C . Mash1 protein (M p ) forms heterodimer (ME) with E47. In ME, Mash1 undergoes CK2 mediated phosphorylation. The effect of CK2 mediated phosphorylation is incorporated in the rate constant k fme1 . The phosphorylated form (ME pc ) forms dimer ME pc2 . ME pc2 gets transported in to the nucleus and in the model it is abbreviated as ME pn , which contributes to the Mash1 gene activation (Supplementary Fig. 1). This implies that Mash1 auto-regulates its own expression in presence of BMP2. Similarly Id1 mRNA (I mn ) is transported from nucleus to cytoplasm. In cytoplasm Id1 protein (I p ) is translated from its mRNA (I mc ). I p sequesters E47 by making an inactive complex (IE) (corresponds the Id1 to E47 repression link in Module-1 of Fig. 1 (main text)). I p also activates the degradation of M p at k dmpid rate ( Supplementary Fig. 1) (corresponds the Id1 to Mash1 repression link in Module-1 of Fig. 1 (main text)).
Thus Id1 negatively regulates Mash1 expression in two possible pathways mentioned above. The equations for module-1 (Fig. A1a) are given in A1 Table (parameters values used are given in SI Text).
Hes6 also promotes proteolytic degradation (corresponds the Hes6 to Hes5 repression link in Module-2 of Fig. 1 (main text)) of Hes5 protein similar to Hes1 7,12 . Hes5 forms non-functional heterodimer with E47 and inhibits (corresponds the Hes5 to E47 repression link between Module-2 to Module-1 of Fig. 1 (main text)) its functional activity 1 . We assume Hes6 suppresses Hes5 to form complex with E47 similar to Hes1 25 . Pro-neural protein Ngn positively regulates Hes6 expression in committed neuronal precursors 6,7 . Pro-neural genes not only promote neuronal differentiation but also inhibit astrocyte-specific gene expression 11 . Therefore, our model also accounts for the transcriptional repression of Hes5 gene by Ngn (Details are shown in Supplementary Fig. 1).
We consider Hes5 transcription occurs in nucleus from both inactive gene (G 5i ) and active gene (G 5a ). Hes5 mRNA (H5 mn ) is transported from nucleus to cytoplasm. In cytoplasm Hes5 protein (H5 p ) is translated from Hes5 mRNA (H5 mc ) in cytoplasm.
H5 p forms a dimer (H5 p2 ). H5 p2 binds with another protein Gro/TLE and forms a complex (HP 2 ). Then hyper-phosphorylation of Gro/TLE takes place in the complex.
We consider four times phosphorylation of Gro/TLE in HP 2 . Then the hyperphosphorylated complex is transported to the nucleus. It is termed as HP 4n . Hes5 represses its own transcription (corresponds the self repression link of Hes5 in Module-2 of Fig. 1 (main text)) through HP 4n similarly as Hes1. HP 4n also inhibits (corresponds the Hes5 to Ngn repression link in Module-2 of Fig. 1 (main text)) Ngn gene activation. Ngn mRNA (N m ) is transported from nucleus to cytoplasm. In cytoplasm Ngn protein (N pc ) is translated from Ngn mRNA (N c ). N pc is transported from cytoplasm to nucleus. The Ngn protein in nucleus is labeled as N pn . N pn promotes Hes6 gene activation and inhibits (corresponds the Ngn to Hes5 repression link in Module-2 of Fig. 1 (main text)) Hes5 gene activation. Hes6 mRNA (H6 mn ) is transported from nucleus to cytoplasm. In cytoplasm Hes6 protein (H6 p ) is translated from Hes6 mRNA (H6 mc ). H6 p inhibits the formation of HP 2 . We consider H6 p also induces the degradation (corresponds the Hes6 to Hes5 repression link in Module-2 of Fig. 1 (main text)) of H5 p similarly as it does in case of Hes1. H5 p forms a nonfunctional dimer (EH 5C ) with E47. H6 p inhibits Hes5 to form EH 5C . We consider BMP2 stimulates Hes5 gene expression at k bmp22 rate. The effect of Notch stimulation on Hes5 is incorporated in the rate constant k hdel .
After constructing the gene interaction network for module-2, we couple module-2 with module-1 with the following additional interactions. ME pn induces gene activation of Ngn and Hes6, whereas HP 4n suppresses (corresponds the Hes5 to Mash1 repression link between Module-2 to Module-1 of Fig. 1 (main text)) gene activation of Mash1. The equations for the extended model (Fig. A2a) are given in A2 Table (parameters are given in SI Text). All the interactions shown in module-1 and module-2 are modeled as mass-action kinetic terms as shown in A2 Table except the Hes6 mediated degradation of Hes5.
Interestingly, the bifurcation diagram shown in Fig. A2b contains only the second bistable region in the steady state level of total Mash1 protein as a function of BMP2. Fig. A2c shows that at a pre-defined low BMP2 (BMP2=2 s.u.), the steady state level of total Mash1 remains in the high expression state and total Hes5 remains in the low Hes5 expressing state indicative of a neuronal state. Now at a pre-defined high BMP2 (BMP2=20 s.u.) dose, the steady state level of total Mash1 deterministically stays at the low expressing state and total Hes5 reaches the higher expression level representative of a gliogenic like state. Thus increase in BMP2 causes a decrease in the total Mash1 level with consequent up-regulation of total Hes5 level, which indicates acceleration of gliogenesis at high BMP2. Further, one can observe that the system can even show an oscillatory dynamics at the low Hes5 expressing state with a period of nearly 2 hr (Fig. A2c) but in the high Hes5 expressing state the oscillations are not present (Fig. A2c inset). The overall dynamics of the NSC's now is quite close to what has been observed for NSC's in CNS.

The complete model after including module-3
Dynamical regulation of negative HLH factor Hes1 is illustrated in module-3 ( Fig. 1 This whole module-3 is added to the model already generated by coupling module-1 and module-2 to obtain the full model ( Supplementary Fig. 1) by incorporating the following additional interactions. In our model, P 4n inhibits Mash1 (corresponds the Hes1 to Mash1 repression link between Module-3 to Module-1 of Fig. 1 (main text)) and Ngn (corresponds the Hes1 to Ngn repression link between Module-3 to Module-2 of Fig. 1 (main text)) gene activation as Hes1 known to inhibit pro-neural genes.
Similarly, Ngn protein (N pn ) inhibits (corresponds the Ngn to Hes1 repression link between Module-2 to Module-3 of Fig. 1 (main text)) Hes1 gene activation. H6 p inhibits (corresponds the Hes6 to Hes1 repression link between Module-3 to Module-2 of Fig. 1 (main text)) the formation of GP 2 . H6 p also induces the proteolytic degradation (corresponds the Hes6 to Hes1 repression link between Module-3 to Module-2 of Fig. 1 (main text)) of H1 p . H1 p forms a non-functional dimer (EH 1C ) with E47 and H6 p inhibits Hes1 to form EH 1C .
The governing equations, abbreviated names of different species involved in the network and parameters for the full model are depicted in SI Text and it seems to satisfactorily describe the WT as well as the mutant and over-expression phenotypes for NSC's in CNS (described in the main text). In all the gene knockout cases (Fig. 3,   Supplementary Fig. 3) we use total gene=0 for the corresponding gene that is knocked out and for the over-expression cases (Fig. 3, Supplementary Fig. 3) we over-express the corresponding total gene five times in comparison to the WT situation. Other parameters are same as provided in SI Text.

The model describes a number of experimental phenotypes for NSC's in PNS
To characterize BMP2 driven neuronal differentiation in PNS we have used the same network architecture (Supplementary Fig. 1  The model predicts the possible routes to reconcile the developmental

features of WT NSC's in PNS
Our systematic sensitivity analysis revealed that for few other pairs of parameters, we could also get a PNS like behavior (Other than shown in Fig. 4). The model predicts that for the following parameter values provided in A4 Table one can reconcile a PNS like situation as shown in Supplementary Fig. 7. Other parameters are same as SI Text (mentioned in the description of the parameters and their values and sources for CNS section). Supplementary Fig. 7a (left and middle panels) demonstrates that increase in k mp (translation rate of Mash1 protein) has more effect on the saddle node SN 1 than the saddle node SN 4 . The fact that k bmp22 and k mp act in an opposite way towards SN 1 and SN 4 , immediately suggests that both of them can be tuned to get the PNS like feature as shown in the Supplementary Fig. 7a (right panel). Thus, lowering the values for both k bmp22 and k mp reinforces neurogenesis (high Mash1 expression) at high BMP2 level (BMP2=20 s.u.) and ensures gliogenesis (low Mash1 expression) at lower values (BMP2=2 s.u.) of BMP2 doses. Supplementary Fig. 7b (left panel and middle panels) seems to suggest that both the saddle nodes SN 1 and SN 4 are more sensitive towards increase in k dmpid (Id1 mediated degradation rate of Mash1 protein) in comparison to k bmp22 . This again indicates that by decreasing k bmp22 to a sufficient extent and by slightly increasing k dmpid , one can enforce neuronal differentiation at high BMP2 (i.e., at BMP2=20 s.u.) and consequently onset of gliogenesis can happen at much lower (BMP2=2 s.u.) BMP2 levels (PNS like situation in Supplementary Fig. 7b (right panel)). In similar note to what we observe in Supplementary Fig. 7b, Supplementary Fig. 7c (left panel and middle panels) evidently depicts that both the saddle nodes SN 1 and SN 4 are relatively more sensitive towards increase in k hsynp (translation rate of Hes5 protein) in comparison to k bmp22 . Thus a sufficient decrease in k bmp22 and a moderate increase in k hsynp create a PNS like feature as shown Supplementary Fig. 7c (right panel). More experiments in this direction will be helpful to decipher this issue more clearly.

A4
On the other hand, one can easily tune the parameters k mp (translation rate of Mash1 protein), k dmpid (Id1 mediated degradation rate of Mash1 protein) and k hsynp (translation rate of Hes5 protein) experimentally and verify these model predictions.