Influence of inhibitory synapses on the criticality of excitable neuronal networks
F S Borges, P R Protachevicz, V Santos, M S Santos, E C Gabrick, K C Iarosz, E L Lameu, M S Baptista, I L Caldas, A M Batista
aa r X i v : . [ q - b i o . N C ] A ug Influence of inhibitory synapses on the criticality of ex-citable neuronal networks
F S Borges , P R Protachevicz ,V Santos , M S Santos , E C Gabrick , K CIarosz , E L Lameu , M S Baptista , I L Caldas , A M Batista Center for Mathematics, Computation, and Cognition, Federal University of ABC, 09606-045, S˜ao Bernardodo Campo, SP, Brazil. Graduate Program in Science - Physics, State University of Ponta Grossa, 84030-900,Ponta Grossa, PR, Brazil. Institute of Physics, University of S˜ao Paulo, 05508-900, S˜ao Paulo, SP, Brazil. Faculdade de Telˆemaco Borba, 84266-010, Telˆemaco Borba, PR, Brazil. Graduate Program in ChemicalEngineering, Federal University of Technology Paran´a, 84016-210, Ponta Grossa, PR, Brazil. Cell Biol-ogy and Anatomy Department, University of Calgary, AB T2N 4N1, Calgary, AB, Canada. Institute forComplex Systems and Mathematical Biology, University of Aberdeen, AB24 3UE, Aberdeen, Scotland, UK. Department of Mathematics and Statistics, State University of Ponta Grossa, 84030-900, Ponta Grossa, PR,Brazil.Corresponding author: [email protected], [email protected].
Abstract
In this work, we study the dynamic range of a neuronal network of excitableneurons with excitatory and inhibitory synapses. We obtain an analytical expres-sion for the critical point as a function of the excitatory and inhibitory synap-tic intensities. We also determine an analytical expression that gives the criti-cal point value in which the maximal dynamic range occurs. Depending on themean connection degree and coupling weights, the critical points can exhibit ceas-ing or ceaseless dynamics. However, the dynamic range is equal in both cases.We observe that the external stimulus mask some effects of self-sustained activity(ceaseless dynamic) in the region where the dynamic range is calculated. In theseregions, the firing rate is the same for ceaseless dynamics and ceasing activity. Fur-thermore, we verify that excitatory and inhibitory inputs are approximately equalfor a network with a large number of connections, showing excitatory-inhibitorybalance as reported experimentally.
Keywords: inhibitory synapses, dynamic range, self-sustaining dynamics
The relation between stimuli and sensation is one of the main research topics in Psy-chophysics [1]. Stimulus of different sources and intensities can cause different re-sponses in the sensory system [2]. In the early 19th century, Weber and Fechner pro-posed that stimuli-response relation correspond to a logarithmic function [3, 4]. In the1950s, Stevens proposed that stimuli-response relation is given by a power law [5]. Dueto physiological and anatomical limitation, the relation between stimuli and responsehave upper and lower limits. The stimuli difference, between the smaller and biggersensation, define the dynamic range (DR) associated with its sense [6]. In the contextof biological systems, e.g. neuronal networks, the DR corresponds to the ability todifferentiate the intensity of external stimulus [7].The DR is proportional to the logarithm of the ratio between the largest value of theexternal applied stimulus in which the response is close to saturation of the firing rateand the smallest value of the external applied stimulus in which it is weak to modify1he firing rate. The human sense of sight can perceive changes in about ten decadesof luminosity and the hearing covers twelve decades in a range of intensities of soundpressures [5, 8]. The DR of the human vision plays an important role in the design ofdisplay devices [9], where as the hearing case it is relevant to cochlear implants [10].The DR of a neuronal network increases with the network size until it reaches asaturation value [11]. The increase of the DR value is also associated with the increaseof the number of excitatory chemical synapses [12, 13]. Borges et al. [14] reportedthe complementary effect of chemical and electrical synapses on the enhance of theDR. Protachevicz et al. shown that chemical synapses can enhance DR of the neuralnetwork submitted to external stimuli [15].The mammalian brain is composed of excitatory and inhibitory neurons [16]. Thebalance between excitation and inhibition plays a crucial role in the transmission ofinformation, signal propagation, and regular firing of neurons in many brain areas [17,18]. Neuronal networks with excitatory and inhibitory neurons have been consideredto describe the dynamics of primary visual cortex [16, 19], cortical firing patterns [20,21, 22, 23], and synaptic plasticity mechanisms [24, 25, 26, 27].Kinouchi and Copelli [3] proposed a model of an excitable network based on Erd¨os-R´enyi (ER) random graphs [28]. They demonstrated that the DR is maximised at thecritical point of a non equilibrium phase transition. A theoretical approach to study theeffects of network topology on the DR was presented by Larremore et al. [29, 30], inwhich was considered only excitatory nodes. Pei et al. [31] investigated the collectivedynamics of excitatory-inhibitory excitable networks in response to external stimuli.They found that the DR is maximised at the critical point of phase transition whichdepends only on the excitatory connections.The spiking dynamics of a network of excitable excitatory nodes resulting from aninitial stimulus ceases after a typically short time at a critical point [3]. However, wheninhibitory nodes are considered the collective dynamic can become self-sustaining asshown by Larremore et al. [32]. They showed this behaviour considering an additiveprobabilistic model, where excitatory nodes increase the probability of activation oftheir neighbours, and inhibitory nodes decrease the probability. In addition, at criticalpoint the collective dynamics can become self-sustainable (ceaseless dynamic) if afraction of inhibitory nodes is greater than a threshold. However, in their model they didnot consider a refractory period and, for this reason, the neuronal firing rate obtained ishigher than the experimentally observed. When refractoriness is included in the model,it is possible to obtain the critical point leading to realistic firing patterns [33, 34].In this work, we investigate the criticality and dynamic range of a cellular automa-ton modelling a neuronal network in which the neurons are connected by means ofexcitatory and inhibitory chemical synapses [12, 14]. In order to understand the re-lationship between maximisation of the DR and the critical self-sustainable activity,we consider a refractory period in the model like the one proposed by Larremore etal. [32]. With the refractory period, the model exhibits more realistic firing rates andcritical self-sustained activity. In our simulations, we observe a transition from cease-less dynamics to ceasing activity when the mean connection degree of the network isincreased. We observe that the external stimulus mask effects of self-sustained activityin the region where the DR is calculated, and the firing rate is the same for the cease-less dynamics and ceasing activity. Furthermore, we obtain an analytical expression2or the DR as a function of the mean excitatory and inhibitory synaptic intensities. Ina network with a large number of connections, we show that the maximal DR valueoccurs in the critical points where excitatory and inhibitory inputs are approximatelyequal. In this situation, the neuronal network is in a balanced state. Shew et al. [35]showed experimentally that the DR is maximised when the excitatory and inhibitorysynaptic input are balanced. Our work thus provides theoretical explanations for thisexperimental result.The paper is organised as follows. In Section 2, we introduce the model. Section3 presents our analytical results about the dynamic range. In the last Section, we drawour conclusions.
We consider a n states cellular automaton model composed of N excitable elements.The state of each neuron i is described by the variable s i ( i = , ..., . . . n ). In this rep-resentation, each neuronal state is associated with the neuronal activity [3, 36]. Theresting state is given by s i =
0, the excited state by s i =
1, and s i = , ..., n − G ( x i ) [32] G ( x i ) = G N ∑ j = A i j δ ( s j ( t ) , ) ! , (1)where G ( x i ) = x i ≤ G ( x i ) = x i for 0 < x i <
1, and G ( x i ) = x i ≥ G ( x i ) is apiecewise linear function, known as transfer function, with three pieces. The weightedmatrix A has elements A i j > A i j < δ ( a , b ) is equal to 1 when a = b and zero otherwise.The dynamics of both excited and the refractory states are deterministic. If s i =
1, in thenext time steps the state is updated to s i =
2, and so forth, until s i = n −
1, returning tothe resting state s i = f ex and f in , respectively, and the condition f ex + f in = i indexes as 1 ≤ i ≤ f ex N for excitatory nodes, whereas f ex N + ≤ i ≤ N for inhibitory ones. Fig. 1 displaysan schematic illustration of (a) the neuronal dynamics for n = G ( x i ) as a function of the sumof all synaptic inputs.The neuronal response at a given time t can be quantified using the density ofspiking neurons p ( t ) = N N ∑ i = δ ( s i ( t ) , ) , (2)which is interpreted as the probability for a random neuron to be in the excited state at3 G ( ) x i x x i = Σ Nj =1 Aij δ s j i ( ( ),1) t s i i V ( m V ) i time (ms) (c)(b)(a) Figure 1: Representation of the neuronal activity by a cellular automaton with n = i , where s i represent the rest( s i = s i = s i =
2) states. (b) Chemical synaptic inputsarriving in the neuron i . Red triangles and blue squares represent the excitatory andinhibitory inputs, respectively. (c) The neuronal activation probability ( G ( x i ) ) is givenby a function of all chemical inputs arriving in the neuron i at time t .4ime t . With the time series of p ( t ) , we calculate the average firing rate F = p ( t ) = T T ∑ t = p ( t ) , (3)where T is the time window chosen to calculate the average.In this work, we consider random networks and for this case the update equationsare the same for both excitatory and inhibitory nodes [31]. Our networks are built ac-cording to the Erd¨os-R´enyi random graphs with probability equal to K / ( N − ) , where K is the average degree of connections of the network. Assuming that the events of theneighbours of an excited node are statistically independent for large t , we obtain thefollowing mean field map for the density of spiking neurons p ( t + ) = [ − ( n − ) p ( t )]( η + G ( x ) − η G ( x )) , (4)where the external stimulus η = − exp ( − r ∆ t ) is a Poisson process with mean pertur-bation rate r in the time interval ∆ t [3]. In our simulations, we use ∆ t = A i j = S ex for the excitatory connections and A i j = − S in for theinhibitory ones, when the network reaches a stationary state, the mean value of x i isgiven by h x i = f ex KS ex p ( t ) − f in KS in p ( t ) . (5)Defining σ ex = KS ex and σ in = KS in , we obtain h x i = ( f ex σ ex − f in σ in ) p ( t ) . (6)In the stationary state we have p ( t + ) = p ( t ) = p ∗ and F ≈ p ∗ . Substituting in Eq.(4), and considering the case of no external perturbation ( η = F = ( − ( n − ) F ) G ( x ) . (7)In the regime 0 < ( f ex σ ex − f in σ in ) F <
1, the model implies G ( x ) = x , and therefore F = ( − ( n − ) F ) ( f ex σ ex − f in σ in ) F . (8)Solving for F we get F = − ( f ex σ ex − f in σ in ) − n − . (9)There is a phase transition from ceasing activity ( F =
0) to ceaseless activity ( F > F → σ in = f ex σ ex − f in . (10)This relation shows that the critical point in the model is given by f ex σ ex ≥
1, implyingthe necessity of a minimum fraction of excitatory neurons. In addition, we observethat f ex σ ex ≈ f in σ in for σ ex ≫
1. Then, for a highly connected network ( σ ∝ K ), weobtain approximately the same amount of excitatory and inhibitory mean inputs from5 σ in F ( r = ) σ ex = 1.5 σ ex = 2.0 σ ex = 2.5 p ( t ) σ in = 0.5 σ in = 1.0 σ in = 1.5 time (ms) p ( t ) σ in = 4.5 σ in = 5.0 σ in = 5.5 σ ex = 1.5 σ ex = 2.5 (b)(a)(c) Figure 2: Time series of the density of spiking neurons for subcritical (black line),critical (red line), and supercritical (blue line) values of σ in for (a) σ ex = . σ ex = .
5. In (c), we plot the average firing rate as a function of σ in for σ ex = . σ ex = . σ ex = . N = , K = , r = n =
3, and f ex = . f ex , σ ex , f in , and σ in . (i) if F < F = F > σ in in the subcritical, critical, andsupercritical regime for different values of σ ex . We choose randomly 0 .
4% of neuronsto be active ( s i =
1) at t =
0. In Fig. 2(c), we show the relation between F and σ in for some values of σ ex . We verify that the theoretical results given by Eq. 9 are inagreement with our numerical simulations.When a great number of neurons presents x i <
0, even at the critical point, in thenumerical simulations F can be positive for a large time span. In Fig. 2, we seethat the spiking activity ceases rapidly when σ ex = .
5, whereas it is persistent at thecritical point when σ ex = .
5. We verify that the activity is not persistent if we increasethe average degree of connections. Fig. 3(a) exhibits the density of spiking neuronsconsidering K = × , for subcritical, critical (three different initial conditions) andsupercritical values of σ in . In Fig. 3(b), we plot the distribution of x i for 1000 timesteps, N = , K = (blue), and K = × (red). In both cases, we find h x i and F ≈ . .
18% of x i present negative values. Inthe second case, about 5 .
00% of x i are less than zero. We observe that greater values of S ex = σ ex / K and S in = σ in / K contribute for the persistent activity at the critical point. time (ms) -6 -4 -2 p ( t ) σ in = 4.5 σ in = 5.0 σ in = 5.5 -0.005 0 0.005 0.010 0.015 x i -5 -4 -3 -2 -1 fr ac . ( % ) K =2.10 (2.18% x i <0) K =1.10 (5.00% x i <0) σ ex = 2.5 K = 2.10 (a)(b) Figure 3: (a) Time series of the density of spiking neurons for the subcritical (blackline), critical (red line), and supercritical (blue line) values of σ in with σ ex = .
5. (b)Distribution of x i values for the average degree of connections K = × (red) and K = × (blue). Parameters are N = , f ex = . σ ex = .
5, and σ ex = . Dynamic Range (DR)
The behaviour of the average firing rate ( F ) as a function of the external stimulus ( r )shows a minimum and a maximum saturation ( F and F max , respectively) for a range of r values, as shown in Fig. 4. The DR is defined as ∆ =
10 log r high r low , (11)where ∆ is the stimulus interval (measured in dB) in which changes in r can be per-ceived as changes in F , and it is between the disregarding stimuli that cause a responsesmall to be distinguished from F and the saturation F max [3]. The interval [ r low , r high ]is found from its correspondent in F , [ F low , F high ], where F high = F + . ( F max − F ) and F low = F + . ( F max − F ) . -4 -3 -2 -1 r F FFF low r r high F low highmax0 Figure 4: Mean firing rate as a function of intensity stimuli.For 0 < ( f ex σ ex − f in σ in ) F < F = [ − ( n − ) F ][ η + ( f ex σ ex − f min σ in ) F − ( f ex σ ex − f in σ in ) η F ] . (12)Rearranging the terms, we obtain [( n − )( f ex σ ex − f min σ in )( − η )] F + [ + ( n − ) η − ( f ex σ ex − f min σ in )( − η )] F − η = . (13)As η depends on r , by solving Eq. 13, we are able to determine the dependence of theaverage firing rate on the mean perturbation rate r , as well as its dependence on all theparameters of the network.In Fig. 5(a), we plot F as a function of r for subcritical, critical and supercriticalvalues of σ in . The lines represent the theoretical values from the solution of expression13 and the symbols are obtained through numerical simulations. In the inset of Fig.5(a), we show a magnification to demonstrate that there are diferences between thetheoretical and the numerical values of F for r values out of the region where DR iscalculated (green). 8or a cellular automaton with n states, the maximum average firing rate is given by F max = / n . Deriving F in Eq. 9, F low and F high can be obtained. Then, η low and η high can be calculated directly by η low , high = λ F low , high − λ F low , high (cid:20) λ − ( n − ) λ F low , high − (cid:21) , (14)where we substitute λ = f ex σ ex − f in σ in for convenience. Now we calculate r low and r high according to r low , high = − ln | − η low , high | . (15)Using the equations 9, 14, 15, and the expressions for F max , F low , and F high , we calculatethe dynamic range. σ e x σ in
20 25 30 35 ∆ -4 -2 r F σ ex = 1.5 , K = 1·10 σ ex = 2.5 , K = 1·10 σ ex = 2.5 , K = 2·10 -5 -4 -3 -2 -4 -3 -2 -1 F low = 0.05· F max σ in (b)(a)(c) Figure 5: (a) Average firing rate ( F ) as a function of the mean perturbation rate ( r ).The black, red, and blue symbols correspond to subcritical, critical and supercriticalvalues of σ in , respectively. (b) Dynamic range as a function of σ in for σ ex = .
5. Thecoloured circles are obtained by means of simulations and the black lines representthe theoretical results from the analytical expression. Dynamic range of Fig. (b) isindicated on the σ ex × σ in parameter space by a dash line in Fig. (c). We consider N = , K = , and f ex = . σ in =
1, which is the critical point for the considered9arameters ( σ ex = . ( σ ex , σ in ) . The dashed line indicates the range taken in (b). From the figure,we see that the maximum value of DR follows the line given by the critical pointexpression σ in = σ ex − σ ex value increases. Therefore, the model shows boththe critical and balanced state in a network with a great number of connections, wherethe weights are not small. For instance, for K = × and S ex = .
5, we obtain σ ex = KS ex = and the critical σ in = × −
5. The ration between excitatory andinhibitory input is × σ ex σ in ≈ . The firing dynamics of a network of excitable excitatory nodes resulting from an initialstimulus ceases after a typically short time at a critical point. However, when inhibitorynodes are considered the collective dynamic can become self-sustained. In this work,we build a cellular automaton model with excitatory and inhibitory connections. Inour network, we consider that the connections have different weights. We find anexpression that relates the mean of excitatory and inhibitory weights at the criticalpoint. We also calculate an expression for the dynamic range and show that at thecritical point it reaches its maximal value.Depending of mean connection degree and coupling weights, the critical points canexhibit ceasing or ceaseless dynamics (self-sustained activity). However, the dynamicrange is equal in both cases. We observe that the external stimulus mask some effectsof self-sustained activity in the region where the DR is calculated. In these regions,the firing rate is the same for ceaseless dynamics and ceasing activity. Furthermore,we show that at the critical point the amount of excitatory and inhibitory inputs can beapproximately equal in a densely connected network. This result showing excitatory-inhibitory balanced was experimentally reported by Shew et al. [35].In future works, we plan to consider other network topologies, such as small-worldand scale-free, to study the influence of inhibitory synapses on the criticality of ex-citable neuronal networks.
Acknowledgment
This study was possible by partial financial support from the following Brazilian gov-ernment agencies: Fundac¸˜ao Arauc´aria, National Council for Scientific and Techno-logical Development, Coordination for the Improvement of Higher Education Person-nel, and S˜ao Paulo Research Foundation (2015/07311-7, 2017/18977-1, 2018/03211-6,2020/04624-2). 10 eferences [1] G A Chescheider. Psychophysics: The Fundamentals 3th ed., Psychology Press(2013)[2] S Stevens. Psychophysics. Transaction Publisher (1975)[3] O Kinouchi, M Copelli. Nature Physics e1000402 (2009)[8] D R Chialvo. Nat. Phys. , 819-827 (2012)[14] F S Borges, E L Lameu, A M Batista, K C Iarosz, M S Baptista, R L Viana.Physica A , 58-64 (2017)[25] R R Borges, F S Borges, E L Lameu, P R Protachevicz, K C Iarosz, I L Caldas,R L Viana, E E N Macau, M S Baptista, C Grebogi, A. M. Batista. Braz. J. Phys. , 678-688 (2017)[26] E L Lameu, S Yanchuk, E E N Macau, F S Borges, K C Iarosz, I L Caldas, PR Protachevicz, R R Borges, R L Viana, J D Szezech Jr, A M Batista, J Kurths.Chaos , 085701 (2018)[27] E L Lameu, E E N Macau, F S Borges, K C Iarosz, I L Caldas, R R Borges, PR Protachevicz, R L Viana, A M Batista. Eur. Phys. J. Spec. Top. , 673-682(2018)[28] P Erd¨os; A R´enyi. On Random Graphs. Publicationes Mathematicae65