How Memory Conforms to Brain Development

Nature exhibits countless examples of adaptive networks, whose topology evolves constantly coupled with the activity due to its function. The brain is an illustrative example of a system in which a dynamic complex network develops by the generation and pruning of synaptic contacts between neurons while memories are acquired and consolidated. Here, we consider a recently proposed brain developing model to study how mechanisms responsible for the evolution of brain structure affect and are affected by memory storage processes. Following recent experimental observations, we assume that the basic rules for adding and removing synapses depend on local synaptic currents at the respective neurons in addition to global mechanisms depending on the mean connectivity. In this way a feedback loop between “form” and “function” spontaneously emerges that influences the ability of the system to optimally store and retrieve sensory information in patterns of brain activity or memories. In particular, we report here that, as a consequence of such a feedback-loop, oscillations in the activity of the system among the memorized patterns can occur, depending on parameters, reminding mind dynamical processes. Such oscillations have their origin in the destabilization of memory attractors due to the pruning dynamics, which induces a kind of structural disorder or noise in the system at a long-term scale. This constantly modifies the synaptic disorder induced by the interference among the many patterns of activity memorized in the system. Such new intriguing oscillatory behavior is to be associated only to long-term synaptic mechanisms during the network evolution dynamics, and it does not depend on short-term synaptic processes, as assumed in other studies, that are not present in our model.


INTRODUCTION
A complex interrelation between "form" and "function" is known to play an important role in nature (Gross and Blasius, 2008;Vazquez et al., 2008;Sayama et al., 2013). The idea has been efficiently developed in the field of adaptive networks, in which a sort of coupling feedback loop sets in between the network dynamic activity and its topological structure. Outstanding phenomena then emerge, including self-organization into complex topologies that exhibit robust dynamics, spontaneous differentiation of the nodes, or complex mutual dynamics in both activity and topology, in any case mimicking many different conditions in nature (Bullmore and Sporns, 2009;Sayama et al., 2013;Millán et al., 2018a). This framework has revealed quite useful to understand fundamental questions concerning mammal brains, e.g., how structural and functional properties relate to each other both at the level of models involving sets of neurons and synapses and at the coarse-grained scale of connectomes and functional nets which is captured by imaging techniques (Bullmore and Sporns, 2009). A main question that we can thus address is how an efficient brain develops by synaptic pruning after a sort of "wild" proliferation of synaptic connections between neurons following conception (Chechik et al., 1998;Iglesias et al., 2005;Santos and Noggle, 2011;Presumey et al., 2017). In humans, for example, synaptic density at birth is about twice that at puberty, and certain brain disorders, such as autism spectrum disorder (ASD) and schizophrenia, have been related to details of synaptic pruning (Keshavan et al., 1994;Geschwind and Levitt, 2007;Faludi and Mirnics, 2011;Kolb et al., 2012;Fornito et al., 2015). In particular, ASD has been associated with a defect of synaptic pruning in certain brain areas (Tang et al., 2014), whereas schizophrenia could be related to an excessive pruning (Sekar et al., 2016). In any case, it now seems clear that such synaptic pruning involves in some way an optimization process, probably aimed at minimizing both energy consumption and the genetic information that otherwise would be needed to build an efficient and robust network (Chechik et al., 1999;Chklovskii et al., 2004;Johnson et al., 2010;Knoblauch et al., 2010;Navlakha et al., 2015). In particular, recent studies on associative memory have shown that this process could greatly improve memory retrieval under a noisy environment, such as it is the case in biological systems (Millán et al., 2018a). Moreover, ongoing structural plasticity in the adult brain has also been suggested to improve substantially the storage capacity (Chklovskii et al., 2004;Knoblauch et al., 2010), and has been related to graded amnesia, catastrophic forgetting, and the spacing effect (Knoblauch et al., 2014;Knoblauch and Sommer, 2016). These results are based on the fact that the number of potential synapses a neuron could develop, i.e., its potential connectivity, is much greater than the actual number of synapses, and structural plasticity allows the system to explore different wiring possibilities (Stepanyants et al., 2002;Fares and Stepanyants, 2009).
Here, we use an adaptive-sometimes also called co-evolvingbrain network model, which has already been used by us to describe synaptic pruning in humans (Millán et al., 2018a,b), to analyze how the dynamical processes of adding and removing synapses during brain development can affect the ability of the network to store and optimally retrieve a given set of memories. Our system combines the auto-associative Amari-Hopfield neural network (Amari, 1972;Hopfield, 1982) with a preferential-attachment dynamics for the network evolution in a way that has been shown to accurately reproduce the observed variation of neuron connectivity data on human brains during infancy (Johnson et al., 2010;Millán et al., 2018b). As empirically observed-see (Holtmaat and Svoboda, 2009) and references therein-this model assumes that the probabilities of growth and death of synapses depend on both the mean connectivity in the system and the neural activity. Previous studies have analyzed the effect of thermal noise in the system and its emergent behavior, and they have shown that the coupling between neuronal activity and connectivity creates a feedback loop between form and function since the system activity influences its topology and, in turn, it is affected by the network structure through the synaptic currents the neurons receive (Millán et al., 2018a). As a matter of fact, depending on parameters, this system is then able to produce heterogeneous networks with the presence of hubs, similar to the ones observed in actual neural systems (Van Den Heuvel and Sporns, 2011;Crossley et al., 2014;Oh et al., 2014;Stafford et al., 2014), with high memory retrieval and noise tolerance. Another recent work has also studied the effect of a transient period of high connectivity before synaptic pruning begins (Millán et al., 2018b), as observed in mammal brains (Huttenlocher and Dabholkar, 1997;Navlakha et al., 2015), demonstrating that it has beneficial effects for memory recovery and the emergence of an organized stationary state in the system.
Here we develop on the effect that synaptic (or quenched) disorder resulting from the interference among many patterns of activity-stored by Hebbian learning on the synaptic weightshas on the emergent behavior of the system. We show that, as a consequence of the interplay between structural (i.e., pruning), thermal and quenched disorder, oscillations can emerge in the activity of the model which imply visiting different memorized patterns, an emergent behavior that had not been reported before in this model. This intriguing behavior is precisely due to long-term synaptic mechanisms associated with the network evolution dynamics, and not to short-term synaptic processes, such as synaptic depression and facilitation (Pantic et al., 2002;Marro et al., 2007;Torres et al., 2007Torres et al., , 2008Torres and Marro, 2015) or spike frequency adaptation (Knoblauch and Palm, 2002;Ha and Cheong, 2017), which are not present in our model. These have already been described to induce oscillations among stored patterns of network activity, however the biophysical mechanisms behind them are different from the topological rewiring process considered here, and in particular they act on shorter time-scales-on the order of ms as opposed to the time scale of hours or days in which synaptic rewiring operates. It would be straightforward to extend the present study by adding short-term mechanisms, and we hypothesize that the interplay between different neuron and synaptic processes during learning and brain evolution could give rise to other types of oscillatory phenomena associated with non-equilibrium phases not yet reported, a fact that we glimpse could have strong computational implications.

MODEL AND METHODS
Our system consists in a time-dependent, symmetric, undirected, N-node complex network (Boccaletti et al., 2006) of neurons, defined at time t by the adjacency matrix e e e(t), with elements e ij (t) = {0, 1}, in which each node represents a neuron and each edge [e ij (t) = 1] stands for a synapse. The degree of node i at time t is defined as and the mean degree of the network is Following a familiar (Amari-Hopfield) prescription (Amari, 1972;Amit, 1989), each neuron i is modeled as a stochastic binary unit, s i (t) = {0, 1} (representing respectively a silent and a firing neuron), whose state evolves in time according to the probabilistic dynamics is the local field at neuron i quantifying the incoming input from neighbor neurons and is the neuron's threshold for firing. This definition of the threshold is typically considered, in the case of static networks, when the more biologically plausible {0, 1} code is used instead of the canonical {±1} one, since it allows one to recover the phase diagram of the canonical, fully connected Amari-Hopfield model (Amit, 1989). Therefore, we maintain it when extending the model to a time dependent topology, and it naturally leads to a dynamic threshold. This is not a strong assumption since dynamic or adaptive thresholds have been widely described in several neural systems. For instance, they have been shown to create a nontrivial motion between the attractors of the system (Horn and Usher, 1989;Itskov et al., 2011) and to have a major role in stochastic resonance (Mejias and Torres, 2011) and in the functioning of sensory systems (Fricker et al., 1999;Gray, 2000, 2003;Cardin et al., 2008;Kobayashi et al., 2009). Mechanisms of threshold adaptation have been found to help to avoid saturating activity during developmental changes (Turrigiano et al., 1998), and to be related to homeostatic regulation mechanisms observed in cortical neurons (Abbott and LeMasson, 1993;Turrigiano et al., 1998), and to the emergence of self-organized criticality in neural systems (Uhlig et al., 2013;Hobbiss et al., 2018). In our context, θ i (t) depends only on the existing synapses, which can be seen as a means of homeostasis since the response of a neuron is regulated by the number and strength of its synaptic contacts, thus avoiding silencing lowdegree neurons and saturation of hubs. Furthermore, in our model the term e ij w ij in Equation (5) characterizes the intensity of the synaptic transmission between neurons i and j, so that the threshold dynamics depends indirectly on the neural activity.
On the other hand, the noise parameter or temperature T (T > 0) sets the level of stochasticity on the activity of the neurons, so that if T = 0 the evolution of the system is deterministic and the state of a neuron at time t is completely determined by the states of its neighbors at time t − 1. For T > 0, however, the evolution is stochastic and, as T is increased, the thermal noise has a stronger effect. The strength of each synapse, or its synaptic weight, w ij , is a real variable defined by means of a set of P binary patterns of neural activity, ξ µ i ∈ {0, 1}, µ = 1, ..., P, according to the Hebbian learning prescription (Amit, 1989), where κ 0 = κ(t = 0) and a 0 is the mean activation of the patterns, i.e., a 0 = (NP) −1 P µ=1 N i=1 ξ µ i . This definition of the synaptic weights makes the patterns ξ µ i attractors of the activity dynamics of the system, and therefore it constitutes the final step of a process of "learning" or "storing" of a set of activity patterns by the system in the synaptic weights. Notice also that w ij = w ji by construction so that the network is symmetric, in the spirit of previous studies (Sompolinsky and Kanter, 1986). This is for simplicity and also as a reference to compare with the canonical Amari-Hopfield model (Amari, 1972;Hopfield, 1982).
The overlap of the network state with each of these patterns determines the global state of the system, and it is defined as It follows from this definition that −1 ≤ m µ (t) ≤ 1. We say that the system is in a memory state or, equivalently, that it has retrieved pattern µ, if m µ > 2/3 and m ν → 0 ∀ν = µ. This indicates that the activity state of the network strongly resembles that of pattern µ. In the case of a non-trivial topology, it is also of interest the degree dependent overlap, m µ (k, t), defined as where p(k, t) is the degree distribution of the network, which indicates the probability that a node has degree k at a certain time t (p(k, t) ≥ 0 ∀k, t, N k=1 p(k, t) = 1 ∀t). Therefore, m µ (t) = N k=1 p(k, t)m µ (k, t). Notice also that if the patterns of activity are not homogeneously distributed through the neurons, m µ (k, t) is not bounded by ±1.
The "canonical" setting of the Amari-Hopfield model, in the case of a fully connected network and random orthogonal patterns, exhibits three characteristic phases. In the absence of thermal noise, T = 0, the patterns ξ µ i are stable attractors of the dynamics of the system for P < P c = 0.138N, and the system is in what is called the memory phase. P c defines the maximum storage capacity of the network (Amit, 1989), that is, the maximum quantity of information-or number of patternsthat can be stored and effectively retrieved from the network. This phase is (mathematically) equivalent to the ferromagnetic or ordered phase of interacting spin networks (as in the Ising model). The storage of a large number of different patterns in the network gives rise to quenched noise as a consequence of the interference between them in w ij , which can destabilize such memory phase. Therefore, if P is further increased above P c there is a discontinuous phase transition to a spin-glass (SG) phase, in which there appear metastable states which are combinations of the stored patterns (therefore called mixed states) that trap the dynamics of the system. Similarly, in the case of P = 1 when T > T c = 1, there is a continuous phase transition from the memory phase to a noisy or paramagnetic phase (also called disordered phase) in which there are no stable attractors, and the dynamics of the system is driven by noise (Amit, 1989). In the more general case in which both T > 0 and P > 1, the location of the phase transitions depends both on T and P. The emergent behavior of the Amari-Hopfield model has also been studied on non-trivial network topologies, such as scale-free and small-world networks (Torres et al., 2004;Boccaletti et al., 2006;Oshima and Odagaki, 2007). Such systems have been shown to present the same phases as the canonical fully connected model, with transition lines that depend on the topology. In particular, it has been reported that, for heterogeneous networks and a single stored pattern, the overlap reduces for T < T c , so that memory is recovered but with more errors than in a fully connected network. However, the critical temperature diverges, T c → ∞ as N → ∞, due to the presence of hubs that retain pattern information. Therefore, the memory phase expands to much higher values of thermal noise. On the other hand, the capacity of the network is known to decrease as the mean connectivity of the network decreases (Torres et al., 2004).
We here consider an evolving network whose structure, contrary to the canonical model above, changes constantly in time subjected to the pruning dynamics, as we shall describe. Moreover, we consider a highly sparse network, with values of κ/N ∈ 10 −3 , 10 −2 , which can be homogeneous (i.e., every node having roughly the same connectivity degree), or heterogeneous, with the formation of hubs. Both sparseness and heterogeneity damage severely the memory retrieval ability of the neural network that, for such cases, diminishes fast with P compared with the case of highly connected and homogeneous neural networks (Stauffer et al., 2003;Castillo et al., 2004;Morelli et al., 2004;Torres et al., 2004;Oshima and Odagaki, 2007;Akam and Kullmann, 2014) However, there is experimental evidence that the configurations of neural activity related to particular memories in the animal brain involve many more silent neurons, ξ µ i = 0, than active ones, ξ µ i = 1 (Chklovskii et al., 2004;Akam and Kullmann, 2014). Notice that in this case there is a positive correlation between different patterns due to the sparseness, since a 0 = 0.5, which is also known to improve the storage capacity of a neural network (Knoblauch et al., 2014;Knoblauch and Sommer, 2016), and in particular that of heterogeneous and sparse neural networks (Morelli et al., 2004). Consequently, we consider here this kind of activity patterns, and we further define them as non-overlapping regions of active neurons, each consisting of N/P neurons, so that they cover the whole network (and therefore the mean activity of the patterns is a 0 = P −1 ). This corresponds to a particular definition of sparse or biased patterns, which in other works have been considered to be randomly distributed with a given a 0 (Knoblauch et al., 2014;Knoblauch and Sommer, 2016), what allows for a good visualization of the activity of the network by means of the raster plots.
Moreover, this scheme allows us to define another measure of the overlap between the state of the system and the memorized patterns, considering only the corresponding active neurons as with m µ 1 ∈ [0, 1]. If is also of interest its binearized extension, m µ B , defined as m µ B (t) ≡ 1 if m µ 1 (t) > m th and 0 otherwise, so that m B (t) = m 1 B (t), m 2 B (t), ..., m P B (t) indicates, in a binary code, which combination of patterns is recovered at time t. Equivalently, the decimal variable d s can be defined, which one can interpret as a one-dimensional variable indicating the global memory state of the system. Interestingly, the activity patterns defined here are such that when a number P r of them are recovered at the same time, in a SG-like state, the maximum overlap that they can have is less than one. In order to see this, one can decompose Equation (7) in P sums, each over the neurons corresponding to the region associated with each of the activity patterns, as where N = a 0 N = N/P is the size of each region and the time dependency has been dropped for clarity. Here, the first sum is over the P patterns stored in the network, whereas the second one goes over the N neurons in the region associated with each pattern. If the pattern µ is recovered together with other P r − 1 patterns, then the sum over ν can be split in three terms: the region associated with the pattern µ, the ones corresponding to the other retrieved patterns, and finally those of the non-retrieved patterns (which do not contribute to the sum). Therefore, the overlap corresponding to this pattern is which only meets the equality in the case P r = 1, that is, if only the pattern µ is retrieved. The network structure changes in time following a preferential attachment process. This is characterized by the probability each node i has to gain or lose an edge at each time t -namely, where I i = h i − θ i is the scaled input that each neuron receives as a consequence of the coupling with its neighbors, a sort of recurrent current in the network, and the time dependence has been dropped for clarity. Here, u and d account for global factors that affect synaptic growth and death, such as the diffusion of different molecules through large areas of tissue, for which the mean degree κ is taken as a proxy. The second terms π and η introduce a dependence on the pre-synaptic activity of the nodes, closing the activity-topology coupling. This creates a feedback loop between the evolution of the structure of the network (form), mediated by the local currents, and the neural activity on the network (function).
Taking the local probabilities to be normalized over the network, the number of edges that are added and removed at each time t depends only on the global probabilities u(κ(t)) and d(κ(t)). In this way, they determine the temporal evolution of the mean connectivity κ(t), whereas the local probabilities π(I i ) and η(I i ) characterize the second order statistics of the network structure, such as the variance of the degree distribution or the degree-degree correlations, as we show below (see also Millán et al., 2018a). These definitions allow us to simulate the dynamics of the system via a Monte Carlo method (in particular, we make use here of the BKL algorithm Bortz et al., 1975) as follows. First, the number of edges to be created and destroyed at time t is sorted according to the global probabilities u(κ(t)) and d(κ(t)). Then, we select as many nodes as indicated with these draws, independently of each other, according to π(I i ) and η(I i ). This process is done in a serial manner, and the same node can be selected more than once. Notice that for each node i that gains or looses an edge e ij , the degree of the second node j to which that edge links also changes accordingly. Therefore, there are in fact two paths that can lead to the change of a node's degree: either through the primary process with probability π(I i ) for a gain (or η(I i ) for a loss), or when it is randomly connected to (or disconnected from) an already chosen node. Therefore, the effective values of the second factors in Equation (13) are where the 1/2 factor is included to assure normalization. Following our previous work, we consider which are normalized over the network, The power-law relation in π i allows us to explore both suband super-linear responses by just modifying a single parameter, namely α. The probability η i , on the other hand, is fixed in a linear response, which corresponds to edges being chosen at random for removal, which can be seen as a first order approximation to the pruning dynamics (Millán et al., 2018a). Therefore, α is the control parameter for the pruning dynamics. If α < 1, high degree nodes are more likely to lose edges than to gain new ones, thus creating a homogeneous network structure. On the other hand, if α > 1, high degree nodes are more likely to continue to gain edges than to lose them, which gives rise to a highly heterogeneous, bimodal structure. Finally, the case α = 1 corresponds to the critical case in which networks develop a scale-free topology as shown in previous works (Johnson et al., 2010) that reproduces the scaling behavior observed in the longrange connections of the human brain (Gastner and Ódor, 2016) and in protein interaction networks (Albert, 2005), which decay as a power-law with exponent µ ≈ 2.5.
The local probabilities are then given by We impose further restrictions on the network. First of all, e ij is a binary matrix, so that only one edge per pair of nodes is allowed and the strength of the connection between two neurons, resembling the number of multiple contacts between actual neurons (Fares and Stepanyants, 2009), is considered to be given by w ij . Moreover, we set the minimum degree of the network, k i = 1, so that there cannot be any disconnected nodes, and we forbid self-connections, e ii = 0 ∀i. The maximum degree a node can have is therefore N − 1. We do not impose a hard bound on it as other works have done (Knoblauch et al., 2014;Knoblauch and Sommer, 2016). This would exclusively affect hubs, which only appear for α > 1, as discussed above (Johnson et al., 2010;Millán et al., 2018a), reducing their connectivity. This might affect the memory capabilities of the network in the limit P → ∞ but, since we do not work on this limit, we do not expect any changes on the qualitative behavior and main findings of our model.
Under this framework, the evolution of the mean degree is For a careful derivation of this equation, we direct the reader to Johnson et al. (2010). Intuitively, u κ(t) and d κ(t) set the number of edges that are created and destroyed at every time step, so that u κ(t) − d κ(t) gives the net change in the number of edges. Since for each of these edges two nodes change their degree, there is a factor 2 in the variation of the mean degree. The simplest way to approximate the pruning dynamics is to consider an exponential decay of κ(t) from κ 0 to κ ∞ , where κ 0 = κ(t = 0) is the initial mean degree of the network and κ ∞ = κ(t → ∞) the stationary mean degree after synaptic pruning has occurred, so that κ 0 ≥ κ ∞ . This is achieved by defining where the parameter n sets the timescale for the pruning dynamics. Notice also that Equation (18) assures that u(κ) ≥ 0 ∀κ. By substituting these definitions into Equation (17), we obtain the time evolution of κ(t), where τ p = Nκ ∞ /(2n). This set-up has been previously used to reproduce experimental data on the connectivity of the human pre-frontal cortex using values of κ 0 ∈ (60, 80) and κ ∞ ∈ (30, 50), depending on the region, and also of the mouse somatosensory cortex, with κ 0 = 3.10 and κ ∞ = 1.64 (for the other parameters see Johnson et al., 2010;Millán et al., 2018a). The definitions in Equation (18) take into account that synaptic growth and death relay in some way on the concentrations of various molecules (that can have an important role in synaptogenesis, as axonal growth factors), which can diffuse through large areas of tissue and therefore cannot in general be considered local (Klintsova and Greenough, 1999), and here we consider κ(t) as a proxy for the amount of resources consumed by the existing synapses in the network. In an environment with a finite presence of nutrients, it is reasonable to think that there is a competition for the existing resources, and that neurons are sensitive to the amount of nutrients available to them, so that synapses are less likely to grow, and more likely to atrophy, when the connectivity is high, and viceversa, as assumed by Equation (18). Finally, the network macroscopic state is described via the degree distribution p(k, t) and its homogeneity, defined as where σ 2 (t) = k 2 (t) − κ 2 (t). For homogeneous networks, in which all nodes have similar degrees, σ (t) is small and the homogeneity approaches one, with the trivial case of g = 1 if p(k, t) = δ k,k 1 . For heterogeneous networks, on the other hand, there are big fluctuations in the degrees of the nodes and g(t) → 0.
The timescale for structure changes is set by the parameter n, whereas the time unit for activity changes, h s , is the number of Monte Carlo Steps (MCS) that the states of all neurons are updated according to the Amari-Hopfield dynamics between each structural network update. Our studies show a low dependence on these parameters in the cases of interest, so we only report here results for h s = 10 MCS and n = 10. Measures on the stationary state of the system are carried out by temporal averages of the macroscopic variables,f = A recent work (Millán et al., 2018a) showed that, within this framework, three phases emerge: a homogeneous memory phase when both α and T are low (T, α < 1), in which the network is capable of memory retrieval and the topology dynamics keeps a homogeneous configuration; a heterogeneous memory phase for high α (α > 1) in which the dynamics leads to bimodal networks (with the appearance of hubs or highly connected nodes); and a homogeneous noisy phase for high noise T. However, as we will depict in the next section, the combination of thermal noise together with the introduction of a larger number of patterns of activity-which induces interference among them-induces other non-reported non-equilibrium phases characterized by the emergence of complex oscillations among the activity associated with the stored patterns.

RESULTS
Previous preliminary analysis of the storage capacity of developing brains under the present framework revealed that the capacity of the network can be greatly improved if a feedback loop between structure and function is considered (Millán et al., 2018a). This is because the interplay between form and function gives rise to a topological structure that enhances the stability of the memory attractors which are recovered during the evolution of the system. In order to explore this interesting picture under other conditions, here we analyze in detail the phase diagram of the system with respect to four relevant parameters in the model, namely, α, κ ∞ , T, and P. The first two characterize the network structure dynamics, whereas the temperature, T, and the number of stored patterns, P, account respectively for thermal and quenched disorder. As already said, the latter is a consequence of the interference among many stored patterns, and it can affect the recall process. Other parameters, such as the initial connectivity κ 0 or the speed of the pruning, n, where shown to have little or no effect on the dynamics (Millán et al., 2018a,b).

Steady State Solutions for T = 0
We first analyze the behavior of the system at T = 0, that is, in absence of thermal fluctuations that can affect the stability of the fixed point solutions of the system dynamics. As stated above, there are, however, other sources of noise in our system which can have a prominent influence in its behavior. One is the interference among stored patterns, which can significantly reduce the memory retrieval ability of the system (Amit, 1989). Another is the pruning dynamics that adds a second source of noise; this is an intrinsic, structural noise that emerges due to the stochastic adding and removal of synapses associated with the network dynamics during brain development, and which can dynamically affect the performance of the system during memory acquisition and consolidation.
In Figure 1 we show the corresponding phase diagrams of the system (depicting different phases or kinds of behavior) for different values of κ ∞ = 20, 40, and 60, respectively from left to right. These depict some non-equilibrium phases associated with different computational abilities during memory recall. The top panels show, in the steady state, the ratio of patterns that can be retrieved with high overlap (m µ ≥ 0.66), namely g P ≡ P r /P (where P r is the number of retrieved patterns), as a function of α and P. A value g P = 1/P indicates a pure memory state, whereas larger values correspond to mixtures and SG-like states (Amit, 1989), and g P = 0 corresponds to the noisy or non-memory state. Meanwhile, the middle panels show the mean overlap of the recovered patterns during memory recall, namely m P and, finally, the bottom panels show the stationary mean homogeneity,ḡ.
These diagrams show up different types of dynamical behavior. In order to illustrate the characteristics of each one, in Figure 2 we depict the time series m µ (t) (top graph of each panel), raster plots showing the whole activity of the system (bottom graph of each panel) and the steady-state degree distribution (inset of each panel) for some particular values of α and P corresponding to different characteristic behaviors in the phase diagrams in Figure 1. For a given stationary connectivity (e.g., κ ∞ = 60, Figures 1C,F,I) we find that, for P = 1, the system is able to retain memory for almost every value of α, as it can be seen by the yellow region at P = 1 in Figure 1F, which indicates an overlap equal to 1. For small P and small α, SG-like states, or mixture states (in which some of the memories are partially FIGURE 1 | Phase diagrams depicting the steady-state of the system with respect to P and α at T = 0. The first row of panels (A-C) shows the fraction of patterns, g P , retrieved after a given transient; the second one (D-F) shows the average overlap with the recovered patterns, m P ; and finally the third one (G-I) shows the stationary homogeneityḡ. Each column is for a different value of κ ∞ , with κ ∞ = 20, 40, and 60, respectively from left to right. The network size was set here to N = 1, 600 and each point has been averaged over 10 realizations of the system. In this figure a memory phase appears as a blue region in the diagram of g p and a high value of m p , indicated by a yellow color. A SG phase appears as a blue region in g p and a lower value of m p , indicated by a green color, whereas a noisy phase appears as black in g p and m p . Similarly, homogeneous structures take place for high values ofḡ, indicated by a yellow region in the corresponding diagram, whereas heterogeneous structures are for low values ofḡ, indicated by a black region. Finally, the pink stars in diagrams (C,F,I) indicate the (P, α) points corresponding to the time series shown in Figure 2. retrieved at the same time), start to emerge as it is illustrated in Figure 2A, which corresponds to the point α = 0.5, P = 10. As a consequence, both g P and m P take intermediate values; the former since only a finite number of patterns is retrieved, g p < 1 (light blue region of the diagram in Figure 1C), the later because these are retrieved at the same time, and therefore the overlap is reduced, m p < 1 (green and light-blue region of the diagram in Figure 1F). In general, however, the observed SG-like states present high values of the overlap with all the recovered patterns due to the high correlation between memories we have considered in this work. Moreover, in this region the network structure is homogeneous since α < 1, so thatḡ approaches 1 and the degree distribution resembles a Poisson distribution (see Figure 1I and the inset of Figure 2A). In these conditions, when P is increased the memories lose stability until there is a transition from the SG-like state to the noisy one, where the network structure remains homogeneous, as shown in Figure 2B for the point α = 0.5 and P = 30. This is indicated by g P → 0 (black region in Figure 1C), m P → 0 (black region in Figure 1F), andḡ → 1 (yellow region on the bottom-right side of Figure 1I).
On the other hand, for high α (α > 1), just one (or very few) pattern is retrieved, with m P ≈ 1, and the network structure becomes heterogeneous since α > 1 (see inset of Figure 2C). As a consequence, g P → 1/P, (dark-blue region in Figure 1C), m P approaches 1 (yellow region in Figure 1F), andḡ → 0 (black region in Figure 1I). Memory is achieved due to heterogeneity and the presence of hubs, which can maintain the information content of the retrieved pattern even in the presence of the strong noise induced by the interference with other stored patterns and the dynamic changes of the network structure. Therefore, when P is increased the recovered patterns remain stable, so that m P remains close to 1 ( Figure 1F) and g P decreases as 1/P since only one pattern is retrieved (Figure 1C, see also the inset of Figure 2D, showing the appearance of hubs, and Figure 1I, The stationary mean connectivity of the network, κ ∞ , also affects the behavior of the system, as it determines location of the phase transition from the SG phase to the noisy one for α < 1. As the diagrams in Figure 1 show, larger values of κ ∞ increase the tolerance of the system to quenched disorder, so that a bigger number of patterns can be stored. This is in line with the known result that the information is stored in the synaptic weights, and therefore increasing the number of synapses also increases the amount of information that the system can store (Amit, 1989).
Notice also that the qualitative state of the system is approximately independent of P for P > 20, as shown in Figure 1 where one can see that g p , m p andḡ remain essentially constant as P is increased with constant α above P = 20, in agreement with previous studies (Millán et al., 2018a). Therefore, in the following we restrict our analysis to the most interesting region P < 20 and do not analyze the large storage limit of the system (Knoblauch et al., 2014;Knoblauch and Sommer, 2016). This is because our interest here is in characterizing the dynamic behavior arising as a consequence of the interplay between structure and function FIGURE 2 | Time evolution of the system at T = 0 and κ ∞ = 60 in four typical cases corresponding to different values of α and P as marked with pink stars in the phase diagrams of Figure 1. Each composite panel illustrates the overlap time series, m µ (t), (top graph), raster plots of neuron activity (bottom graph) and the steady-state degree distribution of the network (inset), computed at t = 10 6 Monte Carlo Steps and averaged over 10 realizations of the system. The panels correspond respectively to α = 0.5 and P = 10 (A), α = 0.5 and P = 30 (B), α = 1.5 and P = 10 (C), and α = 1.5 and P = 30 (D). In all presented simulations we set N = 1, 600.
under the presence of thermal and quenched noise, rather than its storage capacity. Similarly, the inclusion of a hard bound on the maximum degree of the nodes would primarily affect the degrees of the hubs of bimodal networks. However, these typically form a highly connected core in the network, so the average path length between nodes would not increase heavily. Therefore, we expect that this bound would not have an important effect in the regime P ≪ N in which we set the system here.
In summary, for T = 0, that is, when there are only two sources of noise in the system (structural and quenched disorder), the stationary state for a given P depends strongly on the network structure, determined by α. As so, for α > 1, the network develops heterogeneous structures in which hubs arise. These are very densely connected with the rest of the network, and can maintain information about the memories even when P is very high. For α < 1, on the other hand, the network is always homogeneous, with every node having similar, low degree, and a SG-like phase soon arises, which is then suddenly lost as the quenched disorder becomes too strong and finally the system falls into the noisy state.

Behavior of the System for T > 0
Our previous analysis has determined the phase diagram of the system at T = 0, which characterizes the effect of the dynamical topological structure on the memory capabilities of the system. In this section, we consider the effect of thermal noise in our system's emergent behavior. In order to do so, we analyze in Figure 3 the phase diagrams of the system with respect to α and T, for some representative numbers of stored patterns, namely P = 5, 10, 15, and 20 (each column of the figure corresponds to a different P), and for three values of κ ∞ = 20, 40, and 60, as before. The selected values of P correspond to the left region of diagrams in Figure 1, where the phase transitions from memory to the SG and noisy states takes place. In order to illustrate better the behavior of the system in the cases of interest, we also include in Figure 4 the time series of m µ (t) is some exemplary  Figure 4. Results are for N = 1, 600 and have been averaged over 5 realizations of the system dynamics. In this figure a memory phase appears as a blue region in the diagram of g p and a high value of m p , indicated by a yellow or green color. A SG phase appears as an orange region in g p and a lower value of m p , indicated by a green or blue color, whereas a noisy phase appears as black in g p and dark-blue in m p . Finally, the oscillatory phase appears for high values of g p , (light yellow regions in the corresponding diagrams) and relatively low values of m p (associated blue regions of the corresponding diagrams). Similarly, homogeneous structures take place for high values ofḡ, indicated by a yellow region in the corresponding diagram, whereas heterogeneous structures are for low values ofḡ, indicated by a black or dark blue region.
points for κ ∞ = 20, as indicated in the phase diagrams by a pink star. Each panel corresponds to a given value of P, and each graph on them to a point in the (T, α) space. We find that the combination of thermal and quenched disorder, associated with the interference among patterns, can give rise to oscillations among the memorized patterns for α < 1-that is, when the networks are homogeneous-and T < 1, which are correspond to the yellow regions in the g p panels of Figure 3. Note that the observed oscillations occur at level of the neuronal population as measured by the global network parameter m µ (t), and not on the single neuron level-which appear as small, high-frequency oscillations of m µ (t).
In order to illustrate the emergent behavior of the system, we refer here to Figure 3A, which corresponds to κ ∞ = 20 and P = 5. The top graph in the panel represents g P , the number of patterns visited by the system after the transient  1.1, 1.1). We have selected slightly different points for each P so as to show an example of the oscillatory behavior in each case, and the region of its appearance depends on P. Results are for N = 1, 600. evolution takes place. For α > 1, we find that this number remains finite, and greater than zero, up to very high values of the temperature (T ≈ 2.0) corresponding to the lightblue region in the top panel of Figure 3A, for instance. This indicates that the system is in a memory state (or in a SG-like state in which only a small number of patterns are retrieved), such as the ones depicted in the bottom graphs of Figure 4A. The stability of the memory state for T > 1 is possible due to the emergence of heterogeneous structures (since α > 1), and consequently hubs, as indicated in the black region of the bottom plot of Figure 3A, sinceḡ → 0 for α > 1. Notice also in the middle graph of Figure 3A how, as T is increased, the overlap corresponding to these states, m P , decreases, indicating that these states are also becoming less stable as the thermal noise becomes stronger. In these conditions, only the more densely connected hub nodes are able to maintain information about the memories, and these are the ones contributing the most to the overlap.
As α is decreased, however, the behavior shown in the diagrams becomes more complex and different regions (phases) start to emerge. We find, as expected, that memory is completely lost for T ≫ 1, i.e., due to the strong noise the system falls into the noisy or non-memory state (as indicated respectively by the black region and by the dark-blue region of the top and middle diagrams of Figure 3A). In fact, now networks are homogeneous (α < 1) and there are no hubs that preserve memory (as indicated byḡ → 1 in the bottom plot of Figure 3A, indicating that the degree distribution is homogeneous). A typical time series of m µ (t) for this situation is shown in the top-right graph of Figure 4B.
For small values of T and α (T, α < 1) on the other hand, g P → 1 (orange and yellow region of the top panel of Figure 3A), indicating that a great number of patterns are being retrieved (P r → P) with a moderate value of the overlap m P , as indicated by the green and blue region in the middle panel of Figure 3A. Moreover, results in Figure 2 indicate that, at least in some cases, g P actually increases when T goes from 0 to 1 (see, e.g., the yellow area in the top panel of Figure 3A). That is, as the temperature increases, more memories take place in the state of the system, with a relatively high overlap m P . This is because, for α < 1, there is a wide region of oscillatory behavior between the SG-like and the noisy phases corresponding to the yellow region of the g p diagrams in Figure 3. An exemplary series of oscillations is illustrated in the top-left graph of Figure 4A. This emerges as a consequence of the interplay between structural and thermal noise and the activity of network, since the process of addition and removal of synapses, that creates a dynamical network structure, together with the thermal noise, can make the recovered patterns unstable. Moreover, given that α < 1, the structure of the networks remains homogeneous and no real hubs emerge. Notice however that due to the non-trivial interplay between activity and topology, in the region of oscillatory behavior the networks display a more heterogeneous structure, andḡ < 1. This effect will be discussed in more detail in the following section.
As in the previous section, we have also analyzed the role of the stationary mean connectivity, κ ∞ , on the phase diagram of the system for T > 0 (see Figure 3). This parameter holds physiological interests since it can be taken as a measure of the extension of the process of synaptic pruning. From this point of view, a brain that has undergone a more drastic synaptic pruning would have smaller κ ∞ than one that has been less pruned. This is to be related to recent experiments that have associated an excessive pruning in certain brain areas with schizophrenia (Sekar et al., 2016), whereas ASD has been related to a defect of synaptic pruning (Tang et al., 2014). We find in our model that the area associated with the oscillatory behavior (for α, T < 1) for a given κ ∞ is maximum at intermediate values of P: for very small P there is a dominance of stable SG-like states, whereas for large P the system falls easily on the noisy phase. Similarly, for a given P the greater extension of the oscillatory phase is found for an intermediate κ ∞ . For instance, for P = 10 the noisy phase extends to T < 1 for κ ∞ = 20 ( Figure 3B) and the oscillatory region is small, whereas for κ ∞ = 60 ( Figure 3J) there is a combination of stable SG-like states and oscillations for α, T < 1. Finally for κ ∞ = 40 ( Figure 3F) the oscillatory phase is most robust. Consequently, the absence of dynamical memories in the system could be associated with a defect of the pruning process that causes κ ∞ to be greater than usual, and could be therefore associated with ASD. Interestingly, it has been recently reported that short-term memory and episodic memory are impaired in ASD subjects (Poirier et al., 2011;Lind et al., 2014), which is consistent with our findings here since, in order to be able to recall a sequence of memories, it is first necessary to destabilize the already recalled ones so as to allow the system to remember new ones. On the other hand, schizophrenia is typically associated with erratic behavior (Loh et al., 2007), which could be related to the high frequency memory oscillations found here for smaller values of κ ∞ .

Emergence of Hubs
The appearance of hubs and heterogeneity plays a significant role in the emergent dynamics of the system. In particular, with a given level of noise (T > 0), the topological structure of the network determines whether the system relaxes to a memory state, wanders among different patterns or falls into a noisy state. Therefore, here we discuss in more detail the emergence of hubs during the network evolution and their effect on the emergent state of the system. We first notice that, according to the previous analysis, for α < 1 networks are homogeneous, as evidenced by the homogeneous degree distributions shown in the insets of Figures 2A,B. This is also revealed by the high value of the homogeneity parameterḡ shown in Figures 1, 3 for α < 1, indicating that the variance of k i is small. As a consequence, no real hubs can be defined, since all nodes have similar low degree (given that κ ∞ ≪ N, so that the connectivity of the nodes is bounded). On the contrary, for α > 1 and in the case of memory, networks are heterogeneous as evidenced bȳ g → 0 (black regions in the corresponding diagramsḡ(α, P) and g(α, T) respectively in Figures 1, 3). This indicates that there are nodes with very different degrees and, in particular, the degree distribution p(k, t) is bimodal and it splits in two, as shown in the insets of Figure 2C,D, with the emergence of hubs. Therefore, one can set the connectivity threshold k th -that defines the minimum node's degree to characterize it as a hub-at the value of k at which p(k, t → ∞) presents a local minimum between the two modes. This establishes a clear separation between high and low degree nodes. In particular, in all cases studied here, we find that a threshold k th = 2κ ∞ also suffices to differentiate between homogeneous and heterogeneous structures, since for α < 1 (homogeneous case) the maximum degree of a network is always below 2κ ∞ .
Interestingly, due to the underlying stochastic rewiring process and to the system's finite size, there is always some variability in the degrees of the nodes and, particularly in the region of oscillatory behavior, there is a relative increase in the variability of k i with respect to the SG phase (as evidence by a decreasedḡ in the corresponding diagrams of Figure 3). We argue that this is due to the intrinsic coupling between activity and topology, and to the combination of thermal (since T > 0), topological (due to the ongoing rewiring process) and quenched (due to the learning of different patterns) disorder in the system. In the region of oscillatory behavior, the instability of the memories influences the synaptic currents I i creating variability, thus causing the observed increased heterogeneity. This causes the emergence of relatively-high degree nodes that correspond to the tail of the homogeneous distribution p(k, t → ∞) and which might have an important effect on the system. Therefore, in order to explore as well the dynamics of these relatively-high degree nodes, we have selected a lower threshold, k th = 1.75κ ∞ , for the analysis of hub dynamics.
Hubs (and relatively-high degree nodes for the homogeneous case) dynamics is investigated in Figure 5, where we compare two different cases for P = 5 and κ ∞ = 20. The first one, shown in Figures 5A,B, corresponds to the region of oscillatory behavior for homogeneous networks (α < 1) and it is for T = 0.7 and α = 0.3, corresponding to the bottom-left graph of Figure 4A. The second one, shown in Figures 4C,D, corresponds to the heterogeneous-memory phase for α > 1, and it is for T = 0.3 and α = 1.1, corresponding to the top-left graph of Figure 4A. In particular, we analyze in Figures 4A,C the temporal evolution of the system as given by the overlap m µ (t), the homogeneity g(t) and the hub raster plots, where we represent the existing hubs at each time t, in different colors according to their active or inactive state (respectively pink and green). Furthermore, in Figures 4B,D we show the degree-dependent overlap m µ (k, t 0 ) [defined in Equation (8)] and the degree histogram N(k, t 0 ) = Np(k, t 0 ), for a particular time, t 0 = 5 · 10 6 MCS, corresponding to the systems respectively depicted in Figures 4A,C. We observe, for α > 1 (Figure 5C), that a great number of hubs emerge in the system, and that almost all hubs correspond to the active nodes of the retrieved pattern. Moreover, in this case m µ (k, t 0 ) of the recovered pattern µ is larger for highdegree nodes (Figure 5D), indicating that they contribute most to the overlap m µ (t 0 ). On the contrary, for the non-recovered patterns ν, m ν (k, t 0 ) remains small for all k. On the other hand, for α < 1, no real hubs emerge and only transient relativelyhigh degree nodes are observed in Figure 5A. These do not FIGURE 5 | Emergence and effect of hubs in the system. (A,C) Show the temporal evolution of the system in two representative cases of the dynamics for P = 5 and κ ∞ = 20, corresponding to the emergent oscillatory behavior for T = 0.7 and α = 0.3 (A) and to the heterogeneous memory phase for T = 1.1 and α = 1.1 (C). In these panels, the top plots represent m µ (t), the middle ones the homogeneity parameter g(t), and finally the bottom ones show the existing hubs in the network at each time t, where active and inactive hubs are plotted in different colors (respectively, pink and green). (B,D) Show a snap-shot of the state of the system in (A,C), respectively, at time t 0 = 5 · 10 5 , as represented by the degree-dependent overlap m µ (k, t 0 ) and the number of nodes with degree k, N(k, t 0 ) = Np(k, t 0 ). Results are for N = 1, 600.
only correspond to the recovered patterns but are scattered throughout the network, and no significant correlation can be measured between the pattern oscillations and the hubs dynamics. This causes instabilities that ultimately lead to the oscillatory behavior (see Figure 5B, indicating that relatively high-degree nodes contribute more to m µ (k, t) of the recovered patterns but not only).
In summary, Figure 5 shows that for α > 1 there are active hubs in the system that correspond to the recovered pattern, making it stable. On the other hand, for α < 1 no real hubs can emerge in the system, and the transient relativelyhigh degree nodes are scattered throughout the network, not only corresponding to the recovered pattern, thus inducing the observed oscillatory behavior.

Analysis of the Oscillatory Behavior
In the previous sections we have shown the emergence of oscillations for α < 1 and T > 0 and their relation to the existence of transient relatively-high degree nodes on the network. Here, we develop further on the structure and patterns of these oscillations. For simplicity, we focus on the case of κ ∞ = 20 and P = 5 as before, and in Figure 6 we show a long time series corresponding to this oscillatory phase (T = 0.7 and α = 0.3 as in the top graph of Figure 4A and in Figure 5A). Plots of the active-overlap parameter m µ 1 (t) (Figure 6A) defined in Equation (9), its binearized version m µ B (t) (Figure 6B) and the global memory state parameter d s (t) defined in Equation (10) (Figure 6C) indicate that the state of the system corresponds to oscillations between SG-like states in which either 2 or 3 patterns are transiently retrieved. These plots also evidence that the oscillations do not follow any clear periodic or regular pattern.
In order to analyze the pattern of oscillations, we show the power spectra of m µ B (t), S µ (f ), and of d s (t), S s (f ), in Figure 6F, which displays a power-law decay with an exponent equal to −0.9, indicating that there is not a dominant frequency of the oscillations, but that jumps between different patterns occur at all time scales. This is in accordance with previous studies that have repeatedly reported 1/f -type noise in brain activity under healthy conditions. It has been reported, for instance, in electroencephalogram (EEG) and functional magnetic resonance (fMRI) measures of human brain activity (Linkenkaer-Hansen et al., 2001;Voytek et al., 2015) and also in behavioral processes related to human cognition and motion as well as animal motion (Chialvo, 2010). 1/f noise indicates the existence of temporal correlations within the data, and has been related to emerging self-organized criticality in the brain (Chialvo, 2010).
Moreover, we also investigate the frequency of appearance of each global state, as seen in Figure 6E, which indicates that global states have different probabilities of occurrence in each realization of the system. However, when averaged over realizations, the mean probability of each state,p s , converges to a uniform distributionp s → 1/N s , where N s = 20 is the total number of possible 2− and 3−pattern states (p s = 0.054(3) when averaged over 20 realizations). Similarly, individual patterns may have different probability of appearance, p µ , in each realization of the dynamics, but when averaged over realizationsp µ converges to 1/P [p µ = 0.21 (1)]. Finally, we also computed the transition matrix between global SG-like states ( Figure 6F) which indicates that, in a given realization of the system, some transitions are preferred by the system depending on the emergent coupling between activity and topology.
In summary, these results show that the oscillations are not periodic, but occur at all time scales, and that all SGlike states are visited in a non-periodic order. Interestingly however, in a given realization of the system not all transitions are allowed, but only some of them occur. It could be interesting to analyze in more detail in further studies whether the coupling between structure and activity induces a particular pattern of oscillations, and how the scaling of the frequency of oscillations depends on the parameters of model.

DISCUSSION
We report here on recent studies of the emergent behavior of developing brain models in which structure and function cooperate and influence each other through a feedback loop, thus affecting the system's memory storage and retrieval abilities. This is a prominent example of how inter-synaptic factors at the network level can affect the processing of information in developing brains in a nontrivial way. In particular, our study focuses on the analysis of the conditions under which such feedback loop can enhance the storage and retrieval of a set of correlated patterns. Our work also pays attention to the emerging dynamics of the system, which is a consequence of the interplay between structural, quenched, and thermal disorder during its maturation. The results presented here demonstrate that a heterogeneous network can greatly improve the stability of the memory patterns, since its structure is optimized to preserve information about them in the network hubs which, as we have shown, correspond to the active neurons of the retrieved memories during the recall process. Moreover, due to the structural plasticity, once a pattern is retrieved, the ability of the system to recall it again increases thanks to pruning optimization. This illustrates the constructive role of synaptic pruning to consolidate memories in the "memory phase" of the system.
Our study also shows that the interplay between thermal noise, the interference among stored patterns and the dynamics driving the evolution of the topology creates instabilities on the memory attractors, which can make the system wander among different configurations for certain values of the parameters. In our model, these oscillations among stored patterns are caused by the topological synaptic plasticity due to the death and birth of synapses, which change the energy landscape of the system. In fact, in the absence of this rewiring process, this oscillatory phase is not present and the model would reduce to an Amari-Hopfield model on top of a non-trivial fixed topology. This has been shown to present the same phases as the canonical fully connected version of the model, with transition lines that depend on the topology, so that for instance the critical temperature diverges, T c → ∞ as N → ∞, due to the presence of hubs that retain pattern information. Interestingly, in our model the oscillatory behavior that takes place on the homogeneous networks phase of the system is also associated with an increased transient heterogeneity of the underlying structure, in which transient relatively-high degree nodes (whose degree is however smaller than typical hubs) emerge and disappear in time. Moreover, these relatively-high degree nodes do not correspond in general to the active nodes of the transiently recovered patterns, but appear distributed throughout the whole network, corresponding also to active nodes of the rest of non retrieved patterns. This creates a non-trivial time-dependent competition among the different patterns which, together with the subsequent removal of some synapses during brain development, can make the currently recalled attractor less stable, thus inducing the observed wandering among the memories.
We have also analyzed the characteristics of the oscillations and shown that the oscillatory pattern is not periodic but presents a power spectrum following a power law scaling decay with an exponent of −0.9, so there are not any preferred frequencies. This in accordance with previous studies repeatedly reporting 1/f noise in brain activity under healthy conditions (Linkenkaer-Hansen et al., 2001;Voytek et al., 2015) and also in behavioral processes (Chialvo, 2010), which is related to the existence of temporal correlations within the data, and has been related to emerging self-organized criticality (SOC) in the brain (Chialvo, 2010).
Interestingly, the appearance of an oscillatory phase characterized by dynamical memories could be useful to enhance the learning and recalling of sequences of patterns of activity, as in episodic memories, without the necessity of any external input or current forcing the retrieval of the memories in the sequence. This type of oscillations has already been reported for brain models with synapses enduring short-term synaptic plasticity (STP) (Pantic et al., 2002;Cortés et al., 2006;Marro et al., 2007;Torres et al., 2007Torres et al., , 2008. This occurs at the synapse level and depends on the activity of the pre-synaptic neuron (which therefore depends closely on the synaptic current I i used in our model Amit, 1989). However, STP is caused by biophysical mechanisms controlling the release and recycling of neurotransmitters at the synapses during synaptic transmission and operates at short time scales of the order of ms (Tsodyks and Markram, 1997). The activity dependent topological plasticity reported here, however, is the result of the interplay between form and function in a developing brain, and the ongoing synaptic rewiring in mature brains, which happens at the time scale of hours or days (Holtmaat and Svoboda, 2009). Moreover, topological plasticity allows the system to explore more efficiently its dynamical phase space and it has been shown to improve the capacity of neural networks by allowing them to organize in a more efficient structure. Both mechanisms could happen at the same time in actual systems, together with neuron level phenomena such as spike adaptation (Knoblauch and Palm, 2002;Ha and Cheong, 2017). We hypothesize that the combination of these mechanisms could lead to the extension of the oscillatory behavior to other regions of the phase diagram, although results would strongly depend on the relative time scale between structural plasticity and STP, and it could be an interesting approach for future works.
It is also worth noting that the reported oscillations in our system are for the overlap function that is a measure of the activity of the whole neuron population during memory recall processes. These occur in actual neural systems at a long time scale-normally days or even years-as it is the case in our model. Temporal changes at the single neuron level appear in our system as high frequency fluctuations in the time dependent value of the overlap parameter. If the level of stochasticity is low (low T) and the network size is large enough (N ≫ 1), such single neuron fluctuations are very unlike to be significant on m µ (t). In any case, the model output could be easily tuned up to obtain faster of slower oscillations in the overlap function to match more realistically actual experiments during learning and recalling. This could be done by varying some model parameters to make the recall process more or less efficient in time, or to allow the system to recall dynamic memories-such as episodic memories-that are learned and recalled at different stimuli input frequencies.
We have analyzed in detail how the dynamical behavior of the system depends on the synaptic factors affecting the addition and removal of synapses and on the number of stored patterns. In particular, the stationary mean connectivity of the network, κ ∞ , has been shown to have a great effect on the emergent behavior of the system. For instance, we have found that the absence of dynamical memories in the system, or the presence of memory oscillations with long periods, is associated with a defect of the pruning process. Similarly, we have shown that high frequency oscillations among patterns and more tendency to noisy behavior occur when there is a pruning excess. In particular, the destabilization of recovered memories is necessary for instance to recall a sequence of memories, each during a short period of time, so as to allow the system to remember new ones. One may argue that the induced instability and the associated oscillatory behavior observed in our system could be positive for information processing, since it would allow neuronal media to explore different memories or attractors, for instance following hetero-clinic orbits, and consequently, to process more complex information, such as spatio-temporal patterns of information (see Rabinovich et al., 2006 and references therein). Such emergent behavior could also be useful to respond more efficiently to changing external stimuli, as in episodic memory tasks, as it has been widely stated in previous works in different neural systems (Cortés et al., 2006;Marro et al., 2008;Torres et al., 2008).
The different emergent behavior in the model with varying connectivity could perhaps be associated with cognitive abilities related to autism spectrum disorders (ASD) and schizophrenia. In the former, a pruning defect has been observed in some brain areas (Tang et al., 2014), and it has recently been reported that short-term memory and episodic memory are impaired in ASD subjects (Poirier et al., 2011;Lind et al., 2014). Our results are consistent with these observations since, if the brain is less pruned, the mean connectivity of the corresponding network is higher, thus making the memory attractors more stable. This implies a lower ability of the brain to remember sequences of patterns as described in episodic memory tasks in ASD patients because it is harder for the brain to forget the already recalled pattern due to its strong stability. Results in our model indicate that a lightly pruned brain could be forced out of the memory phase into an oscillatory regime with an increase in the number of stored patterns (see Figure 3). These observations might provide an interesting insight for experimental psychologists to design a cognitive strategy or therapy to learn and recall sequences of patterns, that might improve the cognitive abilities of patients with ASD. On the other hand, we have demonstrated that high frequency oscillations among patterns occur when there is a pruning excess, and this could be perhaps associated with the erratic behavior observed in schizophrenia (Loh et al., 2007), in which case the brain seems to present some areas with an excess of pruning (Sekar et al., 2016). In this case our results here suggest that a learning therapy based on increasing the number of stored memories would not be useful but may in fact be detrimental, as it would make the memory activity patterns more unstable. A learning therapy that moves the patient brain state near to its stable memory phase, for instance, by stabilizing a few old useful memories, could therefore be more convenient.
Finally, we note also that some drastic assumptions have been made in order to simplify the relevant scenario. Firstly, our study is for sparse correlated patterns, as suggested by experimental studies (Chklovskii et al., 2004;Akam and Kullmann, 2014), which are also known to improve the memory retrieval capabilities of the network (Knoblauch et al., 2014;Knoblauch and Sommer, 2016) and particularly so in the case of highly sparse and heterogeneous networks (Morelli et al., 2004). Moreover, we have selected the patterns of activity to be non-overlapping regions of activity, following previous works (Torres and Marro, 2015). This set up corresponds to a particular case that allows for a better visualization of the network dynamics and that has proven out to be useful to investigate the interplay between structure and dynamics, i.e., between form and function, together with the presence of thermal and quenched disorder, on a developing neural network. Similarly, results are for the low storage regime of the neural network, P ≪ N, what allows us to study in detail the dynamical behavior of the system that gives rise to memory wandering. However, given that our qualitative results depend little on P for P > 20, we expect them to hold when P is increased.
Further extensions of this work could also include the consideration of different details of the synaptic pruning process, including for instance the growth of synapses taking place after birth (Millán et al., 2018b), multiple synaptic contacts between neurons (Knoblauch et al., 2014;Knoblauch and Sommer, 2016), or a hard bound on the maximum degree of the nodes (Stepanyants et al., 2002;Fares and Stepanyants, 2009). Moreover, more elaborated definitions of the probabilities of growth and death of synapses (Equations 16 and 18) could also be considered, such as a mechanism of self-organization toward the stationary mean connectivity (Chechik et al., 1999;de Arcangelis et al., 2006;Lewis and Todd, 2007;Tetzlaff et al., 2010) or by explicitly including a dynamics for the available nutrients (Tetzlaff et al., 2010). However these definitions would still need to reproduce the basic characteristics of brain development and synaptic pruning, that is, an initial fast decay of connectivity and an ongoing rewiring of edges after the stationary mean connectivity has been reached. We expect that our main results (existence of a feed-back loop between structure and activity, bistability, and emergence of oscillations) would still hold, at least qualitatively, with these modifications, in accordance with previous studies (Millán et al., 2018b). Similarly, the local probabilities in Equation (16) could also consider more detailed functions, to characterize for instance a specific dependence on the concentration of different proteins and growth factors controlling synaptic growth. These could be obtained experimentally, although to the best of our knowledge it has not yet been done. Our work could thus motivate neurobiologists to design experiments to describe the exact probabilities involved in synaptogenesis and pruning, information that could be easily incorporated in our theoretical framework. Similarly, it could also be particularly interesting to include a learning dynamics that is also coupled to the development of the neural network, thus modeling learning during infancy, or to include an external current on the system, that could certainly be time dependent, to analyze the effect of external inputs on associative memory and memory wandering.

AUTHOR CONTRIBUTIONS
APM, JJT and JM designed the analyses, discussed the results, and wrote the manuscript. APM also wrote the codes and performed the simulations.

FUNDING
We are grateful for financial support from the Spanish Ministry of Science and Technology, and the Agencia Española de Investigación (AEI) under grant FIS2017-84256-P (FEDER funds) and from the Obra Social La Caixa (ID 100010434, with code LCF/BQ/ES15/10360004). This study has been also partially financed by the Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía and European Regional Development Fund (ERDF), with reference SOMM17/6105/UGR.