Self-organization using synaptic plasticity Vicenc Gomez1 ¸´ vgomez@iua.upf.edu Hilbert J Kappen1 b.kappen@science.ru.nl 1 Department of Biophysics Radboud University Nijmegen 6525 EZ Nijmegen, The Netherlands Andreas Kaltenbrunner2 andreas.kaltenbrunner@upf.edu ´ Vicente Lopez2 vicente.lopez@barcelonamedia.org 2 Barcelona Media - Innovation Centre Av. Diagonal 177, 08018 Barcelona, Spain Abstract Large networks of spiking neurons show abrupt changes in their collective dynamics resembling phase transitions studied in statistical physics. An example of this phenomenon is the transition from irregular, noise-driven dynamics to regular, self-sustained behavior observed in networks of integrate-and-fire neurons as the interaction strength between the neurons increases. In this work we show how a network of spiking neurons is able to self-organize towards a critical state for which the range of possible inter-spike-intervals (dynamic range) is maximized. Self-organization occurs via synaptic dynamics that we analytically derive. The resulting plasticity rule is defined locally so that global homeostasis near the critical state is achieved by local regulation of individual synapses. 1 Introduction It is accepted that neural activity self-regulates to prevent neural circuits from becoming hyper- or hypoactive by means of homeostatic processes [14]. Closely related to this idea is the claim that optimal information processing in complex systems is attained at a critical point, near a transition between an ordered and an unordered regime of dynamics [3, 11, 9]. Recently, Kinouchi and Copelli [8] provided a realization of this claim, showing that sensitivity and dynamic range of a network are maximized at the critical point of a non-equilibrium phase transition. Their findings may explain how sensitivity over high dynamic ranges is achieved by living organisms. Self-Organized Criticality (SOC) [1] has been proposed as a mechanism for neural systems which evolve naturally to a critical state without any tuning of external parameters. In a critical state, typical macroscopic quantities present structural or temporal scale-invariance. Experimental results [2] show the presence of neuronal avalanches of scale-free distributed sizes and durations, thus giving evidence of SOC under suitable conditions. A possible regulation mechanism may be provided by synaptic plasticity, as proposed in [10], where synaptic depression is shown to cause the mean synaptic strengths to approach a critical value for a range of interaction parameters which grows with the system size. In this work we analytically derive a local synaptic rule that can drive and maintain a neural network near the critical state. According to the proposed rule, synapses are either strengthened or weakened whenever a post-synaptic neuron receives either more or less input from the population than the required to fire at its natural frequency. This simple principle is enough for the network to selforganize at a critical region where the dynamic range is maximized. We illustrate this using a model of non-leaky spiking neurons with delayed coupling for which a phase transition was analyzed in [7]. 1 2 The model The model under consideration was introduced in [12] and can be considered as an extension of [15, 5]. The state of a neuron i at time t is encoded by its activation level ai (t), which performs at discrete timesteps a random walk with positive drift towards an absorbing barrier L. This spontaneous evolution is modelled using a Bernoulli process with parameter p. When the threshold L is reached, the states of the other units j in the network are increased after one timestep by the synaptic efficacy j i , ai is reset to 1, and the unit i remains insensitive to incoming spikes during the following timestep. The evolution of a neuron i can be described by the following recursive rules: jN a (t) + ij HL (aj (t)) + 1 with probability p i =1,j =i ai (t + 1) = if ai (t) < L jN ai (t) + ij HL (aj (t)) with probability 1 - p =1,j =i N ai (t + 1) = 1 + j ij HL (aj (t)) if ai (t) L (1) =1,j =i where HL (x) is the Heaviside step function: HL (x) = 1 if x L, and 0 otherwise. NN Using the mean synaptic efficacy: = i j,j =i ij /(N (N - 1)) we describe the degree of interaction between the units with the following characteristic parameter: = L-1 , (N - 1) (2) which indicates whether the spontaneous dynamics ( > 1) or the message interchange mechanism ( 1) dominates the behavior of the system. As illustrated in the right raster-plot of Figure 1, at > 1 neurons fire irregularly as independent oscillators, whereas at = 1 (central raster-plot) they synchronize into several phase-locked clusters. The lower , the less clusters can be observed. For = 0.5 the network is fully synchronized (left raster-plot). In [7] it is shown that the system undergoes a phase transition around the critical value = 1. The study provides upper (max ) and lower bounds (min ) for the mean inter-spike-interval (ISI) of the ensemble and shows that the range of possible ISIs taking the average network behavior ( = max -min ) is maximized at = 1. This is illustrated in Figure 1 and has been observed as well in [8] for a similar neural model. The average of the mean ISI is of order N x with exponent x = 1 for > 1, x = 0.5 for = 1, and x = 0 for < 1 as N , and can be approximated as shown in [7] with 1 : L 2 L-1-N N -1-N app = 1 + + +1 + . (3) 2p 2p 2p 3 Self-organization using synaptic plasticity We now introduce synaptic dynamics in the model. We first present the dissipated spontaneous evolution, a magnitude also maximized at = 1. The gradient of this magnitude turns to be simple analytically and leads to a plasticity rule that can be expressed using only local information encoded in the post-synaptic unit. 3.1 The dissipated spontaneous evolution During one ISI, we distinguish between the spontaneous evolution carried out by a unit and the actual spontaneous evolution needed for a unit to reach the threshold L. The difference of both quantities can be regarded as a surplus of spontaneous evolution, which is dissipated during an ISI. 1 The equation was denoted min in [7]. We slightly modified it using and replacing by Eq. (2). 2 60 # neuron 0 25 50 50 0 # neuron time 10 20 # neuron 0 time 50 100 40 = max - min 1 clustering 30 = 0.5 full synchronization 20 >1 noisy firing 10 0 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 Figure 1: Number of possible ISIs according to the bound = max - min derived in [7]. For > 1 the network presents sub-critical behavior and is dominated by the noise. For < 1 it shows super-critical behavior. Criticality is produced at = 1, which coincides to the onset of sustained activity. At this point, the network is also broken down in a maximal number of clusters of units which fire according to a periodic pattern. Figure 2a shows an example trajectory of a neuron's state. First, we calculate the spontaneous evolution of the given unit during one ISI, which it is just its number of stochastic state transitions during an ISI of length (thick black lines in Figure 2a). These state transitions occur with probability p at every timestep except from the timestep directly after spiking. Using the average ISI-length o ver many spikes and all units we can calculate the average total spontaneous evolution: Etotal = ( - 1)p. (4) Since the state of a given unit can exceed the threshold because of the received messages from the rest of the population (blue dashed lines in Figure 2a), a fraction of (4) is actually not required to induce a spike in that unit, and therefore is dissipated. We can obtain this fraction by subtracting from (4) the actual number of state transitions that was required to reach the threshold L. The latter quantity can be referred to as effective spontaneous evolution Eef f and is on average L - 1 minus (N - 1) , the mean evolution caused by the messages received from the rest of the units during an ISI. For 1, the activity is self-sustained and the messages from other units are enough to drive a unit above the threshold. In this case, all the spontaneous evolution is dissipated and Eef f = 0. Summarizing, we have that: L - 1 - (N - 1) for 1 Eef f = max{0, L - 1 - (N - 1) } = (5) 0 for < 1 If we subtract (5) from Etotal (4), we obtain the mean dissipated spontaneous evolution, which is visualized as red dimensioning in Figure 2a: Ediss = Etotal - Eef f = ( - 1)p - max{0, L - 1 - (N - 1) }. (6) Using (3) as an approximation of we can get an analytic expression for Ediss . Figures 2b and c show this analytic curve Ediss in function of together with the outcome of simulations. At > 1 the units reach the threshold L mainly because of their spontaneous evolution. Hence, Etotal Eef f and Ediss 0. The difference between Etotal and Eef f increases as approaches 1 because the message interchange progressively dominates the dynamics. At < 1, we have Eef f = 0. In this scenario Ediss = Etotal , is mainly determined by the ISI and thus decays again for decreasing . The maximum can be found at = 1. 3.2 Synaptic dynamics After having presented our magnitude of interest we now derive a plasticity rule in the model. Our approach essentially assumes that updates of the individual synapses ij are made in the direction of 3 L 11 threshold Ediss Spontaneous evolution Messages from other units Surplus of spontaneous evolution (a) spike (b) 25 20 15 10 5 0 0.5 1 1.5 Sim. E diss 9 a1(t) 7 13 (-1)p E total E eff L-1 (c) 100 80 Etotal E eff 5 12 60 (N-1) 12+13 Ediss 3 40 20 0 0.8 0.9 1 1.1 1.2 1 0 2 4 6 8 10 12 t 14 16 Figure 2: (a) Example trajectory of the state of a neuron: the dissipated spontaneous evolution Ediss is the difference between the total spontaneous evolution Etotal (thick black lines) and the actual evolution required to reach the threshold Eef f (dark gray dimensioning) in one ISI. (b) Ediss is maximized at the critical point. (c) The three different evolutions involved in the analysis (parameters for (b) and (c) are N = L = 1000 and p = 0.9. For the mean ISI we used app of Eq. (3)). the gradient of Ediss . The analytical results are rather simple and allow a clear interpretation of the underlying mechanism governing the dynamics of the network under the proposed synaptic rule. We start approximating the terms N and (N - 1) by the sum of all pre-synaptic efficacies ik : N = (N - 1) + (N - 1) = iN k =1 ik /N =i An update of ij occurs when a spike from the pre-synaptic unit j induces a spike in a post-synaptic unit i. Other schedulings are also possible. The results are robust as long as synaptic updates are produced at the spike-time of the post-synaptic neuron. , i i i Eef f Etotal Ediss = - (9) ij = ij ij ij where the constant scales the amount of change in the synapse. We can write the gradient as: L L-1-P k 0 k=i ik 1 1 if (L - 1 - )<0 -2 +2 i k=i ik 2p Ediss 1 )=0 = - - indef if (L - 1 - P 2P k=i ik ij 2 -1- k=i ik -1 k=i ik if (L - 1 - +1 + =i ik ) > 0. 2p 2p This can be done for large N and if we suppose that the distribution of ik is the same for all i. Ediss is now defined in terms of each individual neuron i as: L k k k 2 ik -1- L-1- ik ik =i =i =i i p + +1 + Ediss = 2p 2p 2p k - max{0, L - 1 - ik }. (8) =i k ik . =i (7) (10) For a plasticity rule to be biologically plausible it must be local, so only information encoded in the states of the pre-synaptic j and the post-synaptic i neurons must be considered to update ij . 4 (a) 1 dEtotal dEtotal+1 0.5 dEdiss (b) 0.5 c = 0.05 c = 0.5 c=5 0.25 0 ij -250 0 Li 250 500 0 -0.5 -0.25 -1 -500 -0.5 -500 -250 0 Li 250 500 Figure 3: Plasticity rule. (a) First derivative of the dissipated spontaneous evolution Ediss for = 1, L = 1000 and c = 0.9. (b) The same rule for different values of c. k We propagate =i ik to the state of the post-synaptic unit i by considering for every unit i, an effective threshold Li which decreases deterministically every time an incoming pulse is received [6]. k At the end of an ISI Li (L - 1 - =i ik ) and encodes implicitly all pre-synaptic efficacies of i. Intuitively, Li indicates how the activity received from the population in the last ISI differs from the activity required to induce and spike in i. The only term involving non-local information in (10) is the noise rate p. We replace it by a constant c and show later its limited influence on the synaptic rule. With these modifications we can write i the derivative of Ediss with respect to ij as a function of only local terms: i Ediss = ij -Li - c 2 ( Li + 2c) + 2c(L - 2 + Li ) sgn(Li ) 2 (11) Note that, although the derivation based on the surplus spontaneous evolution (10) may involve information not locally accessible to the neuron, the derived rule (11) only requires a mechanism to keep track of the difference between the natural ISI and the actual one. We can understand the mechanism involved in a particular synaptic update by analyzing in detail Eq. (11). In the case of a negative effective threshold (Li < 0) unit i receives more input from the rest of the units than the required to spike, which translates into a weakening of the synapse. Conversely, if Li > 0 some spontaneous evolution was required for the unit i to fire, Eq. (11) is positive and the synapse is strengthened. The intermediate case (Li = 0), corresponds to = 1 and no synaptic update is needed (nor is it defined). We will consider it thus 0 for practical purposes. i Figure 3a shows Eq. (11) in bold lines together with Etotal / ij (dashed line, corresponding to i < 1) and Etotal / ij + 1 (dashed-dotted, > 1), for different values of the effective threshold Li of a given unit at the end on an ISI. Etotal indicates the amount of synaptic change and Eef f determines whether the synapse is strengthened or weakened. The largest updates occur in the transition from a positive to a negative Li and tend to zero for larger absolute values of Li . Therefore, significant updates correspond to those synapses with post-synaptic neurons which during the last ISI have received a similar amount of activity from the whole network as the one required to fire. We remark the similarity between Figure 3b and the rule characterizing spike time dependent plasticity (STDP) [4, 13]. Although in STDP the change in the synaptic conductances is determined by the relative spike timing of the pre-synaptic neuron and the post-synaptic neuron and here it is determined by Li at the spiking time of the post-synaptic unit i, the largest changes in STDP occur also in an abrupt transition from strengthening to weakening corresponding to Li = 0 in Figure 3a. Figure 3b illustrates the role of c in the plasticity rule. For small c, updates are only significant in a tiny range of Li values near zero. For higher values of c, the interval of relevant updates is widened. The shape of the rule, however, is preserved, and the role of c is just to scale the change in the synapse. For the rest of this manuscript, we will use c = 1. 5 = 0.1 1.8 ' 1.6 ' 1.4 1.2 ' 1 1 200 400 600 0.98 1.02 ' 1 1 = 0.1 ' 0.8 ' ' 0.98 1.02 = 0.01 0.6 0 1000 # periods 1 0.98 0 100 200 # periods 300 0.6 0 0.8 1 0 1 0 1.02 1 = 0.1 1.02 = 0.01 1.2 ' 0.98 1.4 1.8 1.6 = 0.01 1000 2000 3000 2000 1 # periods 2 x 10 4 Figure 4: Empirical results of convergence toward = 1 for three different initial states above (top four plots) and below (bottom four plots) the critical point. Horizontal axis denote number of ISIs of the same random unit during all the simulations. On the left, results using the constant = 0.1. Larger panels shows the full trajectory until 103 timesteps after convergence. Smaller panels are a zoom of the first trajectory 0 = 1.1 (top) and = 0.87 (bottom). Right panels show the same type of results but using a smaller constant = 0.01. 3.3 Simulations In this section we show empirical results for the proposed plasticity rule. We focus our analysis on the time conv required for the system to converge toward the critical point. In particular, we analyze how conv depends on the starting initial configuration and on the constant . For the experiments we use a network composed of N = 500 units with homogeneous L = 500 and p = 0.9. Synapses are initialized homogeneously and random initial states are chosen for all units in each trial. Every time a unit i fires, we update its afferent synapses ij , for all j = i, which breaks the homogeneity in the interaction strengths. The network starts with a certain initial condition 0 and evolves according to its original discrete dynamics, Eq. (1), together with plasticity rule (9). To measure the time conv necessary to reach a value close to = 1 for the first time, we select a neuron i randomly and compute every time i fires. We assume convergence when (1 - , 1 + ) for the first time. In these initial experiments, is set to /5 and is either 0.1 or 0.01. We performed 50 random experiments for different initial configurations. In all cases, after an initial transient, the network settles close to = 1, presenting some fluctuations. These fluctuations did not grow even after 106 ISIs in all realizations. Figure 4 shows examples for 0 {0.58, 0.7, 0.87, 1.1, 1.3, 1.7}. We can see that for larger updates of the synapses ( = 0.1) the network converges faster. However, fluctuations around the reached state, slightly above = 1, are approximately one order of magnitude bigger than for = 0.01. We therefore can conclude that determines the speed of convergence and the quality and stability of the dynamics at the critical state: high values of cause fast convergence but turn the dynamics of the network less stable at the critical state. We study now how conv depends on 0 in more detail. Given N , L, c and , we can approximate the global change in after one entire ISI of a random unit assuming that all neurons change its afferent synapses uniformly. This gives us a recursive definition for the sequence of t s generated by the synaptic plasticity rule: -Lef f (t ) - c sgn(t - 1) ( (t ) = (N - 1) , + 2 2 2 Lef f (t ) + 2c) + 2c(L - Lef f (t )) 6 (a) 10 5 Periods required to reach =1 (b) 10 6 Time-steps required to reach =1 time (number of timesteps) 10 4 10 5 time (periods) 10 4 10 3 10 3 10 2 = 0.01 conv 10 2 = 0.01 conv_steps 10 1 Simulations ----------- = 0.1 conv Simulations 10 1 Simulations ----------- = 0.1 conv_steps Simulations 10 0.5 0 1 1.5 2 10 0.5 0 1 0 1.5 2 0 Figure 5: Number of ISIs (a) and timesteps (b) required to reach the critical state in function of the initial configuration 0 . Rounded dots indicate empirical results as averages over 10 different realizations starting from the same 0 . Continuous curves correspond to Eq. (12). Parameter values are N = 500, L = 500, p = 0.9, c = 1, = /5. 1 1 - t a nd where Lef f (t ) = (L - 1) t+1 = t + (t ). Then the number of ISIs and the number of timesteps can be obtained by2 : conv = min({i : |t - 1| }), conv steps = t onv c =0 app (t ). (12) Figure 5 shows empirical values of conv and conv steps for several values of 0 together with the approximations (12). Despite the inhomogeneous coupling strengths, the analytical approximations (continuous lines) of the experiments (circles) are quite accurate. Typically, for 0 < 1 more spikes are required for convergence than for 0 > 1. However, the opposite occurs if we consider timesteps as time units. A hysteresis effect (described in [7]) present in the system if 0 < 1, causes the system to be more resistant against synaptic changes, which increases the number of updates (spikes) necessary to achieve the same effect as for 0 > 1. Nevertheless, since the ISIs are much shorter for supercritical coupling the actual number of time steps is still lower than for subcritical coupling. 4 Discussion Based on the amount of spontaneous evolution which is dissipated during an ISI, we have derived a local synaptic mechanism which causes a network of spiking neurons to self-organize near a critical state. Our motivation differs from those of similar studies, for instance [8], where the average branching ratio of the network is used to characterize criticality. Briefly, is defined as the average number of excitations created in the next time step by a spike of a given neuron. The inverse of plays the role of the branching ratio in our model. If we initialize the units uniformly in [1, L], we have approximately one unit in every subinterval of length , and in consequence, the closest unit to the threshold spikes in 1/ cases if it receives a spike. For > 1, a spike of a neuron rarely induces another neuron to spike, so < 1. Conversely, for < 1, the spike of a single neuron triggers more than one neuron to spike ( > 1). Only for = 1 the spike of a neuron elicits the order of one spike ( = 1). Our study thus represents a realization of a local synaptic mechanism which induces global homeostasis towards an optimal branching factor. This idea is also related to the SOC rule proposed in [3], where a mechanism is defined for threshold gates (binary units) in terms of bit flip probabilities instead of spiking neurons. As in our model, criticality is achieved via synaptic scaling, where each neuron adjusts its synaptic input according to an effective threshold called margin. 2 The value of app (t ) has to be calculated using an corresponding to t in Eq. (3). 7 When the network is operating at the critical regime, the dynamics can be seen as balancing between a predictable pattern of activity and uncorrelated random behavior typically present in SOC. One would also expect to find macroscopic magnitudes distributed according to scale-free distributions. Preliminary results indicate that, if the stochastic evolution is reset to zero (p = 0) at the critical state, inducing an artificial spike on a randomly selected unit causes neuronal avalanches of sizes and lengths which span several orders of magnitude and follow heavy tailed distributions. These results are in concordance with what is usually found for SOC and will be published elsewhere. The spontaneous evolution can be interpreted for instance as activity from other brain areas not considered in the pool of the simulated units, or as stochastic sensory input. Our results indicate that the amount of this stochastic activity that is absorbed by the system is maximized at an optimal state, which in a sense minimizes the possible effect of fluctuations due to noise on the behavior of the system. The application of the synaptic rule for information processing is left for future research. We advance, however, that external perturbations when the network is critical would cause a transient activity. During the transient, synapses could be modified according to some other form of learning to encode the proper values which drive the whole network to attain a characteristic synchronized pattern for the external stimuli presented. We conjecture that the hysteresis effect shown in the regime of < 1 may be suitable for such purposes, since the network then is able to keep the same pattern of activity until the critical state is reached again. Acknowledgments We thank Joaqu´n J. Torres and Max Welling for useful suggestions and interesting discussions. i References [1] P. Bak. How nature works: The Science of Self-Organized Criticality. Springer, 1996. [2] J. M. Beggs and D. Plenz. Neuronal avalanches in neocortical circuits. Journal of Neuroscience, 23(35):11167­11177, December 2003. ¨ [3] N. Bertschinger, T. Natschlager, and R. A. Legenstein. At the edge of chaos: Real-time computations and self-organized criticality in recurrent neural networks. In Advances in Neural Information Processing Systems 17, pages 145­152. MIT Press, Cambridge, MA, 2005. [4] G. Q. Bi and M. M. Poo. Synaptic modifications in cultured hippocampal neurons: Dependence on spike timing, synaptic strength, and postsynaptic cell type. Journal Of Neuroscience, 18:10464­10472, 1998. [5] G. L. Gerstein and B. Mandelbrot. Random walk models for the spike activity of a single neuron. Biophys J., 4:41­68, 1964. ´ ´ [6] V. Gomez, A. Kaltenbrunner, and V. Lopez. Event modeling of message interchange in stochastic neural ensembles. In IJCNN'06, Vancouver, BC, Canada, pages 81­88, 2006. ´ ´ [7] A. Kaltenbrunner, V. Gomez, and V. Lopez. Phase transition and hysteresis in an ensemble of stochastic spiking neurons. Neural Computation, 19(11):3011­3050, 2007. [8] O. Kinouchi and M. Copelli. Optimal dynamical range of excitable networks at criticality. Nature Physics, 2:348, 2006. [9] C. G. Langton. Computation at the edge of chaos: Phase transitions and emergent computation. Physica D Nonlinear Phenomena, 42:12­37, jun 1990. [10] A. Levina, J. M. Herrmann, and T. Geisel. Dynamical synapses causing self-organized criticality in neural networks. Nature Physics, 3(12):857­860, 2007. [11] N. H. Packard. Adaptation toward the edge of chaos. In: Dynamics Patterns in Complex Systems, pages 293­301. World Scientific: Singapore, 1988. A. J. Mandell, J. A. S. Kelso, and M. F. Shlesinger, editors. ´ ´ i [12] F. Rodr´guez, A. Suarez, and V. Lopez. Period focusing induced by network feedback in populations of noisy integrate-and-fire neurons. Neural Computation, 13(11):2495­2516, 2001. [13] S. Song, K. D. Miller, and L. F. Abbott. Competitive hebbian learning through spike-timing-dependent synaptic plasticity. Nature Neuroscience, 3(9):919­926, 2000. [14] G. G. Turrigiano and S. B. Nelson. Homeostatic plasticity in the developing nervous system. Nature Reviews Neuroscience, 5(2):97­107, 2004. [15] C. Van Vreeswijk and L. F. Abbott. Self-sustained firing in populations of integrate-and-fire neurons. SIAM J. Appl. Math., 53(1):253­264, 1993. 8