Elucidating the prefrontal cortical microcircuit has been challenging, given its role in multiple complex behaviors, including working memory, cognitive flexibility, attention, social interaction and emotional regulation. Additionally, previous methodological limitations made it difficult to parse out the contribution of certain neuronal subpopulations in refining cortical representations. However, growing evidence supports a fundamental role of fast-spiking parvalbumin (PV) GABAergic interneurons in regulating pyramidal neuron activity to drive appropriate behavioral responses. Further, their function is heavily diminished in the prefrontal cortex (PFC) in numerous psychiatric diseases, including schizophrenia and autism. Previous research has demonstrated the importance of the optimal balance of excitation and inhibition (E/I) in cortical circuits in maintaining the efficiency of cortical information processing. Although we are still unraveling the mechanisms of information representation in the PFC, the E/I balance seems to be crucial, as pharmacological, chemogenetic and optogenetic approaches for disrupting E/I balance induce impairments in a range of PFC-dependent behaviors. In this review, we will explore two key hypotheses. First, PV interneurons are powerful regulators of E/I balance in the PFC, and help optimize the representation and processing of supramodal information in PFC. Second, diminishing the function of PV interneurons is sufficient to generate an elaborate symptom sequelae corresponding to those observed in a range of psychiatric diseases. Then, using this framework, we will speculate on whether this circuitry could represent a platform for the development of therapeutic interventions in disorders of PFC function.

Keywords: excitation/inhibition balance, GABA, PV interneurons, prefrontal cortex, cognition, psychiatric disorders

# INTRODUCTION

Cognitive impairment stifles independance by making even the simplest everyday tasks seem challenging and burdensome. Through limiting one's ability to focus, encode, retain and manipulate information to make menial decisions, simple endeavors like traveling from one place to another, or completing a typical work assignment become arduous. This decreases productivity, and forces afflicted individuals to rely partially or wholly on the care of others. For most, this loss of independance can be crippling to one's morale, resulting in increasing isolation

#### Edited by:

Vito Di Maio, Istituto di Scienze Applicate e Sistemi Intelligenti Eduardo Caianiello (CNR), Italy

#### Reviewed by:

Aleksey V. Zaitsev, Institute of Evolutionary Physiology and Biochemistry (RAS), Russia Lucy Maree Palmer, Universität Bern, Switzerland Scott S. Bolkan, Princeton University, United States

#### \*Correspondence:

Brielle R. Ferguson<EMAIL_ADDRESS>Wen-Jun Gao<EMAIL_ADDRESS>
Received: 08 February 2018 Accepted: 17 April 2018 Published: 16 May 2018

#### Citation:

Ferguson BR and Gao W-J (2018) PV Interneurons: Critical Regulators of E/I Balance for Prefrontal Cortex-Dependent Behavior and Psychiatric Disorders. Front. Neural Circuits 12:37. doi: 10.3389/fncir.2018.00037 from family, friends and society. Not surprisingly, disorders with reduced cognitive ability share high comorbidity with decreased sociability, as well as anxiety and depression. What remains unclear is whether these behaviors rely on a common circuitry, and if that circuitry could represent a platform for development of therapeutic interventions.

A large body of research suggests that the prefrontal cortex (PFC) is critical in regulating these behaviors and that the balance between excitatory and inhibitory neurotransmission (E/I balance) plays a fundamental role. The PFC is a unique associative cortical region that receives multimodal inputs from sensory regions, limbic structures and neuromodulatory nuclei (Carmichael and Price, 1996; Hoover and Vertes, 2007), and is taxed with processing these heterogeneous inputs to effectively guide ongoing behavior. The seamless execution of a range of behaviors relies on the integration of past experiences and current goals to select appropriate behavioral programs. Across species, we deem the myriad abilities linked to the PFC as ''executive functioning,'' and this category includes working memory, cognitive flexibility, planning, errormonitoring, decision-making, attention and social cognition. Although we have yet to unravel the precise mechanisms of information representation and processing in the PFC, recent optogenetic and imaging breakthroughs along with decades of lesion, pharmacological, and electrophysiological data highlight certain fundamental principles that may have relevance across PFC-dependent functions. Here, we argue that the E/I balance is essential for the proper execution of a range of PFC-dependent behaviors, and targeting this balance may be an effective approach in disorders that harbor related behavioral impairments.

We will begin with a discussion of GABAergic signaling and E/I balance using foundational research from sensory cortices to inform our interpretation of data from the PFC. Then, we will summarize data supporting our hypotheses, with a primary focus on research in rodents, as currently these types of studies allow for the most experimental flexibility and direct manipulation of specific neuronal types that regulate E/I balance. We acknowledge that the agranular rodent frontal cortex has no apparent structural equivalent of the dorsolateral PFC (dlPFC), a region in humans and primates linked to many of the behaviors that will be discussed in this review. Additionally, the human and non-human primate PFC are substantially more developed and complex relative to the comparable region in the rodent. Nevertheless, using a combination of tracing studies and lesion-symptom mapping, many scientists have concluded that the prelimbic region of the rodent medial PFC (mPFC) serves as the closest functional homolog to the dlPFC in humans and non-human primates (Uylings et al., 2003; Seamans et al., 2008). Human and primate studies will be utilized for comparison and framing purposes, but a detailed review of data from all three species is beyond the scope of this review. For a discussion involving a more comprehensive overview of clinical research discussing E/I balance and implications for psychiatric disease, see the following recent reviews (Foss-Feig et al., 2017; Krystal et al., 2017; Takarae and Sweeney, 2017).

# GABAergic SIGNALING AND E/I BALANCE

Within the cortex, there are two primary types of neurons: glutamatergic excitatory pyramidal neurons and GABAergic inhibitory interneurons. While a network of purely excitatory connections offers little computational complexity, GABAergic interneurons confer a circuit with tremendous flexibility, dynamically modulating the gain of pyramidal neuron responses from simple all-or-none responsiveness (Klausberger et al., 2003; Isaacson and Scanziani, 2011). Through inhibition of neighboring pyramidal neurons, GABAergic interneurons act as an important brake on excitatory signaling, control spike generation and timing (Pouille and Scanziani, 2001), and prevent incoming excitation from afferent structures from causing runaway feed-forward excitation (Dichter and Ayala, 1987).

GABAergic interneurons are a diverse family that can be classified by their morphology, electrophysiological properties, or histological markers (Markram et al., 2004, 2015). The most common nomenclature segregates interneurons into three broad types, parvalbumin (PV), somatostatin (SST) and ionotropic serotonin receptor 5HT3a (5HT3aR) expressing interneurons (Rudy et al., 2011). Although not entirely distinct, this population minimizes overlap between the groups and accounts for nearly the entirety of known interneuronal subtypes.

The most common subtype, PV interneurons, are known for their fast-spiking phenotype, low input resistance, and high-amplitude rapid after-hyperpolarization (AHP; Kawaguchi et al., 1987; Kawaguchi and Kubota, 1997). This combination of properties confers an ability to fire a rapid train of action potentials unlike any other neuron in the cortex. PV interneurons can be further divided into two subtypes: basket cells that innervate the soma and proximal dendrites, and chandelier cells that synapse onto the axon initial segment (Kawaguchi and Kubota, 1997; Petilla Interneuron Nomenclature Group et al., 2008). Relatively less is known about the properties and function of chandelier cells, and data demonstrate their action may not be purely inhibitory at the axon initial segment (Szabadics et al., 2006). PV basket cells are more common, and by virtue of being easier to identify and study, their properties and role in circuit integration in the normal brain have been characterized in greater detail. Data from paired recordings indicate that a single PV interneuron contacts nearly every local pyramidal neuron providing PV basket cells with an unparalleled ability to regulate the activity of nearby pyramidal neurons (Packer and Yuste, 2011). This allows for a tremendous level of feedforward and feedback inhibition that serves several important functions (Hu et al., 2014).

An ability of PV interneurons that appears to be universal across cortical regions is controlling spike timing in neighboring excitatory neurons (Pouille and Scanziani, 2001; Wehr and Zador, 2003). Pyramidal neurons receive constant barrages of excitatory synaptic input or excitatory post-synaptic potentials (EPSPs), and are tasked with deciding whether or not to fire an action potential. Given the density of excitatory connections both from afferent structures and within the local circuit, along with the high level of synaptic divergence, excitation of a sufficient magnitude to generate an action potential will likely recruit a large portion of the cortical network, making it challenging to represent information with any specificity. This creates two computational dilemmas, creating a network state such that all excitatory neurons are not firing constantly, as well as allowing for distinct neurons or groups of neurons to have differential responsiveness to EPSPs. One issue is addressed by feedforward inhibition (FFI), a principle illustrated canonically in thalamocortical circuits. A major source of excitatory cortical input is the thalamus, whose neurons branch and synapse onto both excitatory and inhibitory neurons. Seminal physiological studies demonstrated that when thalamic fibers are electrically stimulated, this causes a direct excitatory response in cortical pyramidal neurons, followed shortly by a strong hyperpolarization due to the activation of neighboring PV basket cells (Agmon and Connors, 1991; Castro-Alamancos and Connors, 1997). This enhances the temporal fidelity of pyramidal neuron responsiveness by limiting the window in which incoming excitation can generate an action potential, and has also been observed in the mPFC (Cruikshank et al., 2012; Delevich et al., 2015). The second issue is mitigated through another property of PV interneuron activity, gain control, or ''the rate at which the firing of a neuron increases in response to increasing excitatory input'' (Isaacson and Scanziani, 2011). This transforms pyramidal neuron responsiveness from a simple go or no-go pattern of activity to having the capacity to represent a broad range of input levels in their firing patterns, a feature that is likely to be crucial in the mPFC, which processes complex multimodal information.

By innervating pyramidal neurons at the soma, PV basket cells are strategically positioned to exert both FFI and gain control. Additionally, PV basket cells have been implicated in the generation of gamma oscillations (30–80 Hz, Buzsáki and Draguhn, 2004), an oscillation range linked to cognition and information processing across species (Howard et al., 2003; Gaetz et al., 2011; Lundqvist et al., 2016).  These interneurons also synapse primarily onto alpha1-gamma-Aminobutyric acid receptors (α1-GABAARs) that are present at high levels on the soma (Fritschy and Mohler, 1995; Nusser et al., 1996), and have a decay constant corresponding to the rising phase of gamma waves (Gonzalez-Burgos and Lewis, 2008). Demonstrating a definitive link, a pair of studies utilized a battery of optogenetic tools to silence or activate cortical PV interneurons, and found them necessary and sufficient in the generation and maintenance of gamma oscillations (Cardin et al., 2009; Sohal et al., 2009). Increasing their activity not only amplified cortical gamma oscillations, but also enhanced information processing by increasing the gain of incoming information and improving behaviors dependent on the manipulated brain regions. Data also indicate the presence of significant gamma alterations in diseases with cognitive impairment, such as schizophrenia and autism spectrum disorders (ASD; Rojas and Wilson, 2014; McNally and McCarley, 2016).

In order for these processes to occur efficiently, levels of inhibition and excitation must remain in the appropriate balance, which in the mPFC has been measured as approximately 20/80% (excitation/inhibition; Le Roux et al., 2008; den Boon et al., 2015). This balance is not static and requires that inhibition be responsive to fluctuations in cortical state and levels of excitatory input (Galarreta and Hestrin, 1998; Shu et al., 2003; Atallah and Scanziani, 2009; Xue et al., 2014). PV interneurons are well suited for helping maintain the proper balance, as their activity is associated with ''divisive'' rather than subtractive inhibition, which would allow PV interneurons to respond proportionally to dynamic excitation within cortical circuits (El-Boustani and Sur, 2014). Within this balance, however, there must be the capacity for change or plasticity. Data suggests that E/I balance helps facilitate the induction of long-term potentiation, which could be critical in allowing the PFC to remain flexible and adapt to new stimulus-response contingencies (Staff and Spruston, 2003; Marder and Buonomano, 2004).

## PHARMACOLOGICAL MANIPULATION OF E/I BALANCE AND BEHAVIOR

Two cognitive abilities that are among the most consistently associated with damage to the PFC are working memory (Kolb et al., 1974; Stokes and Best, 1990; Kesner et al., 1996) and cognitive flexibility (Ragozzino et al., 1999; Birrell and Brown, 2000; Block et al., 2007; Bissonette et al., 2008). Working memory is the ability to maintain, manipulate and recall information to guide behavior (Baddeley, 1992; Goldman-Rakic, 1994), while cognitive flexibility involves the capacity to update strategies in order to obtain a reward (Bissonette et al., 2013). Both represent key building blocks of higher-level cognitive processes and are often severely impacted in the context of psychiatric disease. Social cognition is similar in denotation to working memory, but applies the principles of acquiring, storing, and manipulating information specifically to the context of interacting in a flexible and appropriate manner with conspecifics (Adolphs, 1999). Supporting its reliance on the frontal cortex are lesion data from humans and non-human primates indicating that damage to the ventromedial and orbital frontal cortex consistently result in abnormal social behavior (Ackerly and Benton, 1948; Butter and Snyder, 1972).

PFC infusion of bicuculline, a GABAAR antagonist, disrupts working memory and cognitive flexibility in multiple species (Sawaguchi et al., 1988, 1989; Enomoto et al., 2011; Paine et al., 2011; Auger and Floresco, 2015). In vivo multi-unit recordings in primates revealed that following bicuculline administration (which would increase the E/I balance), putative FS interneurons that were not linked to any specific task epochs of a spatial working memory task begin to show spatial tuning. Additionally, pyramidal neurons that once exhibited directional sensitivity, began to respond to their non-preferred directions, and neurons that had no task-modulated activity respond erroneously to random directions and task epochs (Rao et al., 2000). In concert, these activity alterations would greatly reduce the signalto-noise ratio, severely limiting the ability of the circuit to generate and relay correct motor commands with precision. Circuits supporting normal patterns of social interaction may be regulated similarly, as bicuculline infusion into the rat mPFC reduces the number and duration of interactions as well as weakens the preference for a social stimulus compared to a non-social object (Cochran et al., 2002).

Conversely, studies that utilize enhancement of GABAergic transmission have been informative, as this approach would also disrupt the E/I balance. Muscimol, a GABAA agonist, is used commonly as a means for transient inactivation of brain regions. Direct infusion in the mPFC in rodents has resulted in impairments in working memory (Urban et al., 2014), retention of strategy-switching behavior (Rich and Shapiro, 2007) and social play behavior (van Kerkhof et al., 2013, 2014).

# OPTOGENETICS, IMAGING AND PHYSIOLOGY HIGHLIGHT THE ROLE OF GABAergic SIGNALING IN PFC-DEPENDENT BEHAVIOR

Advances in tools available for bidirectional modulation of neuronal activity with cell-type specificity have expanded our understanding of the role of PV interneurons in regulating PFC E/I balance and PFC-dependent behaviors. Additionally, the development of genetically encoded calcium indicators, like gCaMP (Chen et al., 2013), allows researchers to monitor activity or calcium transients selectively in PV interneurons or other neuronal subtypes. With these technologies, researchers can catalog the entire population of mPFC interneurons from learning and strategy acquisition through successful performance. Many findings confirm and fine-tune long-standing hypotheses about how PV interneurons regulate pyramidal neuron activity, oscillatory dynamics, and cortical information representation, but allow for more detailed investigation of their causal role in behavior with enhanced temporal and cellular control.

One extremely useful feature of optogenetics (Deisseroth, 2011) has been the in vivo identification of different neuronal subtypes, or ''optical tagging'' (Zhao et al., 2011; Roux et al., 2014), allowing researchers to confirm the circuit consequences of PV interneuron activation as well as how these interneurons contribute to particular behaviors using recording or imaging. For example, optical tagging revealed that PV interneurons provide rapid synchronous inhibition of mPFC excitatory neurons in contrast with the more variable and less forceful inhibition provided by SST interneurons (Kvitsiani et al., 2013). A similar approach was used to demonstrate that PV interneurons increase their firing during goal-directed behavior. In this same study, researchers also showed that increases in PV activity were correlated with inhibition of certain pyramidal neurons and enhancement of others, potentially representing ensembles of relevant vs. irrelevant information for successful responding respectively (Kim H. et al., 2016). Increases in PV activity have also been implicated in the extinction of cue-based responding (Sparta et al., 2014), a process that may mediate rule-shifts during cognitive flexibility tasks.

During working memory assays, PV interneurons in the mPFC fire prominently throughout short delays (5 s or less) between sample and choice presentation (Kim D. et al., 2016), and show preferential firing toward go trials in go, no-go task (Kamigaki and Dan, 2017). However, more complicated tasks reveal they may encode information about distinct task epochs. Lagler et al. (2016) utilized a multidimensional delay task, where a differential gustatory cue signaled the location of reward in a Y-maze, and quantified PV interneuron activity during discrete phases of the assay. Interestingly, mPFC PV interneurons exhibited episodic modulation, and a large group showed preferential activation during the ''goal run, '' when the mouse was proceeding through the base of the Y-maze before making a choice (Lagler et al., 2016). PV interneurons in the mPFC also increase their firing in response to typical social interactions, compared to the investigation of a non-social object (Selimbeyoglu et al., 2017).

Modulation of PV activity has also been useful in demonstrating one of the more obvious principles of E/I balance, that too much inhibition can also impair PFC circuit function. Optogenetic or pharmacogenetic activation of PV interneurons has been an effective method of disrupting the function of particular brain regions, suggesting that alterations in E/I balance in either direction can alter behavior. We recently demonstrated that increasing the excitability of mPFC PV interneurons (thereby reducing the E/I balance) significantly impairs working memory, cognitive flexibility and social interaction (Ferguson and Gao, 2018). Similarly, optogenetic activation of PV interneurons or GABAergic neurons more broadly disrupts delay-related activity and working memory performance (Rossi et al., 2012; Liu et al., 2014). These data corroborate pharmacological studies using muscimol, a GABA agonist, as a method of inactivating PFC to elucidate its role in various behaviors (Narayanan and Laubach, 2006; Urban et al., 2014). Correspondingly, working memory and social interaction deficits can be recapitulated by enhancing pyramidal neuron excitability (enhancing E/I balance; Yizhar et al., 2011; Liu et al., 2014; Kamigaki and Dan, 2017).

# BUILDING THE mPFC MICROCIRCUIT

In combination with numerous elegant studies from the visual, auditory and somatosensory systems, we can begin to develop a picture of how the mPFC microcircuit might be organized to represent and process information relevant for executive function. The largest source of afferents originates from the mediodorsal thalamus (MD), but mPFC neurons also receive convergent inputs from limbic structures including the hippocampus and amygdala, that may convey spatial and emotional information. Additionally, the mPFC has dense innervation from neuromodulatory regions, the ventral tegmental area (VTA), locus coeruleus (LC), dorsal raphe (DR) and basal forebrain (BF), which may influence the state of the animal including level of attention, arousal, or the salience of the current task or resulting rewards or punishment (Kuroda et al., 1996; Hoover and Vertes, 2007). Using working memory as an example, early studies in rodents found that neurons representing differential task-related information in a consistent manner were sparse (Jung et al., 1998, 2000; Baeg et al., 2003).

However, more recent work has identified numerous mPFC neurons that change their firing frequency during distinct phases of a delayed-alternation task. Of these neurons, Yang et al. (2014) was able to separate these into units showing preferential firing on left vs. right trials, and encode information such as the choice, delay, and presence of reward. Other studies highlight a segregation of spatial and rule information encoding in mPFC neurons, of which the rule representation depends critically on the activity of the MD (Bolkan et al., 2017; Schmitt et al., 2017). Taken together, this suggests that pyramidal neurons in the mPFC form functional units that represent different features relevant to successful working memory performance and other PFC-dependent behaviors (**Figure 1**).

How might the appropriate E/I balance be maintained in the mPFC? As a general principle, thalamic inputs provide stronger activation to PV interneurons vs. pyramidal cells in sensory cortices (Hestrin, 1993; Kloc and Maffei, 2014). Similarly, enhancing excitability of MD neurons results in enhanced spiking of PV interneurons in mPFC, and not pyramidal cells (Schmitt et al., 2017), while dampening MD activity increases E/I balance by selectively reducing inhibitory currents onto pyramidal cells (Ferguson and Gao, 2018). Thus, the MD seems to function drive PV interneurons in the mPFC to maintain an E/I balance biased towards inhibition. Maintaining a network state with a high level of inhibition, would strengthen computational capacity through enhancing temporal precision and increasing the gain of pyramidal neuron responsiveness, as well as enhancing the signal to noise ratio of neurons encoding task or rule representations (**Figure 1**). Supporting that this large-scale inhibition is mainly PV interneuron mediated, we observed that enhancing PV interneuron excitability can fully normalize shifts in E/I balance with corresponding behavioral improvements (Ferguson and Gao, 2018).

PV interneurons are powerful regulators of pyramidal neuron activity and appear to be the most vulnerable across psychiatric disorders, but research highlights other major inhibitory interneurons that help fine-tune circuits supporting PFC-dependent behaviors (Wang et al., 2004). Axo-somatic inhibition accounts for only up to 10% of total inhibition received by excitatory neurons, while the remaining 90% occurs at proximal and distal dendrites and spines and is mediated by other inhibitory neuronal subtypes (Kubota et al., 2016). For example, SST interneurons exhibit phasic firing during distinct phases of a delayed alternation task, and seem to show segregation of firing among left vs. right trials (Kim D. et al., 2016). Further, in the hippocampus, SST interneurons provide powerful dendritic inhibition that helps to regulate both synaptic integration as well as the gain of pyramidal neuron responsiveness (Lovett-Barron et al., 2012). Additionally, vasoactive-intestinal peptide (VIP, 5HT3aR expressing) interneurons (Férézou et al., 2002), which inhibit both SST and PV interneurons seem to serve an important disinhibitory function in cortical circuits (Pi et al., 2013), that may help modulate the gain of pyramidal neuron responses (Fu et al., 2014). Recent data demonstrates VIP interneuron activity enhances mPFC circuit performance during two delay tasks involving working memory (Kamigaki and Dan, 2017). VIP interneurons may also facilitate forming associations in pyramidal neurons, as VIP interneurons show the greatest increase in activity along with pyramidal neurons in a go/nogo task following punishment (Pinto and Dan, 2015). As they inhibit both PV and SST cells, disinhibition of pyramidal cells by VIP interneurons may be a mechanism for facilitating the activity of neuronal ensembles relevant for optimal task performance in multiple PFC-dependent tasks (**Figure 1**). Further, targeting VIP interneuron activity has shown promising results in an animal model of hypofrontality in schizophrenia (Koukouli et al., 2017).

#### TARGETING E/I BALANCE IN ANIMAL MODELS TO IMPROVE BEHAVIOR

Multiple approaches for enhancing the normal function of prefrontal PV interneurons have been effective in mitigating a constellation of behavioral symptoms animal models of cognitive dysfunction and social abnormalities. For example, altering expression of Dlx5/6, a protein that regulates PV interneuron development leads to cognitive and other behavioral abnormalities in mice. Optogenetic activation of mPFC PV interneurons in Dlx5/6 mice is sufficient to ameliorate cognitive flexibility and social interaction deficits (Cho et al., 2015). Similarly, mice deficient in Contactin-associated protein-like (CNTNAP)–2 display marked social impairments that can be rescued by increasing PV interneuron activity optogenetically (Selimbeyoglu et al., 2017) in agreement with findings from naïve mice with optogenetically disrupted mPFC E/I balance (Yizhar et al., 2011). Dampening activity in the MD (modeling findings of mediodorsal hypofunction in schizophrenia) has been shown to impair working memory, cognitive flexibility, social interaction and alter anxiety-related behavior (Parnaudeau et al., 2013, 2015; Bolkan et al., 2017; Ferguson and Gao, 2018). However, pharmacogenetic activation of mPFC PV interneurons normalizes all mPFC-dependent behavioral deficits (Ferguson and Gao, 2018), highlighting the therapeutic potential of PV interneurons as a strategy for repairing E/I balance and treating a range of behavioral deficits. These results are complemented by pharmacological augmentation of GABAergic signaling, most commonly with benzodiazepines, that has been shown to repair behaviors including social interaction and cognitive flexibility (Wen et al., 2010; Han et al., 2012; Cho et al., 2015).

An ongoing question has been whether altering the properties of PV interneurons is responsible for the initial manifestation of behavioral deficits in psychiatric disease and animal models, in particular, those beyond the cognitive realm. For example altering the excitability of PV interneurons by genetic deletion of N-methyl-D-aspartate receptors (NMDARs), leads to selective deficits in working memory and associative learning (Carlen et al., 2012), while a broader deletion in various subtypes of cortical and hippocampal interneurons, induces cognitive impairments along with anxiety-related behavior, depressive symptoms and social impairments (Belforte et al., 2010). Similarly, reducing the expression of the voltage-gated sodium channel (NaV) 1.1, a channel linked to autism-related behaviors in Dravet syndrome, broadly in forebrain GABAergic interneurons induces spatial memory and social interaction deficits, along with increased repetitive behaviors and anxiety (Han et al., 2012). Optogenetic suppression of PV interneuron activity recapitulates only cognitive impairments, specifically a deficit in extradimensional set-shifting (Cho et al., 2015). Additionally, selective expression of tetanus toxin light chain (TeLC) in mPFC PV interneurons, disrupts working memory and cognitive flexibility, while sparing behaviors representing positive and negative symptoms in schizophrenia (Murray et al., 2015). However, in this study, recording of mPFC local field potentials revealed gamma oscillations did not differ between groups, so without further physiological analysis, the ultimate circuit consequences are somewhat ambiguous. Although it is unclear whether disrupting PV activity can disrupt a broader range of behaviors given different findings with articles employing disparate methodologies, the specificity of the insult seems to be critical. Solely targeting PV interneurons will reliably recapitulate phenotypes of cognitive disruption, while broader interneuronal insults increase the likelihood of impacting other behaviors. Another possibility is that the magnitude of the PV interneuron deficit may correlate directly with the propensity for disrupting behaviors beyond cognition.

### E/I BALANCE IN PSYCHIATRIC DISEASE

The pattern of rodent PFC lesions and deficits mirrors the impairments observed following damage to the homologous structures in primates and humans, suggesting the functions of these brain regions are highly conserved across species. Even more compelling is how myriad findings recapitulate commonly observed morphological alterations and behavioral endophenotypes seen in psychiatric disease. One of the most demonstrative parallels is found in the pathophysiology and symptom sequelae of schizophrenia. Schizophrenia is a debilitating psychiatric disease that affects approximately 1.1% of the world's population, and is characterized by positive, negative, and cognitive symptoms (Regier et al., 1993). Working memory and cognitive flexibility represent core dysfunctions in schizophrenia that remain intractable by treatment with current antipsychotics (Lee and Park, 2005; Insel, 2010). These along with a range of additional cognitive symptoms emerge prior to the onset of psychosis in early adulthood, and are the largest predictor of functional outcome in individual patients (Green, 1996). Negative symptoms involve anhedonia and social withdrawal, and are also minimally responsive to presently available treatments (Remington et al., 2016).

Analogous to experimental observations from animal studies, parallel deficits in PFC GABAergic signaling (Lewis et al., 2012), and E/I balance (Lisman, 2012) are both implicated in the pathophysiology of cognitive dysfunction in schizophrenia. Functional imaging reveals that schizophrenics exhibit patterns of hypofunction in the PFC (Van Snellenberg et al., 2016) as well as reduced volume in frontal cortex postmortem tissue (Selemon et al., 2002). Copious research demonstrates reductions in the marker for the GABA synthesizing enzyme, GAD-67, in PV interneurons in the PFC (Akbarian and Huang, 2006; Lewis et al., 2012). Levels of GAD-67 are activity-dependent (Benson et al., 1989; Sloviter et al., 1996), and coupled with the common finding that PV levels are also decreased within these interneurons (Glausier and Lewis, 2017), it suggests that the function of PFC PV interneurons is diminished in schizophrenia. Correspondingly, patients with schizophrenia show decreases in task-evoked in gamma oscillations, an oscillation band believed to be dependent on PV interneuronal firing, that correlate with the level of functional impairment in working memory (Basar-Eroglu et al., 2007).

Analogous to human pathology, rodent models of cognitive dysfunction in schizophrenia show an extremely high prevalence of reductions in PV or GAD67 expression in the mPFC. As more evidence reveals this disorder likely stems from heterogeneous etiologies, pharmacological, environmental, and transgenic models exhibit this common feature (Cochran et al., 2002; Francois et al., 2009; Carlson et al., 2011), suggesting that GABAergic hypofunction represents a convergence point. The constellation of data from human populations, animal models, experimental disruption of GABAergic signaling, and electrophysiological recordings across these contexts, highlight this system as a pathway to intervention in afflicted individuals.

Intriguingly, this pattern repeats among multiple psychiatric disorders including ASD, depression, and intellectual disability (Gao and Penzes, 2015; Luscher and Fuchs, 2015). ASD in particular harbors significant overlapping behavioral impairments and underlying neurobiological alterations. Autistic patients also exhibit prominent deficits in executive function, including working memory impairments and behavioral inflexibility (Hughes et al., 1994). Social abnormalities are also a behavioral hallmark of autistic pathology, primarily manifested as a deficit in non-verbal communication (Mundy et al., 1986). Autistic patients also suffer from a greater level of anxiety and depression relative to general population (Kim et al., 2000).

E/I balance disruption has also emerged as a prominent hypothesis in ASD (Rubenstein and Merzenich, 2003). Dysregulation of GABAergic signaling has been implicated in the etiology of ASD, and disorders sharing high comorbidity with the disease including anxiety and epilepsy (Coghlan et al., 2012). The chromosomal region 15q11-q13 is comprised of multiple genes encoding subunits of the GABAAR. Microduplications in this region (Cook et al., 1998; Menold et al., 2001) or aberrant expression of the associated gene products (Hogart et al., 2007; Mendez et al., 2013) have been frequently observed in ASD clinical populations. In addition to the finding of copy number variations in the 15qllq13 chromosomal locus, the majority of autism-linked genes are preferentially expressed in interneurons (Xu et al., 2014).

Supporting this are in vivo findings of reduced GABA and GABAAR levels in the frontal cortex of autistic patients (Harada et al., 2011; Mori et al., 2012), along with reductions in gamma oscillations (Sun et al., 2012). Similar deficiencies in inhibitory neurotransmission have been reported in mice with mutations in ASD-linked genes, mice that also exhibit relevant behavioral impairments (Peñagarikano et al., 2011; Han et al., 2012). For example, T(+)Itpr3(tf)/J (BTBR) mice, a model of idiopathic autism, exhibit decreased GABAergic currents and increased excitatory neurotransmission in the hippocampus, indicating a shift in the E/I balance. Concurrent with these physiological changes, these mice display reductions in sociability, cognitive impairments, and alterations in anxiety-related behaviors (Han et al., 2014). Additionally, mouse models with disruption of genes such as methyl-CpG-binding protein-2 (MECP2), Scn1a+/− and CNTNAP–2, all harbor interneuron deficits along with behavioral impairments (Chao et al., 2010; Peñagarikano et al., 2011; Han et al., 2012). Modeling of Fragile X syndrome, a genetic disorder associated with intellectual disability, also reveals cognitive impairments are associated with prominent GABAergic hypofunction in mice (Selby et al., 2007; Curia et al., 2009).

# GABAergic SIGNALING AS A THERAPEUTIC TARGET

Multiple studies suggest that augmenting GABAergic signaling via PV interneuron modulation can be effective in ameliorating deficits in working memory, cognitive flexibility and sociability in animal models of psychiatric disease. This has substantial implications for schizophrenia, given that cognitive symptoms are treatment-resistant with both typical and atypical antipsychotics, and are a significant predictor of quality of life in individual patients (Green, 2006). These results complement a rich history of data and hypotheses surrounding GABAergic hypofunction in the human PFC representing a final common pathway in the cognitive symptoms in schizophrenia (Lewis, 2014). Yet, we still await the discovery of a GABAergic modulator that displays effectiveness in alleviating cognitive dysfunction across multiple disorders.

Of the recent clinical trials targeting cognitive symptoms, very few utilize GABAergic drugs, relative to agonists and modulators of glutamatergic signaling. Although benzodiazepines have shown potential in treating psychosis (Wolkowitz and Pickar, 1991; Carpenter et al., 1999), their potential in reducing cognitive impairments remains largely unexplored with the exception of a few studies. Interestingly in one study, a benzodiazepine, lorazepam, exacerbated working memory impairments in schizophrenic patients and healthy controls and altered activity within networks supporting cognition. Researchers concluded that hyper-inhibition was responsible for cognitive dysfunction, and correspondingly, flumazenil an antagonist of the benzodiazepine site attenuated deficits (Menzies et al., 2007). However, a separate study found iomazenil, a flumazenil analog, increased psychoses in schizophrenic patients (Ahn et al., 2011). Given the high comorbidity of ASD with epilepsy (Canitano, 2007), researchers have gleaned that drugs effective in reducing seizure activity that typically works through elevating GABA levels also mitigate autistic symptoms (Jambaqué et al., 2000). Promisingly, GABA agonists have reached Phase II clinical trials for the treatment of social disability in ASD, but so far have only focused on α2/α3-GABAARs.

The synapse between PV-expressing chandelier cells and the axon initial segment have become frequent but unsuccessful target for therapeutics in schizophrenia. In patients with schizophrenia, chandelier cell axon terminals exhibit decreased levels of the reuptake enzyme, GABA transporter 1 (GAT1) mRNA, along with higher expression of α2-GABAARs on the axon initial segment of pyramidal neurons (Volk et al., 2002). These alterations would result in less GABA reuptake, and increased post-synaptic GABAergic inhibition respectively, which together would serve to augment GABAergic neurotransmission at these particular synapses. This has been interpreted as a compensatory mechanism, implying these synapses may be the site of the initial GABAergic deficit (Volk and Lewis, 2005). However, despite promising results from early studies (Lewis et al., 2008), larger sample sizes have yielded no differences between schizophrenic patients treated with an α2/α3-GABAAR agonists compared to those treated with the placebo (Buchanan et al., 2011). Given the failure of enhancing α2-GABAAR-mediated signaling in improving cognitive dysfunction in clinical trials as well as in the context of compensatory mechanisms in schizophrenia, it is reasonable to conclude these receptors are not a viable therapeutic target.

One difficulty in utilizing GABAergic therapeutics, including benzodiazepines, is their sedative properties, which are mediated by their action at the α1 GABAAR subunit (Löw et al., 2000). This may account for the lack of exploration of α1-GABAAR modulators in clinical trials. However the actions of dopamine across different subjects and contexts provide a compelling example of how optimal levels of neurotransmitter can be associated with normal cognition, attention and alertness, but levels either too high or low, can result in impaired cognition, inability to focus, and drowsiness (Cools and D'Esposito, 2011). Our findings using indiplon, an α1-GABAAR positive allosteric modulator (Ferguson and Gao, 2018), and other studies using benzodiazepines (Han et al., 2012, 2014) suggest that the therapeutic window for influencing cognitive function may be distinct from doses that induce sedation. It is likely that GABAergic signaling in the PFC and its correlated functions follow an inverted-U trajectory, and with better biomarkers for individual differences in GABA levels, indiplon, novel α1- GABAAR modulators or benzodiazepines may harbor significant therapeutic potential.

# FUTURE DIRECTIONS

Additional questions remain that could help inform a circuit model of information processing in the mPFC. For example, how organized are the actions of PV interneurons in the mPFC? Are they providing a general blanket inhibition that helps maintain a cortical state conducive to information processing due to their properties described above? Evidence supporting this is imaging of PV interneurons revealing that PV interneurons show the largest modulation, and have uniform firing patterns that persist throughout delay periods during PFC-dependent tasks in comparison to other interneuron subtypes (Pinto and Dan, 2015; Kim D. et al., 2016). If PV interneurons do provide indiscriminate inhibition throughout the mPFC, to what extent to disinhibitory circuits through VIP and SST interneurons help further enhance information representation, and through what mechanisms and inputs? Other data indicate segregated groups of PV interneurons encode distinct task phases (Lagler et al., 2016), suggesting instead of blanket inhibition, patterns of activity are differentially regulated across different PV interneurons. Future studies should also further explore the actions of PV chandelier cells in the mPFC during behavior. New research indicates this subtype may preferentially inhibit amygdala-projecting pyramidal neurons in the mPFC, suggesting chandelier cells have a distinct role in the mPFC microcircuit from PV basket cells (Lu et al., 2017).

Currently available imaging and physiological methods provide promising approaches for measuring changes in neuronal activity that are correlated with successful behavior. However, given the diversity of neocortical interneuron subtypes and function, how PV interneurons and other cell types contribute to modulation of E/I balance, and whether and how modulation of each subtype impact different PFC-dependent behaviors warrants further investigation. We still lack effective treatments for cognitive deficits, social interaction impairments, and other associated behavioral dysfunctions in numerous psychiatric disorders, underscoring the importance of continuing to unravel how these behavioral processes occur under normal conditions. If researchers can identify how to effectively manipulate PFC circuit activity, we can develop reliable strategies for engineering optimal patterns of cortical activity to ameliorate performance. Ultimately, if successful manipulations can be linked to physiological signatures that can be observed using less invasive recording approaches such as oscillatory patterns, we can potentially determine biomarkers for successful behavioral therapies in humans.

### AUTHOR CONTRIBUTIONS

BF wrote the manuscript and W-JG edited it.

#### REFERENCES

## FUNDING

This work was supported by National Institute of Health (NIH)/National Institute of Mental Health (NIMH), F31MH 111361 to BF and R01MH085666 and R21MH111609 to W-JG and Helen S. Vernik Schizophrenia Pilot Research Project from the Department of Psychiatry, Drexel University College of Medicine to BF, Linda Chamberlin and W-JG (2017).

the treatment of cognitive impairments in people with schizophrenia. Biol. Psychiatry 69, 442–449. doi: 10.1016/j.biopsych.2010.09.052

schizophrenia. Biol. Psychiatry 80, 617–626. doi: 10.1016/j.biopsych.2016. 02.017

for neurogenetic disorders. J. Neurosci. 34, 1420–1431. doi: 10.1523/jneurosci. 4488-13.2014

**Conflict of Interest Statement**: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Ferguson and Gao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# How Memory Conforms to Brain Development

#### Ana P. Millán, Joaquín J. Torres\* and Joaquín Marro

Institute "Carlos I" for Theoretical and Computational Physics, University of Granada, Granada, Spain

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.

#### Edited by:

Vito Di Maio, Institute of Applied Sciences and Intelligent Systems (ISASI), Italy

#### Reviewed by:

Andreas Knoblauch, Hochschule Albstadt-Sigmaringen, Germany Christian Tetzlaff, Max-Planck-Institute for Dynamics and Self-Organisation, Max Planck Society (MPG), Germany

#### \*Correspondence:

Joaquín J. Torres<EMAIL_ADDRESS>
Received: 13 September 2018 Accepted: 26 March 2019 Published: 16 April 2019

#### Citation:

Millán AP, Torres JJ and Marro J (2019) How Memory Conforms to Brain Development. Front. Comput.  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-evolving brain 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 weights has 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., 2007, 2008; Torres 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.

#### 2. 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(t), with elements eij(t) = {0, 1}, in which each node represents a neuron and each edge [eij(t) = 1] stands for a synapse. The degree of node i at time t is defined as

$$k\_i(t) = \sum\_{j=1}^{N} e\_{ij}(t) \tag{1}$$

and the mean degree of the network is

$$\kappa(t) = \frac{1}{N} \sum\_{i=1}^{N} k\_i(t). \tag{2}$$

Following a familiar (Amari-Hopfield) prescription (Amari, 1972; Amit, 1989), each neuron i is modeled as a stochastic binary unit, si(t) = {0, 1} (representing respectively a silent and a firing neuron), whose state evolves in time according to the probabilistic dynamics

$$P\left[s\_i(t+1) = 1\right] = \frac{1}{2}\left\{1 + \tanh\left[T^{-1}\left(h\_i(t) - \theta\_i(t)\right)\right]\right\},\tag{3}$$

where

$$h\_i(t) = \sum\_{j=1}^{N} \omega\_{ij} e\_{ij}(t) s\_j(t) \tag{4}$$

is the local field at neuron i quantifying the incoming input from neighbor neurons and

$$\theta\_i(t) = \frac{1}{2} \sum\_{j=1}^{N} w\_{ij} e\_{i\vec{j}}(t) \tag{5}$$

is the neuron'sthreshold 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; Azouz and 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 eijwij 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 ortemperature 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, wij, 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),

$$\begin{aligned} \boldsymbol{w}\_{i\bar{j}} &= [\boldsymbol{\kappa}\_0 \boldsymbol{a}\_0 \,(1 - \boldsymbol{a}\_0)]^{-1} \sum\_{\mu=1}^p \left(\boldsymbol{\xi}\_i^{\mu} - \boldsymbol{a}\_0\right) \left(\boldsymbol{\xi}\_j^{\mu} - \boldsymbol{a}\_0\right), \; i \neq j \\ \boldsymbol{w}\_{i\bar{i}} &= \boldsymbol{0}, \end{aligned} \tag{6}$$

where κ0 = κ(t = 0) and a0 is the mean activation of the patterns, i.e., a0 = (NP) −1 PP µ=1 PN 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 wij = wji 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

$$m^{\mu}(t) = \left[\text{Na}\_0(1 - a\_0)\right]^{-1} \sum\_{i=1}^{N} \left(\xi\_i^{\mu} - a\_0\right) s\_i. \tag{7}$$

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 ∀ν 6= µ. 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

$$m^{\mu}(k,t) = \left[Np(k,t)a\_0(1-a\_0)\right]^{-1} \sum\_{i=1}^{N} \left(\xi\_i^{\mu} - a\_0\right) s\_i \delta\_{k,k\_i} \tag{8}$$

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, PN k=1 p(k, t) = 1 ∀t). Therefore, mµ(t) = PN 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 < Pc = 0.138N, and the system is in what is called the memory phase. Pc defines the maximum storage capacity of the network (Amit, 1989), that is, the maximum quantity of information—or number of patterns that 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 wij, which can destabilize such memory phase. Therefore, if P is further increased above Pc 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 > Tc = 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 < Tc , so that memory is recovered but with more errors than in a fully connected network. However, the critical temperature diverges, Tc → ∞ 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 a0 6= 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 a0 = 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 a0 (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

$$m\_1^{\mu}(t) \equiv \frac{1}{N} \sum\_{i=1}^{N} s\_i(t) \xi\_i^{\mu},\tag{9}$$

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) > mth and 0 otherwise, so that **mB**(**t**) = m1 B (t), m2 B (t), ..., mP B (t) indicates, in a binary code, which combination of patterns is recovered at time t. Equivalently, the decimal variable ds can be defined,

$$d\_s(t) \equiv \sum\_{\mu=1}^{P} 2^{\mu-1} m\_B^{\mu}(t),\tag{10}$$

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 Pr 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

$$m^{\mu} = \left[\text{Na}\_{0}(1 - a\_{0})\right]^{-1} \sum\_{\upsilon=1}^{P} \sum\_{i=1}^{\mathcal{N}} \left(\xi\_{\mathcal{N}\upsilon - \mathcal{N} + i}^{\mu} - a\_{0}\right) s\_{\mathcal{N}\upsilon - \mathcal{N} + i},\tag{11}$$

where N = a0N = 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 Pr − 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 mµ = 1 − P −1 −1 - 1 − P −1 − (Pr − 1) P −1 . This yields

$$m^{\mu} = 1 - (P\_r - 1) \left(P - 1\right)^{-1} \le 1,\tag{12}$$

which only meets the equality in the case Pr = 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,

$$\begin{array}{lcl}P\_i^{\mathbb{g}} & = \iota(\kappa)\pi(I\_i),\\P\_i^{l} & = d(\kappa)\eta(I\_i),\end{array} \tag{13}$$

where Ii = hi − θ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 π(Ii) and η(Ii) 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 π(Ii) and η(Ii). 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 eij, 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 π(Ii) for a gain (or η(Ii) 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

$$\begin{array}{lcl}\widetilde{\pi}\_{i} &=& \frac{1}{2} \Big[\pi(I\_{i}) + \frac{1}{N}\Big],\\\widetilde{\eta}\_{i} &=& \frac{1}{2} \Big[\eta(I\_{i}) + \frac{k\_{i}}{\kappa N}\Big],\end{array} \tag{14}$$

where the 1/2 factor is included to assure normalization. Following our previous work, we consider

$$\begin{array}{lcl}\widetilde{\pi}\_{i} &= \frac{I\_{i}^{\alpha}}{\langle I^{\alpha} \rangle N},\\\widetilde{\eta}\_{i} &= \frac{I\_{i}}{\langle I \rangle N},\end{array} \tag{15}$$

which are normalized over the network, PN i=1 eπi = PN i=1 eηi = 1. The power-law relation in eπi allows us to explore both suband super-linear responses by just modifying a single parameter, namely α. The probabilityeη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

$$\begin{array}{lcl}\pi\left(I\_{i}\right)&=\max\left\{2\frac{I\_{i}^{a}}{\langle I^{a}\rangle N}-\frac{1}{N},0\right\},\\\eta\left(I\_{i},k\_{i}\right)&=\max\left\{2\frac{I\_{i}}{\langle I\rangle N}-\frac{k\_{i}}{\kappa N},0\right\},\end{array} \tag{16}$$

which hold that π(Ii), η(Ii P , ki) > 0 ∀i and normalization, N i=1 π(Ii) = PN i=1 η(Ii , ki) = 1.

We impose further restrictions on the network. First of all, eij 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 wij. Moreover, we set the minimum degree of the network, ki = 1, so that there cannot be any disconnected nodes, and we forbid self-connections, eii = 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

$$\frac{d\kappa(t)}{dt} = 2\left[\mu\left(\kappa(t)\right) - d\left(\kappa(t)\right)\right].\tag{17}$$

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

$$\begin{array}{l} u(\kappa(t)) = \max\left\{\frac{n}{N}\left(1 - \frac{\kappa(t)}{2\kappa\_{\infty}}\right), 0\right\} \\ d(\kappa(t)) = \frac{n}{N}\frac{\kappa(t)}{2\kappa\_{\infty}}, \end{array} \tag{18}$$

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),

$$\kappa(t) = \kappa\_{\infty} \left[ 1 - \left( 1 - \kappa\_0 / \kappa\_{\infty} \right) e^{-t/\tau\_p} \right],\tag{19}$$

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

$$\mathbf{g}(t) = \exp(-\sigma^2(t)/\kappa^2(t)),\tag{20}$$

where σ 2 (t) = hk 2 (t)i − κ 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,k1 . 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, hs , 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 hs = 10 MCS and n = 10. Measures on the stationary state of the system are carried out by temporal averages of the macroscopic variables, ¯ f = 1 Pt0+1 t=t0 f(t).

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.

#### 3. 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).

# 3.1. 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 gP ≡ Pr/P (where Pr is the number of retrieved patterns), as a function of α and P. A value gP = 1/P indicates a pure memory state, whereas larger values correspond to mixtures and SG-like states (Amit, 1989), and gP = 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 mP and, finally, the bottom panels show the stationary mean homogeneity, g¯.

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

gP, retrieved after a given transient; the second one (D–F) shows the average overlap with the recovered patterns, mP; and finally the third one (G–I) shows the stationary homogeneity g¯. 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 gp and a high value of mp, indicated by a yellow color. A SG phase appears as a blue region in gp and a lower value of mp, indicated by a green color, whereas a noisy phase appears as black in gp and mp. Similarly, homogeneous structures take place for high values of g¯, indicated by a yellow region in the corresponding diagram, whereas heterogeneous structures are for low values of g¯, 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 gP and mP take intermediate values; the former since only a finite number of patterns is retrieved, gp < 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, mp < 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 g¯ 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 gP → 0 (black region in **Figure 1C**), mP → 0 (black region in **Figure 1F**), and g¯ → 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 mP ≈ 1, and the network structure becomes heterogeneous since α > 1 (see inset of **Figure 2C**). As a consequence, gP → 1/P, (dark-blue region in **Figure 1C**), mP approaches 1 (yellow region in **Figure 1F**), and g¯ → 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 mP remains close to 1 (**Figure 1F**) and gP 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**, indicating g¯ → 0).

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 gp, mp and g¯ 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

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 = 106 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.

# 3.2. 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 3 | Phase diagrams of the system with respect to α and T for four different values of P, in particular for P = 5, 10, 15, and 20, respectively from left to right, and for three values of κ∞ = 20, 40, 60, respectively from top to bottom (A–L). In each panel we show three diagrams: gp, mp, and g¯, as indicated in the label of the color bar. Pink stars in (A–D) indicate the (T, α) point of the corresponding time series in 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 gp and a high value of mp, indicated by a yellow or green color. A SG phase appears as an orange region in gp and a lower value of mp, indicated by a green or blue color, whereas a noisy phase appears as black in gp and dark-blue in mp. Finally, the oscillatory phase appears for high values of gp, (light yellow regions in the corresponding diagrams) and relatively low values of mp (associated blue regions of the corresponding diagrams). Similarly, homogeneous structures take place for high values of g¯, indicated by a yellow region in the corresponding diagram, whereas heterogeneous structures are for low values of g¯, 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 gp 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 gP, the number of patterns visited by the system after the transient

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 g¯ → 0 for α > 1. Notice also in the middle graph of **Figure 3A** how, as T is increased, the overlap corresponding to these states, mP, 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 g¯ → 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, gP → 1 (orange and yellow region of the top panel of **Figure 3A**), indicating that a great number of patterns are being retrieved (Pr → P) with a moderate value of the overlap mP, 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, gP 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 mP. 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 gp 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 g¯ < 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 κ∞.

### 3.3. 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 g¯ shown in **Figures 1**, **3** for α < 1, indicating that the variance of ki 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 by g¯ → 0 (black regions in the corresponding diagrams g¯(α, 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 kth—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 kth = 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 ki with respect to the SG phase (as evidence by a decreased g¯ 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 Ii 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, kth = 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, t0) [defined in Equation (8)] and the degree histogram N(k, t0) = Np(k, t0), for a particular time, t0 = 5 · 106MCS, 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, t0) of the recovered pattern µ is larger for highdegree nodes (**Figure 5D**), indicating that they contribute most to the overlap mµ(t0). On the contrary, for the non-recovered patterns ν, mν (k, t0) 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

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.

#### 3.4. 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 ds(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 ds(t), Ss(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 distribution p¯s → 1/Ns , where Ns = 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 realizations p¯µ 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.

# 4. 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, Tc → ∞ 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., 2007, 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 Ii

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.

### REFERENCES

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Millán, Torres and Marro. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.