Modelling the brain as an Apollonian network Gian Luca Pellegrini1, Lucilla de Arcangelis2, Hans J. Herrmann3 and Carla Perrone-Capano4 1 Department of Physical Sciences, University of Naples Federico II, 80125 Napoli, Italy 2 Dept. of Information Engineering and CNISM, Second University of Naples, 81031 Aversa (CE), Italy 3 Computational Physics, IfB, ETH-H¨onggerberg, Schafmattstr. 6, 8093 Zu¨rich, Switzerland 4 Dept. of Biological Sciences, University of Naples ”Federico II”, 80134, 7 Naples, Italy and IGB ”A.Buzzati Traverso”, CNR, 80131 Naples, Italy. 0 0 Networksoflivingneuronsexhibitanavalanchemodeofactivity,experimentallyfoundinorgan- 2 otypic cultures. Moreover, experimental studies of morphology indicates that neurons develop a network of small-world-like connections, with thepossibility of veryhigh connectivity degree. Here n we study a recent model based on self-organized criticality, which consists in an electrical network a withthresholdfiringandactivity-dependentsynapsestrengths. Westudythemodelonascale-free J network, the Apollonian network, which presents many features of neuronal systems. The system 7 exhibitsapowerlawdistributedavalancheactivity. Theanalysisofthepowerspectraoftheelectri- 2 cal signal reproducesvery robustly thepower law behaviour with theexponent0.8, experimentally measured in electroencephalograms (EEG) spectra. The exponents are found quite stable with re- ] C specttoinitialconfigurationsandstrengthofplasticremodelling,indicatingthatuniversalityholds for a wide class of brain models. N . PACSnumbers: 87.19.La,05.65.+b,05.45.Tp,89.75.-k o i b - I. INTRODUCTION a power law behaviour as function of the duration time q normalised by the binning time with an exponent equal [ to -2.0 followed by an exponential cutoff. These results 1 Neuronal networks exhibit diverse patters of activity, have been interpreted relating spontaneous activity in a v including oscillations, synchronization and waves. Dur- cortical network to a critical branching process [3], in- 5 ing neuronal activity, each neuron can receive inputs deed the experimentalbranching parameter is very close 4 by thousands of other neurons and, when it reaches to the critical value equal to one, at which avalanchesat 0 a threshold, redistributes this integrated activity back all scales exist. 1 to the neuronal network. Recently a neuronal activity 0 On the other hand, the dynamics observed in sponta- 7 based on avalanches has been observed in organotypic neousbrainactivityis verysimilar to self-organizedcrit- 0 cultures from coronal slices of rat cortex [1] where neu- icality (SOC) [4, 5, 6, 7]. The term SOC usually refers / ronal avalanches are stable for many hours [2]. More o toamechanismofslowenergyaccumulationandfasten- precisely, recording spontaneous local potentials contin- bi uously by a multielectrode array, has shown that activ- ergy redistribution driving the system toward a critical - state,wherethedistributionofavalanchesizesisapower q ity initiated at one electrode might spread to other elec- law obtained without fine tuning: no tunable parameter : trodesnotnecessarilycontiguous,asinawave-likeprop- v ispresentinthemodel. Thesimplicityofthemechanism agation. Cortical slices are then found to exhibit a new i at the basis of SOC has suggested that many physical X form of activity, producing several avalanches per hour and biological phenomena characterized by power laws r of different duration, in which non-synchronous activity a is spreadoverspace andtime. By analysingthe size and in the size distribution, represent natural realizations of SOC. For instance, SOC has been proposed to model duration of neuronal avalanches, the probability distri- earthquakes[8,9,10],theevolutionofbiologicalsystems bution reveals a power law behaviour, suggesting that [11], solar flare occurrence [12], fluctuations in confined the cortical network operates in a critical state. The plasma [13], snow avalanches[14] and rain fall [15]. experimental data indicate for the avalanche size distri- bution a slope varying between -1.2 and -1.9, depending Moreover, power law behaviour is observed in power on the accuracy of the time-binning procedure, with a spectra of different time series monitoring neural activi- value -1.5 for optimal experimental conditions. Inter- ties. Prominent examples are EEG data which are used estingly, the power law behaviour is destroyed when the byneurologiststodiscernsleepphases,diagnoseepilepsy excitability of the system is increased, contrary to what and other seizure disorders as well as brain damage and expected since the incidence of large avalanches should disease[16,17]. Anotherexampleofaphysiologicalfunc- decrease the power law exponent. The distribution then tionwhichcanbemonitoredbytimeseriesanalysisisthe becomes bimodal, i.e. dominated either by very small or human gait which is controlled by the brain [18]. For all very large avalanches as in epileptic tissue. The power these time series the power spectrum, i.e. the square of law behaviour is therefore the indication of an optimal theamplitudeoftheFouriertransformationdoubleloga- excitability in the system spontaneous activity. More- rithmically plotted against frequency, generally features over the avalanche time duration is also found to follow apowerlawatleastoveroneortwoordersofmagnitude 2 with exponents between 1 and 0.7. Moreover, experi- exponentsisfoundconsideringleakyneuronsorintroduc- mental results show that the neurotransmitter secretion ing a small precentage of inhibitory synapses, indicating rateexhibitsfluctuationswithtimehavingpowerlawbe- that universality holds for a wide class of brain models. haviour [19] and power laws are observed in fluctuations In real brain neurons are known to be able to develop of extended excitable systems driven by stochastic fluc- anextremelyhighnumberofconnectionswithotherneu- tuations [20]. rons, that is a single cell body may receive inputs from On the basis of these observations, recently a model even a hundred thousand presynaptic neurons. One of based on SOC ideas and taking into account synaptic the most fascinating questions is how an ensemble of plasticity in a neural network [21] has been proposed. living neurons self-organizes, developing connections to Plasticityisoneofthemostastonishingpropertiesofthe give origin to a highly complex system. The dynam- brain, occuring mostly during development and learning ics underlying this process should be driven both by the [22, 23, 24], and can be defined as the ability to mod- aim of realizing a well connected network leading to ef- ify the structural and functional properties of synapses, ficient information transmission, and the energetic cost properties which are thought to underlie memory and of extablishing very long connections. The morphologi- learning. Among the postulated mechanisms of synap- calcharacterizationofaneuronalnetworkgrownin vitro tic plasticity, the activity dependent Hebbian plastic- has been studied [28] by monitoring the development of ity constitutes the most fully developed and influential neurites in an ensemble of few hundred neurons from model of how information is stored in neural circuits the frontal ganglion of adult locusts. After few days [25, 26, 27]. Within a SOC approach the four most im- the cultured neurons have developed an elaborated net- portantingredientsforneuronalactivityhavebeenintro- work with hundreds of connections, whose morphology duced,namelythresholdfiring,neuronrefractoryperiod, and topology has been analysed by mapping it onto a activity-dependent synaptic plasticity and pruning. connected graph. The short path length and the high clusteringcoefficientmeasuredindicate thatthe network The system consists in an electrical network on a belongs to the category of small-world networks [29], in- squarelattice,onwhicheachsiterepresentsthecellbody terpolatingbetweenregularandrandomnetworks. How- of a neuron, each bond a synapse. Therefore, each site ever, the system grown in vitro necessarily lacks some is characterized by a potential and each bond by a con- features of in vivo systems, therefore the average node ductance. Whenever at a given time the value of the connectivity is found equal only to few units and the potential at a site is above a certain threshold, approx- imately equal to −55mV for the real brain, the neuron ”scale-free” feature [30] of many real networks was not recovered. Small-worldnetworksarecharacterizedby an fires, i.e. generates an ”action potential”, distributing efficient information transmission with a small number charges to its connected neighbours in proportion to the of long range connections. The activity dependent brain currentflowingthrougheachbond. Afterfiring,aneuron goesbacktotherestingpotentialof−70mV andremains model [21] has been implemented on small world net- works, by rewiring a small percentage of the square lat- inactive during the refractory period, when it is unable tice bonds. Again universal scaling behaviour is recov- to send or receive information from other neurons. This ered for both the avalanche distribution and the power time corresponds for real neurons to the physiological spectra. The simple rewiring procedure, however, only time needed to reset ion channels after the transmission allows long range connections leaving the average node of the action potential through the axon. The conduc- connectivity equalto a few units, as forin vitro systems. tances, on the other hand, represent Hebbian synapses, for which the conjunction of activity at the presynap- Inthispaperweinvestigatethebehaviouroftheactiv- tic and postsynaptic neuron modulates the efficiency of ity dependent brainmodel onscale-free networks,whose the synapse [27]. To this extent, each time a synapse featureareclosertothemorphologyofneuronalnetworks transmitsanactionpotentialbetweenactiveneurons,its inlivingbrains. Scale-freenetworksareindeedcharacter- strength is increased proportionally to the intensity of izedbyapowerlawdistributionofthenodeconnectivity, the transmitted signal,whereassynapses inactive during allowing a high number of connections per neuron. We a neuronal avalanche have their strength decreased, as develop the model on the Apollonian network [31], that for Hebbian rules. Synapses successively weakened may hasthepropertyofbeingsimultaneouslysmall-worldand have their strength finally set to zero, i.e. are ”pruned”, scale-free and therefore exhibits all characteristics found eliminating that particular connection between neurons. for neuronal networks. Analogously to previous stud- Pruning implies that, as activity goes on, the initial reg- ies, we analyse the behaviour of the avalanche size and ular lattice is transformed, some patterns are strength- duration distribution and the power spectra related to ened and the connectivity of some neurons decreased. electrical activity. We also study a system composed by Thesystemexhibitsanavalancheactivitypowerlawdis- both excitatory and inhibitory sysnapses,to be closer to tributed with an exponent close to -1.5, as measured for real brains. The paper is organized as follows: In sect.II spontaneousactivity[1]. The analysisofthe powerspec- the scale-free Apollonian network is described, whereas tra of the electrical signal reproduces very robustly the in section III the activity dependent brain model is pre- power law behaviour with the exponent 0.8, experimen- sented and the results on brain activity are discussed in tally measured in EEG spectra. The same value of both section IV. Concluding remarks are given in section V. 3 has small-world features. This implies [29] that the av- eragelength of the shortestpath l behaves as in random networks and grows slower than any positive power of N, i.e. l ∝ (lnN)3/4. Furthermore the clustering co- efficient C is very high as in regular networks (C = 1) and contrary to random networks. For the Apollonian network C has been found to be equal to 0.828 in the limit of large N. On this basis the Apollonian network appears to have all the features typical of neuronal net- works: small-world property found experimentally [28] and possibility of very high connectivity degree (scale- free). Moreover it also presents bonds connecting sites of all lengths. Also this last feature is characteristic of neuronalnetworksinbraincortex,wherethelengthofan axonconnecting the pre-synaptic with the post-synaptic neuron can vary over several orders of magnitude, from µm to cm. FIG.1: Apollonian network for N =2: iterations n=0,1,2 are symbols (cid:13),(cid:4),•, respectively. III. ACTIVITY DEPENDENT MODEL II. APOLLONIAN NETWORK On a Apollonian network at generation N, we assign The Apollonian network has been recently introduced at each site a neuron at potential vi and at each bond a [31] in a simple deterministic version starting from the synapseofconductancegij. Wheneverattimetthevalue problem of space-filling packing of spheres according to of the potential at a site i is above a certain threshold theancientGreekmathematicianApolloniusofPerga. In vi ≥vmax, the neuron generates an action potential, dis- itsclassicalversionthenetworkassociatedtothepacking tributing charges to connected neurons in proportion to gives a triangulation that physically corresponds to the the current flowing through each bond forcenetworkofthepacking. Onestartswiththezero-th i (t) ij order triangle of corners P1,P2,P3, places a fourth site vj(t+1)=vj(t)+vi(t)P i (t) (1) P in the center of the triangle and connects it to the k ik 4 three corners (n = 0). This operation will divide the where v (t) is the potential at time t of site j, connected j originaltriangle in three smaller ones having in common to site i,i =g (v −v ) and the sum is extended to all ij ij i j the central site. The iteration n = 1 proceeds placing k sitesconnectedtosite ithatareatapotentialv <v . k i one more site in the center of each small triangle and After firing a neuron is set to a zero resting potential. connecting it to the corners (Fig.1). At each iterationn, The conductances can be initially all set equal or else goingfrom 0 to N, then the number ofsites increasesby random between 0 and 1, whereas the neuron poten- a factor 3 and the coordination of each already existing tials areuniformly distributed randomnumbers between site by a factor 2. More precisely, at generation N there v −2 and v −1. In agreement with the SOC sce- max max are nario,theinitialstateforthevoltageisnotrelevantsince the systemevolvestowardthe samecriticalstateregard- m(k,N)=3N,3N−1,3N−2,...,32,3,1,3 lessoftheinitialcondition. Thepotentialisfixedtozero at the three sites 1,2,3 where information can flow out vertices with connectivity degree of the system. The external stimulus can be imposed at k(N)=3,3×2,3×22,...,3×2N−1,3×2N,2N+1+1 one input site choseneither fixed or at random, this last case modelling more closely spontaneous brain electrical respectively where the two last values correspond to the activity. site P4 and the three corners P1,P2,P3. The maximum The firing rate of real neurons is limited by the re- connectivityvalue thenis the one ofthe verycentralsite fractory period, i.e. the brief period after the generation P4,kmax =3×2N,whereasthesitesinsertedattheN-th of an action potential during which a second action po- iteration will have lowest connectivity 3. tential is difficult or impossible to elicit. The practical The important property of the Apollonian network is implicationofrefractoryperiodsisthattheactionpoten- that it is scale-free. In fact, it has been shown [31] tial does not propagate back toward the initiation point that the cumulative distribution of connectivity degree and therefore is not allowed to reverberate between the P(k)=Pk′≥km(k,N)/NN, where NN = 3+(3(N+1)− cell body and the synapse. In the model, once a neu- 1)/2 is the total number of sites at generation N, has ron fires, it remains quiescent for one time step and is a power law behaviour with k. More precisely, P(k) ∝ therefore unable to accept charge from firing connected k1−γ, with γ =ln3/ln2∼1.585. Moreover the network neurons. This ingredient indeed turns out to be crucial 4 for a controlled functioning of the numerical model. In this way an avalanche of charges can propagate far from the input through the system. As soon as a site is at or above threshold v at a max given time t, it fires according to Eq. (1). Then the conductance of all the bonds, connecting to active neu- rons and that have carried a current, is increased in the following way g (t+1)=g (t)+δg (t) (2) ij ij ij where δg (t) = Aαi (t), with α being a dimensionless ij ij parameter and A a unit constant bearing the dimen- sion of an inverse potential. After applying Eq. (2) the time variable of the simulation is increased by one unit. Eq. (2) describes the mechanism of increase of synap- tic strength, tuned by the parameter α. This parameter then representsthe ensembleofallpossible physiological factorsinfluencingsynapticplasticity,manyofwhichare FIG.2: Averagenumberofprunedbondsasfunctionoftime for threedifferentvaluesof αand equalinitial conductances. not yet fully understood. In the inset, asymptotic number of active bonds as function Once anavalancheoffiringscomesto anend, the con- of α. The maximum is for α=0.020 where are active about ductance of all the bonds with non-zero conductance is 80% of bonds. reduced by the average conductance increase per bond, ∆g =P δg (t)/N , where N is the number of bonds ij,t ij b b withnon-zeroconductance. Thequantity∆gdependson sites of the system have always zero potential and rep- αandontheresponseofthebraintoagivenstimulus. In resent open boundary conditions. The input site is ei- this way the network ”memorizes” the most used paths ther chosenat randomorfixed. Itis worthnoticing that ofdischargebyincreasingtheirconductance,whereasthe the case of random input sites simulates more closely less used synapses atrophy. Once the conductance of a the spontaneous activity of the system. Synapses can bond is below an assigned small value σ , it is removed, t be excitatory or inhibitory with probability p . Ini- i.e. is set equal to zero, which corresponds to what is inh tial conductances can be either all equal to g = 0.25 or knownaspruning. Thisremodellingofsynapsesmimicks 0 randomly distributed between 0 and 1. The other pa- thefinetuningofwiringthatoccursduring”criticalperi- rametersinthe simulationare: firingthresholdv =6 ods”inthe developingbrain,whenneuronalactivitycan max and conductance cut-off for pruning σ = 0.0001. Their modify the synaptic circuitry, once the basic patterns of t value does not influence the simulation results. brainwiringareestablished[23]. Thesemechanismscor- respondtoaHebbianformofactivitydependentplastic- Pruning ity, where the conjunction of activity at the presynaptic The strength of the parameter α, controlling both the and postsynaptic neuron modulates the efficiency of the increase and decrease of synaptic strength, determines synapse [27]. To insure the stable functioning of neural theplasticitydynamicsinthe network. Infact,themore circuits,both strengtheningandweakeningrulesofHeb- the system learns strengthening the used synapses, the bian synapses are necessary to avoid instabilities due to more the unused connections will weaken. We apply a positivefeedback[32]. However,differentlyfromthewell sequence of external stimuli and we measure the total known Long Term Potentiation (LTP) and Long Term number of pruned bonds at the end of each avalanche, Depression(LTD)mechanisms,themodulationofsynap- N . This quantity in general could depend on the ini- pb ticstrengthdoesnotdependonthefrequencyofsynapse tial conductance g , therefore the two cases of all initial 0 activation [22, 33, 34]. The external driving mechanism conductances equal to 0.25, and uniformly distributed to the system is imposed by setting the potential of the between 0 and 1, are investigated. input site tothe valuevmax, correspondingto one stimu- First the case of equal initial conductances is anal- lus. This external stimulus is needed to keep the system ysed. For each value of α the average number of pruned functioning andthereforemimicks the livingbrainactiv- bonds, N , is monitored as function of time, where a pb ity. The discharge evolves until no further firing occurs, time unit corresponds to the application of an external then the next stimulus is applied. stimulus. Forinputsitesrandomlychosenateachstimu- lus, Fig.2 shows that pruning starts after a certain time, since all conductances are initially equal to 0.25, and IV. NUMERICAL RESULTS N increases more rapidly with time for larger α. The pb plateau is reached after about 5000 stimuli (for everyα) WeconsideranApolloniannetatthegenerationN =9 afterwhichN increasesonlyoffewunitsintime. From pb (29527 neurons and 177150 synapses). The three corner the asymptotic value of each curve we can evaluate the 5 FIG.3: (Coloronline)Probabilityofpruningforbondsofdif- FIG. 4: (Color online) Connectivity degree distribution n(k) ferentiterationsnasfunctionoftimeforequalinitialconduc- atdifferentpruningstagestforequalinitialconductancesand tances. In the inset, the asymptotic Npb (after 5000 stimuli) α = 0.020. In the inset the corresponding behaviour of the as function of nwith the exponentialfit Npb ≃exp1.2n. numberof prunedbonds. asymptotic number of active bonds as function of α and network (Fig.4) has also been analysed, i.e. the number determine that the value of α maximizing the number of of sites with a given connectivity degree k as function of activebondsisabout0.020. Thiscouldbeinterpretedas k in the initial network and after application of a given an optimal value for the system with respect to plastic number of external stimuli. In order to identify the dif- adaptation: it maximizes the number of active connec- ferent stages in the pruning process, the inset of Fig.4 tions under the competing strengthening and weakening shows the total number of pruned bonds as function of rules. time. After the application of few external stimuli, i.e. for a short plastic training, the distribution n(k) shows Inordertounderstandifpruningactsinthesameway the same scaling behaviour of the Apollonian network. on bonds created at different iterations n,n = 0,...,N, As the pruning process goes on, sites vary their connec- or rather tends to eliminate some particular iteration, tivity and new values of k appear. The result is that the probability to prune bonds of different n is evalu- the scaling behaviour is progressivelylost, as well as the ated, that is the number of pruned bonds over the total scale-free character of the network, since there is a gen- number of bonds for each iteration stage, as function of eralized decrease of connectivity in the network. In the time. Fig.3 shows that the plateau is reached at about analysisof spontaneousactivity itis therefore important the same time and the shape of the curve is similar for to impose a not too extended plastic training in orderto each n. However the probability to prune bonds with avoid an excessive decrease of connectivity degree of the largenis higher: These arethe bondscreatedinthe last network. iterations and therefore embedded in the interior of the network. This suggests that the most active bonds are Spontaneous activity: avalanche distributions thelongrangeones(smalln),thatthereforeoptimizein- Aftertrainingthesystemapplyingplasticityrulesdur- formationtransportthroughthenetwork. Inthe insetof ing N external stimuli, we now submit the system to a p Fig.3 we show the asymptotic number of pruned bonds new sequence of stimuli with no modification of synap- per generation on a semi-log scale, this quantity is well sis strength. The response of the system to this second fitted by the exponential behaviour Npb ≃expn. sequence models the brain spontaneous activity, which The sameanalysishasbeenperformedforrandomini- is analysed by measuring the avalanche size distribution tial conductances between 0 and 1. The results are sim- n(s),thetime durationdistributionn(T),andthepower ilar to the previous case. It can be noticed that pruning spectrum S(f). starts already at t = 1, since conductances close to zero The avalanche size distribution n(s) consistently ex- are present, and the plateau is reached after about 3000 hibits power law behaviour for different values of model stimuli. The value ofαwhichnowoptimizesthe number parameters. Fig.5 shows the avalanche size distribution of active bonds is about 0.030. Finally the pruning be- for different values of N , including also the case N =0 p p haviour for different iterations is similar to the previous (no plasticity training) for random initial conductances. case,with the pruning probability alsoincreasingwith n We notice that, for fixed size s, increasing N decreases p exponentially as Npb ≃expn. the number of avalanches of that size, suggesting that Theeffectofpruningontheconnectivitydegreeofthe strong plasticity remodelling decreases activity. The ex- 6 FIG.5: (Coloronline)Avalanchesizedistributionfordifferent FIG. 6: (Color online) Avalanche size distribution for input values of Np, random initial conductances, α = 0.030 and sitesrandomlychosenamongsiteswiththesameconnectivity random input site. In the inset the corresponding behaviour degreek. Onlydistributionsforsmallk areshown,forhigher of the numberof prunedbonds. k the scaling behaviour is lost (random initial conductances, α=0.030,Np =100). fore the network activity will be damped already at the initial site. Therefore, to reproduce the experimentally ponent appears to be independent of N as long as the p observed scaling behaviour, the fixed input site should numberofprunedbonds,N ,isfarfromtheplateau(see pb be chosen with low connectivity degree (k ≤ 12). The inset in Fig.5). Similar results are found for equal initial avalanche size distribution for fixed input site with con- conductances,The value ofthe exponentis σ =1.8±0.2 nectivity k = 3 or k = 6 exhibits power law behaviour and is stable with respect to variations of the parameter with the same exponents found for random input site: α for both equal and random initial conductance. This σ =1.8±0.2 for equal and randominitial conductances. valueiscompatiblewithinerrorbarswiththevaluefound At time t = 0 a neuron is activated by an external in the experiments of Beggs and Plentz [1], 1.5 ± 0.4, stimulus initiating the avalanche. This will continue un- and with previous results of the model on both regular til no neuron is at or above threshold. The number of and small world lattices. This suggests that the high avalanches lasting a time T, n(T), as function of T ex- levelofconnectivityreducesthe probabilityofverylarge hibits power law behaviour (Fig.7) with an exponential avalanchesbutdoesnotchangesubstantiallythesponta- cutoff. The scalingexponentis found to be τ =2.1±0.2 neousactivity behaviour. For largerN ,the distribution p for equal and random initial conductances. This value exhibits an increase in the scaling exponent and finally is found to be stable with respect to different α (Fig.7) looses the scaling behaviour for very large N values in p and N , provided that the number of pruned bonds N the plateau regime for the number of pruned bonds. p pb is lower than the plateau for that value of α. Moreover It is important to investigate the role of the choice of it does not depend on the choice of the input site, either fixed input site, since in the Apollonian network, con- fixed or random. Finally both values agree within error trary to the regular network, sites may have very dif- bars with the value 2.0, exponent found experimentally ferent connectivity degree. Fig.6 shows the avalanche by Beggs and Plentz [1]. size distribution for input sites randomly chosen among sites with given connectivity degree k. In this way it is Power spectra for spontaneous activity possible to detect solely the effect due to the connectiv- In order to compare the results of the Apollonian net- ity of the input site, eliminating all other effects due to workwithEEGmedicaldata,thepowerspectrumofthe the particular position of the input site in the network. resultingtime seriescanbe calculated. Forthis purpose, Power law behaviour is found for connectivity degree of the number of activeneuronsis monitoredasfunction of the input site up to k = 12. The scaling exponent de- time during spontaneous activity. Fig.8 shows an exam- creaseswithincreasingconnectivitydegreekoftheinput ple of neuronal activity where avalanches of all sizes can site, that is for larger k larger avalanches become more be generated in response to the external stimulus. The probable. However, if the connectivity degree increases powerspectrumiscalculatedasthesquaredamplitudeof further, the scaling behaviour is lost. This is due to the the Fourier transform as function of frequency, averaged fact that an input site with very high connectivity must over many sequences of external stimuli. distribute its charge to many connected sites and there- Fig.9showsthespectrumforequalinitialconductances 7 FIG. 7: (Color online) Avalanche duration distribution for FIG. 9: (Color online) Power spectra for different Np, equal different values of α (random initial conductances, random initial conductances, α=0.020 and random input sites. input sites, Np =500). The dotted line has slope 2.1. the power spectrum for fixed input site shows a scaling exponent β = 0.8±0.1 over two orders of magnitude. The measuredvalue forthe powerspectraexponentis in agreementwiththeexpectedrelationwiththescalingex- ponent of the avalancheduration distribution β =3−τ, being −τ <−1 [5]. The scaling behaviour of the power spectrum can be interpreted in terms of a stochastic process determined by multiple random inputs [38]. In fact, the output sig- nal resulting from different and uncorrelated superim- posed processes is characterized by a power spectrum with powerlaw behaviour anda crossoverto white noise at low frequencies. The crossoverfrequency is related to the inverse of the longest characteristic time among the superimposed processes. The value of the scaling expo- nentdependsontheratiooftherelativeeffectofaprocess of given frequency on the output with respect to other processes. 1/f noise corresponds to a superposition of FIG. 8: Total current flowing in the system as function of time. Avalanchesof all sizes can beobserved. processes of different frequency having all the same rela- tive effect on the output signal. In our case the scaling exponentissmallerthanunity,suggestingthatprocesses withhighcharacteristicfrequencyaremorerelevantthan and different values of N . For N = 0, i.e. when no p p processes with low frequency in the superposition [38]. plasticity mechanism is applied, the spectrum has a be- haviour 1/f,characteristicof SOC. For values of N dif- Inhibitory synapses p ferent from zero but before N reaches the plateau, one In the mature living brain synapses can be excitatory pb can distinguish two different regimes: a power law be- or inhibitory, namely they set the potential of the post- haviour with exponent β = 0.8±0.1 at high frequency, synaptic membrane to a level closer or farther, respec- followed by a crossover toward white noise at low fre- tively, to the firing threshold. This ingredient can be quency. However, for N = 2000 (close to the plateau indroduced by considering each synapse inhibitory with p value for N ) the scaling behaviour with exponent 0.8 probability p and excitatory with probability 1−p . pb in in is detected over a wider frequency range. The difference The avalanche size and duration distributions show that between β = 1 for N = 0 and β ≃ 0.8 for higher N , the exponents σ and τ increase for increasing p , there- p p in suggeststhatthe existence ofplasticityrulesreducesthe foreforahighpercentageofinhibitorysynapsestheprob- power spectrum exponent reaching agreement with ex- abilityoflargeavalanchesdecreases(Fig.10). Onthereg- perimental EEG spectrum [36, 37]. The stability of the ular lattice for p = 0.5 no longer power law but expo- in exponent with respect to α has also been verified, find- nential behaviour is found [21]. In the present case scal- ing consistently β =0.8±0.1 at high frequency. Finally ing behaviour persists due to the very high connectivity 8 free Apollonian network. The results are compared with previous simulations on regular and small world lat- tices and with experimental data. We first find the striking result that an optimal value of of the plastic- ity strength α exists with respect to the pruning pro- cess. Moreover, it appears that synapses of later gen- erations, deeply embedded in the network, are pruned with higher probability with respect to bonds of the early generations, mostly long range, that optimize in- formationtransmission. Moreovertheavalanchesizedis- tribution shows a power law behaviour with an expo- nent σ = 1.8±0.2 for equal and random initial conduc- tances. This value is compatible with 1.5±0.4, exper- imentally found for neuronal avalanches and recovered by the model on the square lattice and small world net- works. Also the avalanche durationdistribution exhibits powerlawbehaviourwithanexponentialcutoff,inagree- FIG.10: (Coloronline)Avalanchesizedistributionsfordiffer- ment with experimental results of Beggs and Plentz [1]. ent pin, equal initial conductances, α=0.020, random input The exponenthasvalue: τ =2.1±0.2for equalandran- sites and Np =1700. dominitialconductance,inagreementwith2.0foundex- perimentally. Furthermore the power spectrum exhibits powerlawbehaviourathighfrequencywithβ =0.8±0.1, degree, suggesting that the Apollonian network is better in agreement with experimental data [36, 37]. At inter- suited to model the neural connections in the brain. mediate frequency the slope becomes greaterthan unity, The power spectra for different values of p exhibit a in crossing over to white noise at low frequencies. None of complex behaviour. In fact, for a small fraction of in- thescalingexponentsforspontaneousactivityinthecase hibitorysynapses(p ≤0.05)the powerlawexponentβ in of excitatory synapses depends on the particular choice increases with respect to the case where synapses are all for the length or strength of the plasticity training and excitatory up to a value 1.2. Then, for p ∼ 0.10, the in are quite stable with respect to the initial conductance exponentdecreasestowardvaluescompatiblewithexper- configurations. These results suggest that also on Apol- imental results, i.e. between 0.7 and 1.0. By increasing loniannetworkuniversalbehaviourfoundforregularand further the percentage of inhibitory synapses, to values small world networks [21] holds. Furthermore, the scale- close or greater than 0.2, the spectrum becomes the one freeApolloniannetprovidesanexcellentdescriptionboth of white noise. of the morphology and the electrical activity properties of the brain. V. CONCLUSIONS Acknowledgements. This work was supported by MIUR-PRIN 2004, MIUR-FIRB 2001, CRdC-AMRA Extensivesimulations havebeen performedfor the ac- and EU Network Number MRTN-CT-2003-504712. tivity dependent brain model implemented on the scale- H.J.H. acknowledges the Max Planck prize. [1] J. M. Beggs, D. Plenz, J. Neurosci. 23, 11167 (2003). [11] P. Bak, K. Sneppen,Phys.Rev. Lett. 71, 4083 (1993). [2] J. M. Beggs, D. Plenz, J. Neurosci. 24, 5216 (2004). [12] E.T. Lu, R.J. Hamilton, Astrophys.J. 380, L89 (1991). [3] S. Zapperi, K.B. Lauritsen, H.E. Stanley, Phys. Rev. [13] P. A. Politzer, Phys. Rev.Lett. 84, 1192 (2000). Lett.75, 4071 (1995). [14] J. Faillettaz, F. Louchet, J.R. Grasso, Phys. Rev. Lett. [4] P. Bak, How nature works. The science of self-organized 93, 208001 (2004). criticality, Springer, New York,1996. [15] O. Peters, C. Hertlein and K. Christensen, Phys. Rev. [5] H.J. Jensen, Self-Organized Criticality, Cambridge Uni- Lett. 88, 018701 (2002) versity Press, Cambridge, 1998. [16] A. Gevins et al, Trends Neurosci. 18, 429 (1995). [6] S. Maslov, M. Paczuski, P. Bak, Phys. Rev. Lett. 73, [17] G. Buzsaki, A. Draguhn,Science 304, 1926 (2004). 2162 (1994). [18] J. M. Hausdorff et al., Physica A. 302, 138 (2001). [7] J. Davidsen, M. Paczuski, Phys. Rev. E 66, 050101(R) [19] S.B. Lowen, S.S. Cash, M. Poo , M.C. Teich, J. Neuro- (2002). science, 17, 5666 (1997). [8] P.Bak, C.Tang, J. Geophys. Res. 94 , 15635 (1989). [20] D.R. Chialvo, G.A. Cecchi, M.O. Magnasco, Phys. Rev. [9] A.Sornette, D.Sornette, Europhys.Lett.9, 197 (1989). E 61, 5654 (2000). [10] E. Lippiello, L. de Arcangelis, C. Godano, Europhys. [21] L. de Arcangelis, C. Perrone-Capano, H.J. Herrmann, Lett.72 , 678 (2005). Phys. Rev.Lett. 96, 028107 (2006). 9 [22] T. D. Albraight et al, Neuron, Review supplement to [31] J.S. Andrade, H.J. Herrmann, R.F.S. Andrade, L.R. da vol.59 (February 2000). Silva, Phys. Rev.Lett. 94, 018702 (2005). [23] T. K. Hensch,Ann.Rev.Neurosci. 27, 549 (2004). [32] N.S. Desai, J. Physiol. Paris 97, 391 (2003). [24] L.F. Abbott, S.B. Nelson, Nature Neurosci. 3, 1178 [33] O. Paulsen, T. J. Sejnowski, Curr. Opin. Neurobiol. 10, (2000). 172 (2000). [25] D.O. Hebb, The organization of behaviour, New York: [34] K. H. Braunewell, D.Manahan-Vaughan, Rev. Neurosci. John Wiley, 1949. 12, 121 (2001). [26] J.Z.Tsien,Curr.Opi.nNeurobiol. 10,266(2000); G.-Q. [35] L. F. Lago-Fernandez, R. Huerta,F. Corbacho,J.A. Bi, M.-M. Poo, Ann.Rev.Neurosci. 24, 139 (2001). Siguenza, Phys.Rev.Lett. 84, 2758 (2000). [27] S.J. Cooper, Neurosci. Biobehav. Rev. 28, 851 (2005). [36] W. J. Freeman et al, J. Neurosci. Meth. 95, 111 (2000). [28] O. Shefi, I. Golding, R. Segev, E. Ben-Jacob, A. Ayali, [37] E. Novikov, A. Novikov, D. Shannahoff-Khalsa, B. Phys.Rev.E. 66, 021905 (2002). Schwartz, J. Wright,Phys. Rev.E. 56, R2387 (1997). [29] D.J. Watts, S.H. Strogatz, Nature 393, 440 (1998). [38] J.M. Hausdorff, C.K. Peng, Phys. Rev. E. 54, 2154 [30] L.A.N. Amaral,A. Scala, M. Barthelemy, E.H. Stanley, (1996). Proc. Natl. Acad.Sci. U.S.A., 97, 149 (2000).