Developmental Emergence of Sparse Coding: A Dynamic Systems Approach

During neocortical development, network activity undergoes a dramatic transition from largely synchronized, so-called cluster activity, to a relatively sparse pattern around the time of eye-opening in rodents. Biophysical mechanisms underlying this sparsification phenomenon remain poorly understood. Here, we present a dynamic systems modeling study of a developing neural network that provides the first mechanistic insights into sparsification. We find that the rest state of immature networks is strongly affected by the dynamics of a transient, unstable state hidden in their firing activities, allowing these networks to either be silent or generate large cluster activity. We address how, and which, specific developmental changes in neuronal and synaptic parameters drive sparsification. We also reveal how these changes refine the information processing capabilities of an in vivo developing network, mainly by showing a developmental reduction in the instability of network’s firing activity, an effective availability of inhibition-stabilized states, and an emergence of spontaneous attractors and state transition mechanisms. Furthermore, we demonstrate the key role of GABAergic transmission and depressing glutamatergic synapses in governing the spatiotemporal evolution of cluster activity. These results, by providing a strong link between experimental observations and model behavior, suggest how adult sparse coding networks may emerge developmentally.


Results
To study developmental changes in postnatal neuronal activity, we used simulations of a spatially localized recurrent neural network with short-term synaptic plasticity (STP-RNN) 13 , where STP renders synaptic efficacies dynamic over time 14,15 ; see Methods. This network is a mean-field Wilson-Cowan-type model of one excitatory (E) and one inhibitory (I) neuronal population with recurrent dynamic synaptic connections (Fig. 1a). This model has the advantage of being biophysically interpretable and mathematically accessible. The activity rates (E r and I r ) can be properly scaled to represent locally the average recorded activities in the populations. Here, we use experimentally reported, postnatal changes in neurobiologically plausible parameters (Table 1) to assess how the network's spontaneous behavior and sensory processing properties are refined toward the postnatal onset of sensory transduction. To study the developmental states of visual cortical networks from before to after eye-opening, we selected four postnatal days (P) for modelling: P3 (period of physiological blindness), P10 (a few days before eye-opening), P14 (the day after eye-opening) and P20 (a few days after eye-opening) 3 . During this period, immature networks not only undergo the sparsification process but also a dramatic development in intrinsic neuronal and STP properties 10,12 ; see Table 1 and Supplementary Methods. To derive the results, we used system dynamics methods described in the Methods section (see also Supplementary Methods). For quick reference, all important technical terms are explained in Table 2.

Transient unstable state hidden in firing activities.
We first address what mechanisms may underlie the emergence of cluster activity in immature networks, by analyzing model behavior at P3 as a representative early postnatal stage.
Previous experimental 16,17 and modeling studies 18,19 usually considered adult networks to be bi-or multi-stable where, e.g., spontaneous (stimulus-absent) network activity can transition between two stable activity states: A relatively quiet activity state, and a higher activity state. For early postnatal stages, our stability analysis and simulations indicate that the proposed developing network model is actually mono-stable whose only stable state is a quiescent, spontaneous state (Fig. 1b). To demonstrate this, we plotted in Fig. 1b the so-called phase plane of network dynamics in terms of the average activity in the E-and I-populations (E r -I r -plane; Table 2). Figure 1b shows that the network's fixed point (FP; Table 2) is located at E r = I r = 0 which is essentially a stable activity state of the network (see Supplementary Methods). In the following, we will call this specific FP the rest state of the network. This analysis result fits well with experimental findings of near-zero hertz spontaneous activity of immature populations 3,20 , and the absence of spontaneous network's persistent activity states mainly during periods of physiological blindness 7,21 .
How can a quiescent, immature, mono-stable network generate large cluster activity? By analyzing network perturbations, we found strikingly that while a P3 network model was relaxed at its stable rest state, threshold-crossing perturbations of E-population were amplified profoundly (Fig. 1c, Supplementary Fig. 1c). This amplification resulted in population spikes (PSs , Table 2), which, in the model, we assume is the analogue to experimentally observed cluster activity. Note that this perturbation can be due to, e.g., a single-shock electrical stimulation 22 , the onset of a longer-lasting external input ( Supplementary Fig. 2c), or a random input driven by spontaneous retinal waves or the thalamus 3,20 .
While previous studies showed the underlying mechanism of PSs generation mainly in adult bi-stable network models 18,19 , here we sought to address this mechanism in an immature mono-stable network (Fig. 2). As the initial phase of PSs is known to be governed by the network's fast (i.e. firing activity) dynamics 18 , we re-plotted the E r -I r -plane after freezing the network's slow (i.e. STP) dynamics at the rest state (Fig. 2d); thus, it is turned to a Static-RNN (Table 2) with frozen synaptic efficacies (Frozen STP-RNN). Surprisingly, we found that in addition to the stable rest state, there is an unstable FP hidden in the network's fast dynamics (i.e. in the respective Frozen STP-RNN, Fig. 2e), which does not exist in the (non-frozen) STP-RNN (Fig. 1b). This unstable FP (black dot, Fig. 2e), which is located close to the rest state, can confine the attraction domain (grey region, Table 2) of the stable rest state in the E r -I r -plane of the corresponding Frozen STP-RNN (Fig. 2d); but, clearly, not in the E r -I r -plane of the STP-RNN itself (Fig. 2c). We call the outer border of this domain in Fig. 2d (before perturbation) the amplification-threshold. Almost exactly the same border exists for activity-perturbation domains of the STP-RNN (Fig. 2c). Cluster activity is initiated by an activity perturbation exceeding this threshold, where the activity will be moved far from the rest state in both corresponding networks with frozen (Frozen STP-RNN) or dynamic (STP-RNN) synapses. Strikingly, when freezing synaptic efficacies at different times during a cluster activity (Fig. 2a), we can see that this unstable FP disappears roughly after the peak of cluster activity (Fig. 2e), thereby allowing network activity to converge back to its stable rest state in the STP-RNN (Fig. 2a). This disappearance is mainly because of the weakening of the excitatory synapses at high activity rates (Fig. 2b); these synapses provide positive feedback to the network. Importantly, as it can be seen in Fig. 2e  network is again at the origin after cluster activity, sufficient recovery time is required for the re-emergence of the unstable FP.
Our analyses further show that the experimentally observed variability in size (i.e. the number of network neurons recruited by cluster activity) and duration of spontaneous cluster activities at each postnatal day 2,3 is effectively governed by the relatively long recovery time of immature excitatory synapses ( Supplementary Fig. 1). In addition, in our model, cluster activity was generated only if the perturbation strength was sufficiently large (Supplementary Fig. 1c and d) and was able to push network activity beyond the amplification-threshold ( Fig. 2c  and d). This is consistent with experimental reports of immature networks 7,23 , and may underlie a relatively all-or-none characteristic of postnatal cluster activity.  Table 2  The route towards sensory processing. The sparsification of activity patterns sets in around eye-opening 3,7 . We now use our model at the stages P10 to P20 to address the question what refinements of the developing network enable sparsification. Firstly, similarly as for P3, at all later stages up to P20 in the model we found that for a developing network relaxed at its stable rest state (Fig. 3a), again an unstable FP is hidden in the network's fast (i.e. firing activity) dynamics (Fig. 3b). Consequently, a threshold-crossing perturbation of the E-population at rest state led to cluster activity (Fig. 3c). In addition, our analysis revealed that the amplification-threshold tends to move toward higher E-activity rates, from P3 to P20 (Fig. 3b, Supplementary Fig. 3). This result suggests that a stronger input may be required to trigger (large) cluster activities at late development stages.
Experimental in vivo findings show that the cluster activity size starts to drop dramatically around the end of the second postnatal week 2,3 . Here, we consider this developmental decrease in cluster activity size as an indicator for the developmental transitioning from dense to sparse coding (sparsification process); similarly to experimental studies 2,3 . In our model, we approximate the cluster activity size by a measure of the average activity of both E-and I-populations (PS net amp ; see Table 2 and Methods). Using this measure, we observed, similarly as in experiments, that cluster activity size undergoes developmental reduction between P10 to P20 ( Fig. 3c and  ). Secondly, the results also show that immature networks probably lack any spontaneous, non-quiescent stable FP (i.e. a spontaneous attractor or persistent activity state) up to a few days prior to eye-opening. In Fig. 3a this is indicated by the existence of only one green dot at P3 and P10. However, after eye-opening, new spontaneous FPs with higher activity rates than the rest state emerged (note the non-origin FPs at P14 and P20 in Fig. 3a). This is compatible with experimental data showing the existence of spontaneous persistent cortical activity for stages mostly after the second postnatal week 3,21,24 and their absence during early development 7,21 . In our model, in the absence of a stimulus, these FPs emerge mainly because of the developmental increase and decrease in background activity 2,3,7 and membrane resistance 2,10 , respectively. Their combined effects determine the increase in population activity thresholds θ E and θ I (see also Supplementary Fig. 4). While in the model the network's FP domains (FP-domains; Table 2) are not affected by these thresholds (see "Characterization of operating regimes" section in Supplementary Methods and also Supplementary equation (17)), the developmental reduction in the network instability (yellow regions, Fig. 3a) allows for these FPs to be stable. This means that attractors start emerging, under a wide range of E-activity; note the non-origin green dots at P14 and P20 in Fig. 3a.
Thirdly, the emergence of spontaneous attractors not only turns the developing mono-stable network (at P3 and P10) into a bi-stable (at P14 and P20) network (Fig. 3a), but also has a striking effect on the dynamical trajectory of cluster activity. That is, in the model, we found a fundamental difference between spontaneous postnatal cluster activities and those in more mature networks (Fig. 3c): before eye-opening (P3 and P10), the cluster activity is of the so-called mono-stable type (see Supplementary Methods), while after eye-opening the cluster activity can be also of the bi-stable type. The difference is that a mono-stable cluster activity (P3 and P10 in Fig. 3c) is initiated and terminated at the stable rest state, while a so-called bi-stable cluster activity is initiated at the stable rest state but converges to a spontaneous attractor (P14 and P20 in Fig. 3c  interpreted as a first expression of information processing where the newly emerged attractors, which a bi-stable cluster activity converges to, can be seen as representative states that are the basis of perception. In the model, the time required for an effective transition to the attractor, and thus an effective representation of perceptual stimuli, is probably equal to the duration of cluster activity (see dashed grey lines, Fig. 3c): ~240 ms at P14 and ~180 ms at P20; these transition times are consistent with experimental observations at these stages 17,24 .
Moreover, our simulations showed that the spatiotemporal characteristics (i.e. size and duration) of mono-stable cluster activities are in general considerably more robust to interfering perturbations (P3 and P10, Supplementary Fig. 5), as compared to the bi-stable ones (P14 and P20, Supplementary Fig. 5). This suggests that the trajectories underlying mono-stable cluster activities (P3 and P10) can be seen as signals with high signal-to-noise ratio, which may act as a reliable neuronal communication mechanism in neonatal networks. The decreased robustness of bi-stable cluster activities after eye-opening (P14 and P20) can in principle augment the network's information processing capability by making it more flexible in responding to sensory input.
In addition, we found that at the new spontaneous attractors (non-origin green dots; Fig. 3a) the network operates as a so-called inhibition-stabilized network (ISN, light-yellow regions, see also Table 2) 25, 26 . An ISN regime is thought to allow the cortex to process complex computations 27 , and some experimental evidence indicates that adult cortical networks operate as ISNs 25,28 . Our results show that the ISN regime becomes effectively accessible after eye-opening, as the unstable FP-domain, which is confined to low E-activity ranges, is developmentally substituted by an ISN FP-domain. This indicates a developmental increase in the ratio of areas of the ISN (AOD ISN ) and unstable FP-domains (AOD Unstable ), which we quantified by AOD ISN/Unstable (Table 2) in Fig. 3e. In practice, this means that in parallel to the strengthening of the sensory inputs after eye-opening, a larger set of potential stimulus-evoked ISN attractors will be accessible for the network, presumably to perform complex sensory computations. Overall in our model, the process of sparsification is translated not only to a potent reduction of postnatal network instability but also to a potential emergence of new attractors and information processing capabilities.
Spatiotemporal evolution of postnatal cluster activity. What mechanisms govern the spatiotemporal evolution (i.e., size and duration) of cluster activity during development in vivo? We addressed this question by analyzing how a postnatal cluster activity, once on a trajectory away from the rest state, is brought back to either the rest state or another FP with an activity level higher than the rest state.
In our model, there are two important factors: (i) the inhibitory transmission, which is believed to play a critical role in the stabilization of adult networks 25, 28 , and (ii) STP, which can dynamically control the gain of neuronal responses 14 and acts strongly depressing in immature networks 10,12 .
Firstly, we found that cluster activity initiated at the rest state can still converge back to this state even when blocking GABAergic receptors, at all stages from P3 to P20 (Fig. 4a). This finding suggests that inhibitory E r -I r -plane A 2D phase plane of E-and I-population activity rates (I r vs. E r ) ‡ .

FP
A fixed point (FP) is the steady state of a network, and is determined as the intersection of (quasi) E rand I r -nullclines in the E r -I r -plane ‡ .
Spontaneous activity Network activity in the absence of stimulus (or sensory input)*.
Attractor A non-quiescent stable FP (in this paper); also known as memory or persistent activity state.
Attraction domain For a stable FP (in this paper), a region in phase plane comprising all initial conditions that lead to that FP.

Amplification-threshold
The outer border of the attraction domain of the rest state in the Frozen STP-RNN, beyond which network perturbations undergo an overall continuous growing (in the Frozen STP-RNN) or will effectively trigger cluster activity (in the STP-RNN).

Nullcline
For example, an excitatory nullcline is a set of points (a curve) in E r -I r -plane for which

FP-domain
For example, an unstable FP-domain is a domain of all potential FPs in the E r -I r -plane, at which the network will be unstable ¶ . The scaled amplitude of PS net , which provides a qualitative approximation to the cluster activity size # . processing may not be necessary for stability in the developing network with depressing excitatory synapses. However, blockage of GABAergic receptors yielded an increase in cluster activity size ( Fig. 4a and c), consistent with in vivo reports 20,29 . This blockage effect on cluster activity size became more pronounced during the course of development (Fig. 4c), in spite of the developmental reduction in the network instability (e.g. see AOD ISN/Unstable in Fig. 3e). This finding may underlie a developmental enhancement of an effective contribution of inhibitory transmission to network activity (see also Supplementary Fig. 2d showing the ability of inhibitory inputs to turn off activity at the attractor 22 , after eye-opening). In addition, we found that the blockage of GABAergic receptors shortened the cluster activity duration, for example, from approximately 330 to 320 ms at P3 (Fig. 4a and b) and from approximately 265 to 210 ms at P10 (Fig. 4a). The range of these cluster activity durations in our model is compatible with previous experimental reports 7,30 . Moreover, we found that this blockage effect on cluster activity duration increased from P3 to P20 (Fig. 4d). Secondly, in contrast to GABAergic transmission, we found that at P3 the simulated removal of the STP effect led to run-away excitation (Fig. 4e). After freezing the synaptic efficacies (see Supplementary Methods) at the stable rest state, a threshold-crossing perturbation of the E-population (similarly as in Fig. 1c) caused a seizure-like surge of network activity (Fig. 4e). Importantly, this means that the rest state in the corresponding Frozen STP-RNN is a locally, but not globally, stable FP (see Methods). Similar results were observed for P10 to P20 (data not shown). Accordingly, we conclude that the strongly depressing immature STP, but not inhibitory synaptic transmission, guarantees the stabilization of postnatal cluster activities during the course of development. This finding also holds, when GABAergic transmission is considered as excitatory at the network level during the first postnatal week (see Supplementary Methods and Supplementary Fig. 6b). Under this assumption, blocking GABAergic receptors led to a decrease in cluster activity size, which is not compatible with in vivo data 20,29 . Therefore, our modeling results (Fig. 4b) are in agreement with recent in vivo findings that while GABA  r max r max ) at [10,10] Hz. Moreover, the developmental decrease in − AOD Non ISN just means that the E-activity rate after which the Non-ISN FPdomain starts in the E r -I r -plane, was shifted to higher levels. See Table 2 for details about technical terms.
SCiEntifiC REPORTS | 7: 13015 | DOI:10.1038/s41598-017-13468-z acts as a mainly depolarizing transmitter during the first postnatal week, it exerts an inhibitory effect on the underlying network activity in the immature cortex and hippocampus 20,31 .
Moreover, cluster activity was abolished by simulated blockage of glutamatergic receptors (data not shown), consistent with experimental in vivo data 20,29 . This is because the main source of the network instability in our model is determined by glutamatergic synapses which provide positive feedback to the network. Therefore, blocking glutamatergic receptors removed all unstable FPs as well as the attractors in Fig. 3a and forced spontaneous network activity to the rest state (as observed experimentally 17,22 ), and removed the unstable FP hidden in the network's fast dynamics in Fig. 3b, throughout P3 to P20 (data not shown). In addition, note that the GABAergic transmission can play an important modulatory effect on this instability (see Supplementary Methods), where an increase in the efficacies of inhibitory synapses can effectively attenuate the instability (not shown). This implies that relatively weak GABAergic inhibition (Table 1) is permissive for the generation of the unstable states (e.g. see Fig. 3a and b) in developing networks; note that in our model these synapses strengthen during development (Table 1).
In sum, we found that strongly depressing excitatory synapses have a key role in the termination of postnatal cluster activities, whereas GABAergic transmission mainly regulates their spatiotemporal evolution.
Key maturational processes mediating sparsification. The most influential, maturational processes mediating sparsification are still not clearly understood 2,3 . By using our modelling approach, we mechanistically quantified the impact of different network parameters on the emergence of sparse coding.
To do this, we replaced single parameters or parameter combinations of the P10 network model by their respective values at P20 and measured to what extent this modification can account for the normal decrease in PS net amp and normal increase in AOD ISN/Unstable , when transitioning from P10 to P20. To quantify effects we computed ratio PS and ratio AOD as the ratios of the modification-induced changes in PS net amp and AOD ISN/Unstable , relative to their respective normal developmental changes during this period (Fig. 5). We found that only three parameters caused a substantial decrease in PS net amp (Fig. 5a): the two absolute synaptic efficacies (J E and J I ) caused the largest decrease close to that in the normal transition and, to a minor degree, the excitatory release probability U E . Strikingly, for the excitatory synaptic time constant (τ E ) the modification led to an increase in PS net amp , i.e. when this parameter is changed on its own, it tended to reverse sparsification (see also Supplementary Figs 6c and 7). We obtained qualitatively similar conclusions when analyzing the effect of replacing single parameters or parameter combinations of the P20 network model by their values at P10 (Supplementary Fig. 8a). When testing for changes in AOD ISN/Unstable , we found that the inhibitory synaptic depression time constant τ r I had the strongest contribution to the normal transition from P10 to P20, whereas all other single parameter substitutions had relatively weak effects (Fig. 5b). Note that the change in τ r I appears to be a necessary but not a sufficient condition for the normal transition from P10 to P20, as substituting J E and J I in the P20 model by their respective P10 values virtually abolished the developmental change in AOD ISN/Unstable (Supplementary Fig. 8b). Figure 5 further shows that the modification effect of the population activity thresholds (θ E and θ I ) was virtually negligible for ratio PS and zero for ratio AOD . However, as shown above, their maturation plays an important role for the emergence of spontaneous attractors (see also Supplementary Fig. 4).

Discussion
We modelled the in vivo activity of a developing cortical network during the first postnatal month by combining an extended Wilson-Cowan model with experimentally reported trajectories of neuronal and synaptic parameters. We revealed mechanistically that a particular combination of a transient, hidden unstable state in firing activities and strong synaptic depression enables an immature network to generate large cluster activity while otherwise being mostly silent. We further found that the normal developmental transition from dense to sparse coding is strongly dependent on an elaborate, parallel refinement of absolute synaptic efficacies, both short-term synaptic plasticity (STP) and intrinsic membrane properties, and background activity. Strikingly, in our model, sparsification translates not only to a reduction of postnatal instability of network activity but also to an effective availability of the inhibition-stabilized network (ISN) regime as well as the emergence of spontaneous attractors, providing a novel mechanistic explanation for how the network's information processing is refined towards eye-opening.
How can immature networks be quiescent for relatively long periods and occasionally generate large cluster activity 3,32 ? Surprisingly, we found that while the developing networks (P3 to P20) are operating at their rest state (Fig. 3a), an unstable state is formed in their fast (i.e. firing activity) dynamics (Figs 2e and 3b). This may be a key to this biphasic behavior ( Fig. 3c and Supplementary Fig. 2c). In addition, we found that at early stages prior to eye-opening (P3 and P10) immature networks are mono-stable, where the only FP of the network is its stable rest state (Fig. 3a). Accordingly, the underlying mechanism of cluster activity emergence in these networks is in stark contrast to models of cluster activity proposed for adult networks, based on the usual assumption of a bi-stable (or multi-stable) network 18,19 . Our finding about the strong effect of the hidden unstable state (Fig. 3b) during the initial phase of network activity at the rest state (see Figs 2 and 3c) may provide an explanation why neonatal networks, e.g. at P3 and P10, are more susceptible to seizures than mature networks 32,33 . Moreover, the lack of any spontaneous attractor in these networks (Fig. 3a) might also contribute to this susceptibility, since the attractors can, in principle, aid in stabilizing a network's activity.
Sparsification coincides with the peak of GABAergic and glutamatergic synaptogenesis 11,12 and overlaps with a developmental reduction in release probabilities, particularly at glutamatergic synapses 10,12 . Using a computational study, we here revealed how these changes are suited to drive the transition from dense to sparse coding during network maturation (e.g. see Figs 3 and 5). Since sensory deprivation affects the total synapse numbers only modestly 34 , our results might also explain previous experimental findings: While sensory experience may have a modulatory effect on sparsification during the first days after eye-opening 3 , this process is largely mediated internally, i.e. independent of sensory inputs 2,3,7 .
How does the sparsification process prepare in vivo developing networks for effective sensory processing? While the answer to this question remains poorly understood by experimental studies, our computational study provides two mechanistic insights: Firstly, we found that Non-ISN and instability may play key roles as the dominating operating regimes prior to eye-opening (Fig. 3a), possibly reflecting an immaturity of inhibitory transmission. That is, our modeling results imply that during stimulation (e.g. in response to long-lasting external inputs) the immature networks will operate under a Non-ISN regime, rather than an ISN regime (P3 and P10; Supplementary Fig. 2b), due to the lack of an effective availability of the ISN regime during early development (P3 and P10; Fig. 3a). The Non-ISN regime will enable the immature networks to maintain their stability without requiring inhibitory transmission effects (see Methods). For adult networks, however, some experimental evidence indicates that cortical networks operate as ISNs 25,28 . Strikingly, in the model, during sparsification, the unstable FP-domain is effectively replaced by that of an ISN regime (Fig. 3a and e). This may enable cortical networks to process complex sensory computations 27 in parallel to the developmental strengthening of sensory inputs.
Secondly, we found that spontaneous attractors start emerging around eye-opening, which, in our model, is mainly due to the combined effect of developmental increase in background activity 2,3,7 and the developmental decrease in membrane resistance 2,10 . In general, attractors can represent the solution to a specific neural computation 35 . They are also thought to be the substrate of, and used to model, e.g., the working memory 36 , eye position stability 37 , and orientation selectivity 38 . Besides, in our developing model, the attractor emergence phenomenon can also underlie three developmental mechanisms. First, it provides the basis for transitioning between FPs (here, via bi-stable cluster activity at P14 and P20; Fig. 3c) and may be thought of as an early step towards representation of perceptual stimuli. Second, it enables inhibitory transmission to contribute effectively more to network activity, e.g. through providing a sufficiently strong balancing of the intrinsic instability of the E-activity under an ISN regime 26,27 (P14 and P20; Fig. 3c), or by becoming amenable to terminate activity at the attractor 22 ( Supplementary Fig. 2d). This may be important for a presumably more effective processing of sensory inputs. Third, this phenomenon may also initiate the effective interaction of spontaneous activity with sensory cortical responses, as observed experimentally 17 .
Evidence from in vitro studies suggests that GABA depolarizes neonatal neurons, e.g. in rodents during the first postnatal week 20,32 . However, the effect of GABA at the network level is still an open question. While most of the previous in vitro studies reported that GABA is excitatory (or both excitatory and inhibitory) at the network level 32 , recent in vivo studies harnessing new experimental techniques found that GABA inhibits intact neonatal networks in neocortex and hippocampus 20,29,31 . In particular, GABAergic transmission was shown to limit the spatial extent of cluster activity 20 . In our model, we tested these two hypotheses by blocking the GABAergic receptors in two P3 networks with either excitatory (Supplementary Fig. 6b) or inhibitory GABAergic transmission (Fig. 4b). This simulated manipulation led to an increase in the cluster activity size only when GABAergic transmission was inhibitory at the network level. Therefore, our results support the recent in vivo findings 31,33 . Importantly, this finding implies that even when GABA depolarizes immature neurons, GABA still can attenuate the instability effect on network activity, thereby restricting cluster activity size.
To our knowledge, there are only few previous modelling studies covering this developmental period (e.g. see refs 39,40 ), where the authors' main focus was on the potential mechanism of the cluster activity generation under in vitro condition. These studies were focused on a single postnatal stage and this stage was mostly restricted to the first postnatal week, like P5. Moreover, these authors probably did not consider some important biophysical considerations which have only recently been revealed by in vivo studies: the abolishment of cluster activity following the removal of glutamatergic synapses 20,29 , the existence of a profound inhibitory effect of GABAergic transmission at the network level 20,29 and its non-excitatory (but depolarizing) effect at the neuron level 20,31 during the first postnatal week, or that the loss of large cluster activity occurs around eye-opening 2,3 and not around P7 41 . In contrast, here, (i) we used a unified computational model which takes into account most of these recently reported biophysical considerations (see Supplementary Method and Table 1), (ii) our work shows results based on, and mainly in accordance to the recent in vivo studies, where these studies revealed some critical features different from those reported for in vitro data, (iii) we studied the whole developmental period, comprising the time course of sparsification, and not only a single postnatal stage, (iv) we not only address the possible underlying mechanism of the postnatal cluster activity generation, but also the developmental changes in the operating regimes and information processing capabilities, as well as the possible refinements in network properties governing the sparsification.
There are advantages and limitations of our modeling approach. Clearly, as compared to a mean-field model as employed here, a spiking network model can in principle provide more detailed results, e.g. one could model the temporal sparseness of spiking activities, or the specific patterns of sensory inputs at the neuron level. However, for this early developmental period, we are not aware of any established neuron model. For such a model, a considerable amount of unspecified neurobiological parameters is required to be measured experimentally first. Instead, using a mean-field model enabled us to forgo such detailed parameterizations by considering only their average characteristics over the network. Although this procedure comes at the price of removing biological details, this enabled us to use an extended Wilson-Cowan-type model 13 ; a well-established model and extensively used before to study adult networks behavior 18,42 . This model enabled us to incorporate the available, experimentally reported developmental trajectories of intrinsic neuronal and synaptic parameters (see Supplementary  Method and Table 1), including short-term synaptic plasticity. In addition to its biophysical interpretability, this model is in general mathematically tractable which allowed us to derive analytical expressions, e.g. in our stability analyses (see Supplementary Methods). Moreover, this model is readily extendable to incorporate other biophysical mechanisms like spike-frequency adaptation, or to build a spatial graph of multiple homogenous and heterogeneous inter-connected networks to study, e.g. the activity propagation over different brain areas.
In sum, we have shown that by establishing and extending a novel application of an existing computational model to immature networks and integrating recent experimental findings, new mechanistic insights into the development of neural networks can be obtained. We expect that, in the future, this modelling approach can also guide research by providing for concrete predictions that can be tested experimentally.

Methods
Here, we briefly describe the main components of our model and the analyses. A detailed description can be found in Supplementary Methods. STP-RNN model. This model (Fig. 1a) is an extended version of Wilson-Cowan's recurrent neural network (RNN) mean-field model 43 , for which the short-term plasticity (STP) of synaptic connections were also modeled 13 . The equations governing the model dynamics over time are (dots denote the time derivatives): r , j is the index of presynaptic population, E r (=A E ) and I r (=A I ) are the average activity (in hertz) of E-and I-populations which receive, respectively, the external inputs e E and e I , e.g., from other brain regions, G i is the linear input-output gain above population activity threshold θ i (otherwise = G 0 i ), and x ij and u ij are the dynamic variables of short-term synaptic depression and facilitation mechanisms. The parameter definitions and values are listed in Table 1. A STP-RNN with synaptic efficacies frozen at the FP is called Frozen STP-RNN (a Static-RNN-type model).
Phase plane components. E r -I r -plane is the 2D phase plane of E-and I-population activity rates (I r vs. E r ), in which the E r -and I r -nullclines are represented as the set of points for which = dE dt / 0 r and = dI dt / 0 r , respectively. For a STP-RNN, the depicted quasi E r and I r nullclines in E r -I r -plane are based on the reduced STP-RNN (Supplementary equation (7)). Any intersection of these two quasi nullclines (i.e. the FP) can represent the steady state of the whole network (i.e. the 10D STP-RNN).

Stability of FPs.
To determine the stability of any FP of interest we applied linear stability analysis to our 10D STP-RNNs (or 2D Frozen STP-RNNs): We investigated whether all eigenvalues of the network's system of equations linearized around the FP (using the Jacobian matrix) have strictly negative real parts (if so, the FP is stable), or whether at least one eigenvalue with positive real part exists (if so, the FP is unstable). This can be re-stated as: A FP is stable if following a small perturbation from that FP (at which the network was before perturbation), the network dynamics converge back to it. Conversely, if the dynamics move away from the FP or dies out, that FP is unstable.
SCiEntifiC REPORTS | 7: 13015 | DOI:10.1038/s41598-017-13468-z Operating regimes and FP-domains. The stable operating regimes of a RNN at a FP can be classified as an inhibition-stabilized network (ISN) vs. a Non-ISN 25,26 . To discriminate between these two regimes three criteria were defined: (A) Excitatory instability: For the inhibitory activity rate fixed at the FP, the recurrent excitation is strong enough to render the E-population intrinsically unstable. (B) Excitatory stability: In contrast to (A), the E-population is stable per se, i.e. even with a feedback inhibition fixed at its level at the FP. (C) Overall stability: The dynamic feedback inhibition to the E-population is strong enough to stabilize the whole network activity. At a FP, a network operating under the (A) and (C) criteria is an ISN, while a network operating under the (B) and (C) criteria is a Non-ISN. A network, which is neither ISN nor Non-ISN at the FP, operates under an unstable regime.
We partition the E r -I r -plane into different domains of operating regimes (FP-domains). Each FP-domain contains all potential steady states (i.e. FPs) at which the network could operate under the corresponding regime. The area of each FP-domain is computed numerically as AOD (area of domain). Note that the borders between the FP-domains (Fig. 3a) were determined by using numerical simulations, when we plotted the FP-domains of operating regimes based on their stability criteria (A-C) obtained analytically in Supplementary Methods (see "Characterization of operating regimes" section in Supplementary Methods).
Cluster activity size and duration. In a RNN, the cluster activity (network spike; PS net ) usually involves the population spikes 44,45 in both E (PS E ) and I (PS I ) populations. To approximate cluster activity size, i.e. the total number of neurons synchronized during a cluster activity, we calculate ω = × − PS A A ( ) net amp sum amp sum 0 (a dimensionless parameter), where A sum = E r + I r is the sum of the activity of E-and I-populations (in hertz), A sum amp and A sum 0 are the maximal activity during the cluster activity and the preceding activity, respectively, and ω is a scaling factor (in units of [Hz −1 ]) to convert the activity rate to an approximate number of activated neurons during the cluster activity. For simplicity, we set ω = 1, since its veridical value is not defined for the mean-field model. We measure the cluster activity duration as its termination time minus its onset time.
Computation of ratio PS and ratio AOD . To investigate the impact of specific maturational processes on sparsification (as used in Fig. 5), we first substituted single parameters or a small combination of parameters at P10 by their values at P20 ( , where for ratio PS we set γ = PS net amp with ϖ = −1, and for ratio AOD we set γ = AOD ISN/Unstable with ϖ = +1. γ P10 res is the value of γ measured after the parameter(s) value(s) substitution.
Data availability. Custom Matlab and Mathematica codes for our model are available upon request from the corresponding author.