Anticipated Synchronization in a Biologically Plausible Model of Neuronal Motifs
Fernanda S. Matias, Pedro V. Carelli, Claudio R. Mirasso, Mauro Copelli
aa r X i v : . [ q - b i o . N C ] S e p Anticipated Synchronization in a Biologically Plausible Model of Neuronal Motifs
Fernanda S. Matias, Pedro V. Carelli, Claudio R. Mirasso, and Mauro Copelli Departamento de F´ısica, Universidade Federal de Pernambuco, Recife, Pernambuco 50670-901 Brazil Instituto de Fisica Interdisciplinar y Sistemas Complejos, CSIC-UIB,Campus Universitat de les Illes Balears E-07122 Palma de Mallorca, Spain
Two identical autonomous dynamical systems coupled in a master-slave configuration can exhibitanticipated synchronization (AS) if the slave also receives a delayed negative self-feedback. Recently,AS was shown to occur in systems of simplified neuron models, requiring the coupling of the neuronalmembrane potential with its delayed value. However, this coupling has no obvious biological corre-late. Here we propose a canonical neuronal microcircuit with standard chemical synapses, where thedelayed inhibition is provided by an interneuron. In this biologically plausible scenario, a smoothtransition from delayed synchronization (DS) to AS typically occurs when the inhibitory synapticconductance is increased. The phenomenon is shown to be robust when model parameters are variedwithin physiological range. Since the DS-AS transition amounts to an inversion in the timing of thepre- and post-synaptic spikes, our results could have a bearing on spike-timing-dependent-plasticitymodels.
PACS numbers: 87.18.Sn, 87.19.ll, 87.19.lm
I. INTRODUCTION
Synchronization of nonlinear systems has been exten-sively studied on a large variety of physical and biologicalsystems. Synchronization of oscillators goes back to thework by Huygens, and in the past decades an increasedinterest in the topic of synchronization of chaotic systemsappeared [1].About a decade ago, Voss [2] discovered a new schemeof synchronization that he called “anticipated synchro-nization”. He found that two identical dynamical sys-tems coupled in a master-slave configuration can exhibitthis anticipated synchronization if the slave is subjectedto a delayed self-feedback. One of the prototypical exam-ples proposed by Voss [2–4] is described by the equations˙ x = f ( x ( t )) , (1)˙ y = f ( y ( t )) + K [ x ( t ) − y ( t − t d )] .f ( x ) is a function which defines the autonomous dynam-ical system. The solution y ( t ) = x ( t + t d ), which char-acterizes the anticipated synchronization (AS), has beenshown to be stable in a variety of scenarios, includingtheoretical studies of autonomous chaotic systems [2–4],inertial ratchets [5], and delayed-coupled maps [6], as wellas experimental observations in lasers [7, 8] or electroniccircuits [9] .More recently, AS was also shown to occur in a non-autonomous dynamical system, with FitzHugh-Nagumomodels driven by white noise [10–12]. In these works,even when the model neurons were tuned to the excitableregime, the slave neuron was able to anticipate the spikesof the master neuron, working as a predictor [9] . Thoughpotentially interesting for neuroscience, it is not trivial tocompare these theoretical results with real neuronal data.The main difficulty lies in requiring that the membranepotentials of the involved neurons be diffusively coupled.While a master-slave coupling of the membrane poten-tials could in principle be conceived by means of electri- cal synapses (via gap junctions) [13] or ephaptic interac-tions [14], no biophysical mechanism has been proposedto account for the delayed inhibitory self-coupling of theslave membrane potential.In the brain, the vast majority of neurons are cou-pled via chemical synapses, which can be excitatory orinhibitory. In both cases, the coupling is directional andhighly nonlinear, typically requiring a suprathreshold ac-tivation (e.g. a spike) of the pre-synaptic neuron totrigger the release of neurotransmitters. These neuro-transmitters then need to diffuse through the synapticcleft and bind to receptors in the membrane of the post-synaptic neuron. Binding leads to the opening of spe-cific channels, allowing ionic currents to change the post-synaptic membrane potential [13]. This means that notonly the membrane potentials are not directly coupled,but the synapses themselves are dynamical systems.Here we propose to bridge this gap, investigatingwhether AS can occur in biophysically plausible modelneurons coupled via chemical synapses. The model isdescribed in section II. In section III we describe ourresults, showing that AS can indeed occur in “physio-logical regions” of parameter space. Finally, section IVbrings our concluding remarks and briefly discusses thesignificance of our findings for neuroscience, as well asperspectives of future work. II. MODELA. Neuronal motifs
We start by mimicking the original master-slave cir-cuit of eqs. (1) with a unidirectional excitatory chemicalsynapse (M −→ S in Fig. 1a). In a scenario with standardbiophysical models, the inhibitory feedback we proposeis given by an interneuron (I) driven by the slave neu-ron, which projects back an inhibitory chemical synapseto the slave neuron (see Fig. 1a). So the time-delayednegative feedback is accounted for by chemical inhibitionwhich impinges on the slave neuron some time after ithas spiked, simply because synapses have characteristictime scales. Such inhibitory feedback loop is one of themost canonical neuronal microcircuits found to play sev-eral important roles, for instance, in the spinal cord [15],cortex [15], thalamus [16, 17] and nuclei involved withsong production in the bird brain [18]. For simplicity, wewill henceforth refer to the 3-neuron motif of Fig. 1a asa Master-Slave-Interneuron (MSI) system. (a)
M S I g A g A g G (b) M S I g A g A g G g N D g N g N FIG. 1. (a) (Color online) Three neurons coupled by chemicalsynapses in the master-slave-interneuron (MSI) configuration: excitatory AMPA synapses (with maximal conductance g A )couple master (M) to slave (S) and slave to interneuron (I),whereas an inhibitory GABA A synapse (with maximal con-ductance g G ) couples interneuron to slave. (b) Same as (a),except that all three neurons of the MSI circuit receive exci-tatory (NMDA) synapses from a driver neuron (D). As we show in section III below, whether or not theMSI circuit can exhibit AS depends, among other fac-tors, on the excitability of the three neurons. In theMSI, this is controlled by a constant applied current (seesection III A). To test the robustness of the results (andat the same time improve the realism and complexity ofthe model), in section III B we study the four-neuron mo-tif depicted in Fig. 1b, where the excitability of the MSInetwork is chemically modulated via synapses project-ing from a global driver ( D ). From now on, we refer tothe 4-neuron motif as a Driver-Master-Slave-Interneuron(DMSI) microcircuit. B. Model neurons
In the above networks, each node is described by aHodgkin-Huxley (HH) model neuron [19], consisting offour coupled ordinary differential equations associated to the membrane potential V and the ionic currents flowingacross the axonal membrane corresponding to the Na, Kand leakage currents. The gating variables for sodiumare h and m and for the potassium is n . The equationsread [20]: C m dVdt = G Na m h ( E Na − V ) + G K n ( E K − V )+ G m ( V rest − V ) + I + X I syn (2) dxdt = α x ( V )(1 − x ) − β x ( V ) x , (3)where x ∈ { h, m, n } , C m = 9 π µ F is the membrane ca-pacitance of a 30 × × π µ m equipotential patch ofmembrane [20], I is a constant current which sets the neu-ron excitability and P I syn accounts for the interactionwith other neurons. The reversal potentials are E Na =115 mV, E K = −
12 mV and V rest = 10 . G Na = 1080 π mS, G K = 324 π mS and G m = 2 . π mS, respectively. Thevoltage dependent rate constants in the Hodgkin-Huxleymodel have the form: α n ( V ) = 10 − V e (10 − V ) / − , (4) β n ( V ) = 0 . e − V/ , (5) α m ( V ) = 25 − V e (25 − V ) / − , (6) β m ( V ) = 4 e − V/ , (7) α h ( V ) = 0 . e − V/ , (8) β h ( V ) = 1( e (30 − V ) / + 1) . (9)Note that all voltages are expressed relative to the restingpotential of the model at I = 0 [20].According to Rinzel and Miller [21], in the absenceof synaptic currents the only attractor of the system ofequations (2)-(9) for I . .
13 pA is a stable fixed point,which loses stability via a subcritical Hopf bifurcation at I ≃ .
51 pA. For 177 .
13 pA . I . .
51 pA, thestable fix point coexists with a stable limit cycle.
C. Synaptic coupling
AMPA (A) and GABA A (G) are the fast excitatoryand inhibitory synapses in our model [see Fig. 1a]. Fol-lowing Destexhe et al [22], the fraction r ( i ) ( i = A, G )of bound (i.e. open) synaptic receptors is modelled by afirst-order kinetic dynamics: dr ( i ) dt = α i [ T ](1 − r ( i ) ) − β i r ( i ) , (10)where α i and β i are rate constants and [ T ] is the neuro-transmitter concentration in the synaptic cleft. For sim-plicity, we assume [ T ] to be an instantaneous function ofthe pre-synaptic potential V pre :[ T ]( V pre ) = T max e [ − ( V pre − V p ) /K p ] , (11)where T max = 1 mM − is the maximal value of [ T ], K p =5 mV gives the steepness of the sigmoid and V p = 62 mVsets the value at which the function is half-activated [22].The synaptic current at a given synapse is given by I ( i ) = g i r ( i ) ( V − E i ) , (12)where V is the postsynaptic potential, g i the maximalconductance and E i the reversal potential. We use E A =60 mV and E G = −
20 mV.The values of the rate constants α A , β A , α G , and β G are known to depend on a number of different factorsand vary significantly [23–26]. To exemplify some of ourresults, we initially fix some parameters, which are setto the values of Table 1 unless otherwise stated (sec-tion III A 1). Then we allow these parameters (as wellas the synaptic conductances) to vary within the phys-iological range when exploring different synchronizationregimes (see section III A 2 and III B).The slow excitatory synapse is NMDA (N) and itssynaptic current is given by: I ( N ) = g N B ( V ) r ( N ) ( V − E N ) , (13)where E N = 60 mV. The dynamics of the variable r ( N ) is similar to eq. (10) with α N = 0 .
072 mM − ms − and β N = 0 . − . The magnesium block of the NMDAreceptor channel can be modeled as a function of postsy-naptic voltage V : B ( V ) = 11 + e ( − . V )[Mg ] o / . , (14)where [Mg ] o = 1 mM is the physiological extracellularmagnesium concentration.In what follows, we will drop the neurotransmitter su-perscripts A , G and N from the synaptic variables r and I . Instead we use double subscripts to denote the re-ferred pre- and postsynaptic neurons. For instance, thesynaptic current in the slave neuron due to the interneu-ron (the only inhibitory synapse in our models) will bedenoted as I IS , and so forth. MSI DMSIα A (mM − ms − ) 1 . . β A (ms − ) 0 .
19 0 . α G (mM − ms − ) 5 . . β G (ms − ) 0 .
30 0 . α N (mM − ms − ) — 0 . β N (ms − ) — 0 . g A (nS) 10 10 I (pA) 280 160TABLE I. Standard values employed in the model. See textfor details. III. RESULTSA. Master-Slave-Interneuron circuits
1. Three dynamical regimes
Initially, we describe results for the scenario where allneurons receive a constant current I >
280 pA. Thiscorresponds to a situation in which the fixed points areunstable and, when isolated, they spike periodically. Allother parameters are as in Table 1. For different sets ofinhibitory conductance values g G our system can exhibitthree different behaviors. To characterize them, we define t Mi as the time the membrane potential of the masterneuron is at its maximal value in the i -th cycle (i.e. its i -th spike time), and t Si as the spike time of the slaveneuron which is nearest to t Mi .The delay τ i is defined as the difference (see Fig. 2): τ i ≡ t Mi − t Si . (15)Initial conditions were randomly chosen for each com-puted time series. When τ i converges to a constant value τ , a phase-locked regime is reached [27]. If τ < τ > τ does not convergeto a fixed value, the system is in a phase drift (PD)regime [27]. The extent to which the AS regime can be le-gitimately considered “anticipated” in a periodic systemwill be discussed below.In Figure 3 we show examples of time series in the threedifferent regimes (DS, AS and PD). The different panelscorrespond to the membrane potential, fraction of acti-vated receptors for each synapse, and synaptic current inthe slave neuron. For a relatively small value of the in-hibitory coupling [ g G = 20 nS, Fig. 3(a)] the slave neuronlags behind the master, characterizing DS. In Fig. 3(b),we observe that by increasing the value of the inhibitorycoupling ( g G = 40 nS) we reach an AS regime. Finally,for strong enough inhibition [ g G = 60 nS, Fig. 3(c)] thePD regime ensues.In the DS and AS regimes the master and slave neu-rons spike at the same frequency. However, when thesystem reaches the PD regime the mean firing rate of theslave neuron becomes higher than that of the master.The counterintuitive result shown in Fig. 4(a) emerges:the mean firing rate of the slave neuron increases whileincreasing the conductance of the inhibitory synapse pro-jected from the interneuron. For the particular combina-tion of parameters used in Fig. 4(a), the transition turnsout to be reentrant, i.e., the system returns to the DSregime for sufficiently strong inhibition (a more detailedexploration of parameter space will be presented below).Figure 4(b) shows the return map of the inter-spike in-terval of the slave, which forms a closed curve (touchingthe trivial single-point return map of the master). Thisis consistent with a quasi-periodic phase-drift regime. (a) t (ms) V ( m V ) MSI t i - t i+1 τ = t i - t i (DS) S SM M (b) t (ms) V ( m V ) MSI t i - t i -1 τ = t i - t i (AS) SS M M FIG. 2. (Color online) Membrane potential V as a functionof time for an external current I = 280 pA in the master (M),slave (S), and interneuron (I) neurons. The plot illustratestwo regimes: (a) g G = 20 nS leads to delayed synchronization(DS), where τ <
0, and (b) g G = 40 nS leads to anticipatedsynchronization (AS), where τ >
0. Other parameters as inTable 1.
Note that in this simple scenario g G plays an analo-gous role to that of K in Eq. 1, for which AS is stableonly when K > K c (eventually with reentrances) [28].Moreover, the behavior of the synaptic current in theslave neuron is particularly revealing: in the DS regime[Fig. 3(a)], it has a positive peak prior to the slave spike,which drives the firing in the slave neuron. In the ASregime [Fig. 3(b)], however, there is no significant re-sulting current, except when the slave neuron is alreadysuprathreshold. In this case, the current has essentiallyno effect upon the slave dynamics. This situation is sim-ilar to the stable anticipated solution of Eq. 1, when thecoupling term vanishes. (a) DS V ( m V ) MSI00.51 r r MS r SI r IS t (s) -10 I s y n ( n A ) I MS I IS I MS + I IS (b) AS V ( m V ) r t (s) -10 I s y n ( n A ) (c) PD V ( m V ) r t (s) -10 I s y n ( n A ) FIG. 3. (Color online) Time series of the membrane poten-tials ( V ), bound receptors ( r ) and synaptic currents ( I ), withmodel parameters as in Table 1. Note that the system is peri-odic in the DS and AS regimes [(a) and (b) respectively], butnot in the PD regime (c).
2. Scanning parameter space
The dependence of the time delay τ on g G is shown inFig. 5 for different values of the external current I andmaximal excitatory conductance g A . Several features inthose curves are worth emphasizing. First, unlike previ-ous studies on AS, where the anticipation time was hard- g G (nS) F S / F M
11 12 13 14 15 t Si-1 - t
Si-2 (ms) t S i - t S i - ( m s ) g G = 60 nS g G = 80 nS g G = 70 nS g G =100 nS (a) (b) FIG. 4. (Color online) (a) The mean firing rate of the slave( F S ) coincides with the mean firing rate of the master ( F M )for DS and AS regimes, but it is larger for PD. (b) In PD, thereturn map of the inter-spike interval of the slave is consistentwith a quasi-periodic system (the pink star shows the returnmap of the master). wired via the delay parameter t d [see eq.(1)], in our casethe anticipation time τ is a result of the dynamics. Notethat g G (the parameter varied in Fig. 5) does not changethe time scales of the synaptic dynamical variables ( r ),only the synaptic strength.Secondly, τ varies smoothly with g G . This continu-ity somehow allows us to interpret τ > g G = 0, we simply have a master-slave configurationin which the two neurons spike periodically. Due to theexcitatory coupling, the slave’s spike is always closer tothe master’s spike which preceded it than to the mas-ter’s spike which succeeded it [as in e.g. Fig. (2)(a)].Moreover, the time difference is approximately 1 . g G , the time difference between themaster’s and the slave’s spikes eventually changes sign[as in e.g. Fig. (2)(b)]. Even though the ambiguity inprinciple remains, there is no reason why we should notcall this regime “anticipated synchronization” (again aphase-locked regime, but with a phase difference of op-posite sign). In fact, we have not found any parameterchange which would take the model from the situationin Fig. (2)(a) to that of Fig. (2)(b) by gradually increas-ing the lag of the slave spike until it approached the nextmaster spike. If that ever happened, τ would changediscontinuously (by its definition). Therefore, the term“anticipated synchronization” by no means implies vio-lation of causality and should just be interpreted withcaution. As we will discuss in section IV, the relative timing between pre- and postsynaptic neurons turns outto be extremely relevant for real neurons.Third, it is interesting to note that the largest antic-ipation time can be longer (up to 3 ms, i.e. about 20%of the interspike interval) than the largest time for thedelayed synchronization ( ≈ . g G further in an attempt to obtain even larger values of τ ,however, the system undergoes a bifurcation to a regimewith phase drift [which marks the end of the curves inFig. 5]. g G (nS) -2-101234 τ ( m s ) I = 280 pA I = 300 pA I = 320 pA g A = 1 nS g A = 5 nS g A = 10 nS g G = 20 nS g G = 40 nS g G = 60 nS(DS) (AS) (PD) FIG. 5. (Color online) Dependence of the time delay τ withthe maximal conductance g G for different values of the appliedcurrent I and g A . The end of each curve (stars) marks thecritical value of g G , above which the system changes from ASto PD. The number of parameters in our model is very large.The number of dynamical regimes which a system ofcoupled nonlinear oscillators can present is also verylarge, most notably p/q -subharmonic locking structured g G ( n S ) g A (nS)
10 20 30 40 50 0 20 40 60 80 100 -2-1 0 1 2 3 4 (PD) (DS)(AS)
FIG. 6. (Color online) Delay τ (right bar) in the ( g A , g G ) pro-jection of parameter space: DS (blue, right), AS (red, middle)and PD (white, left — meaning that no stationary value of τ was found). β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS) β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS) β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS) -2-1 0 1 2 3 4 β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS)
10 20 30 40 50020406080100 β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS)
10 20 30 40 50 β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS)
10 20 30 40 50 -2-1 0 1 2 3 4
FIG. 7. (Color online) Delay τ (right bar) in the ( g A , g G )projection of parameter space for different combinations of β A and β G . From left to right we have respectively PD, ASand DS regimes, as in Fig. 6. in Arnold tongues [29]. These occur in our model, but notin the parameter region we are considering. In this con-text, an attempt to map all the dynamical possibilities inparameter space would be extremely difficult and, mostimportant, improductive for our purposes. We there-fore focus on addressing the main question of this paper,which is whether or not AS can be stable in a biophysi-cally plausible model.In Fig. 6 we display a two-dimensional projection ofthe phase diagram of our model. We employ the val-ues in Table I, except for g A , which is varied along thehorizontal axis. Note that each black curve with circlesin Fig. 5 corresponds to a different vertical cut of Fig. 6,along which g G changes. We observe that the three differ-ent regimes are distributed in large continuous regions,having a clear transition between them. Moreover thetransition from the DS to the AS phase can be well ap-proximated by a linear relation g G /g A ≈ . β G and β A . We observe that AS remainsstable in a finite region of the parameter space, and thisregion increases as excitatory synapses become faster.Figure 5 suggests that larger values of the input cur-rent I eventually lead to a transition from AS from DS.This effect is better depicted in Fig. 8, where the DS re-gion increases in size as I (and therefore the firing rate)increases. Figures 8(b)-(d) also show that the system canexhibit reentrant transitions as g G is varied. Most impor-tantly, however, is that Figs. 7 and 8 show that there isalways an AS region in parameter space, as synaptic andintrinsic parameters are varied.As we will discuss in section IV, the possibility of con-trolling the transition between AS and DS is in prin-ciple extremely appealing to the study of plasticity in I = 280 pA I = 300 pA I = 320 pA g G ( n S ) g G ( n S ) g A (nS) g A (nS)(a) (b)(c) (d) I = 340 pA I = 280 pA I = 300 pA I = 320 pA g G ( n S ) g G ( n S ) g A (nS) g A (nS)(a) (b)(c) (d) -2-1 0 1 2 3 4 I = 340 pA I = 280 pA I = 300 pA I = 320 pA g G ( n S ) g G ( n S ) g A (nS) g A (nS)(a) (b)(c) (d)
10 20 30 40 50020406080100 I = 340 pA I = 280 pA I = 300 pA I = 320 pA g G ( n S ) g G ( n S ) g A (nS) g A (nS)(a) (b)(c) (d)
10 20 30 40 50 -2-1 0 1 2 3 4 I = 340 pA FIG. 8. (Color online) Delay τ (right bar) in the ( g A , g G )projection of parameter space for different values of I . PD,AS and DS regimes as in Fig. 6. neuroscience. However, in a biological network, the in-put current would not be exactly constant, but rather bemodulated by other neurons. In the following, we test therobustness of AS in this more involved scenario, thereforemoving one step ahead in biological plausibility. B. Driver-Master-Slave-Interneuron circuits
Let us consider the MSI circuit under a constant inputcurrent I = 160 pA. This is below the Hopf bifurca-tion [21], i.e. none of the three neurons spikes tonically.Their activity will now be controlled by the driver neu-ron (D), which projects excitatory synapses onto the MSIcircuit [see Fig. 1(b)]. We chose to replace the constantinput current by a slowly varying current, so that thesynapses projecting from the driver neuron are of theNMDA type (see section II). The driver neuron receivesa current I D = 280 pA, so it spikes tonically. All remain-ing parameters are as in the second column of Table 1.The interest in this case is to verify whether AS holdswhen the excitability of the MSI circuit is modulated bya non-stationary current.As shown in Fig. 9, we found in this new scenario asimilar route from DS to AS, and then the PD regime(compare with Fig. 7). Note that the characteristic time( β N = 6 . − ) for the unbinding of the NMDA receptorsis about ten times larger than the inter-spike interval ofthe driver neuron (which spikes at ≈
67 Hz). As a conse-quence, r DM , r DS , r DI are kept at nearly constant values(with variations of ≈
10% around a mean value — datanot shown). The variations in the NMDA synaptic cur-rent are also small, which in principle should make thesystem behave in an apparently similar way to the pre-vious MSI circuit. However, these small variations areimportant enough to increase the AS domain in parame-ter space, in some cases even eliminating the PD region β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS) β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS) β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS) -2-1 0 1 2 3 4 5 6 7 β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS)
10 20 30 40 50020406080100 β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS)
10 20 30 40 50 β G = . m s - β G = . m s - β A = 0.19 ms -1 β A = 0.30 ms -1 β A = 0.60 ms -1 g G ( n S ) g G ( n S ) g A (nS) g A (nS) g A (nS)
10 20 30 40 50 -2-1 0 1 2 3 4 5 6 7
FIG. 9. (Color online) DMSI circuit (see Fig. 1b). Delay τ (right bar) in the ( g A , g G ) projection of parameter space fordifferent combinations of β A and β G . PD, AS and DS regimesas in Fig. 6. g G (nS) F S / F M
10 11 12 13 14 15 t Si-1 - t
Si-2 (ms) t S i - t S i - ( m s ) g G = 80 nS g G = 120 nS g G = 200 nS g G = 300 nS (b)(a) FIG. 10. (Color online) DMSI circuit (see Fig. 1b). (a) Themean firing rate of the slave ( F S ) coincides with the meanfiring rate of the master ( F M ) for DS and AS regimes, but itis larger for PD. (b) In PD, the return map of the inter-spikeinterval of the slave is consistent with a quasi-periodic system. (see e.g. Fig. 9 for β G = 0 .
30 ms − ). Therefore, at leastin this case, the use of more biological plausible parame-ters does not destroy AS, but rather enhances it.In fact, the three regions in the MSI diagrams seemto retain their main features in the DMSI circuit. WhenPD occurs, for example, the slave again spikes faster thanthe master [Fig. 10(a)], like in the MSI circuit [comparewith Fig. 4(a)]. Another signature of the robustness ofthe PD phase against the replacement of a constant bya slowly-varying synaptic current appears in the returnmap shown in Fig. 10(b). It can be seen that it has thesame structure of its three-neuron counterpart shown in Fig. 4(b). IV. CONCLUDING REMARKS
In summary, we have shown that a biologically plausi-ble model of a 3-neuron (MSI) motif can exhibit an at-tractor in phase space where anticipated synchronizationis stable. The transition from the DS to the AS regime isa smooth function of the synaptic conductances. Typi-cally, a further increase in the inhibitory conductance g G leads to a second transition from AS to PD, a quasiperi-odic regime in which the slave firing frequency is largerthan that of the master.We have varied synaptic decay rates ( β ), synaptic con-ductances ( g ) as well as input currents ( I ) within well ac-cepted physiological ranges [23–26]. In all the scenariosthere is always a continuous region in parameter spacewhere AS is stable. Replacing the constant current bya global periodic driver (arguably a more realistic sit-uation), we obtain a model of a 4-neuron (DMSI) mo-tif which exhibits the same three regions of the simplermodel. The synaptic rise constants ( α ) were also var-ied, but have a lesser effect on the transitions among thedifferent regimes (data not shown). Therefore the phe-nomenon seems to be robust at the microcircuit scale.It is important to emphasize that our AS results differfrom those obtained from eq. (1) at a fundamental level.In our model the delayed feedback that leads to AS isgiven by biologically plausible elements (an interneuronand chemical synapses). Hence, the anticipation time isnot hard-wired in the dynamical equations, but ratheremerges from the circuit dynamics. Moreover, the par-ticular circuit we study is a neuronal motif ubiquitouslyfound in the brain [15–17]. We are unaware of other ASmodels in which every parameter has a clear biologicalinterpretation.We believe that our results can be extremely rele-vant for modeling studies of synaptic plasticity. Recentdecades have witnessed a growing literature on spike-timing dependent plasticity (STDP), which accountsfor the enhancement or diminution of synaptic weight[long term potentiation (LTP) and long term depression(LTD), respectively] depending on the relative timing be-tween the spikes of the pre- and post-synaptic neurons(see e.g. [30–32]). Experimental data strongly suggestthat if the pre-synaptic neuron fires before (after) thepost-synaptic neuron, the synapse between them will bestrenghtened (weakened) [33, 34]. STDP is supposedto take place in a window of time differences betweenpost- and pre-synaptic spikes in the order of ten millisec-onds, within which the delay and anticipation times ofour models fall. Since the DS-AS transition amounts toan inversion in the timing of the pre- and post-synapticspikes, then by appropriately controlling this effect onecould dynamically toggle between synaptic strengtheningand weakening. This could be potentially linked withmodeling of large-scale ascending feedback modulationfrom reward systems.Our results, therefore, offer a number of possibilitiesfor further investigation. Including effects from microcir-cuit dynamics (such as the ones we have presented here)in models of synaptic plasticity is a natural next step, onewhich we are currently pursuing. Once we have verifiedAS in a biologically plausible model, one could considerusing simplified models [35, 36] (e.g. by replacing the HHequations and/or the synaptic kinetics) and the influenceof noise [10, 20]. We are also investigating whether thestructure of the phase diagram can be qualitatively repro-duced via a phase-response-curve analysis [37, 38] of theneuronal motifs studied here. Results will be published elsewhere. ACKNOWLEDGMENTS
We thank CNPq, FACEPE, CAPES and special pro-grams PRONEX, INCeMaq and PRONEM for financialsupport. MC is grateful to the hospitality of the IFISC-UIB group at Palma de Mallorca, where these ideas werefirst developed. CM acknowledges support from the Min-isterio de Educaci´on y Ciencia (Spain) and Fondo Eu-ropeo de Desarrollo Regional (FEDER) under projectFIS2007-60327 (FISICOS). [1] S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares, andC. S. Zhou, Phys. Rep. , 1 (2002).[2] H. U. Voss, Phys. Rev. E , 5115 (2000).[3] H. U. Voss, Phys. Rev. E , 039904 (2001).[4] H. U. Voss, Phys. Rev. Lett. , 014102 (2001).[5] M. Kostur, P. H¨anggi, P. Talkner, and J. L. Mateos,Phys. Rev. E , 036210 (2005).[6] C. Masoller and D. H. Zanette, Physica A , 359(2001).[7] S. Sivaprakasam, E. M. Shahverdiev, P. S. Spencer, andK. A. Shore, Phys. Rev. Lett. , 154101 (2001).[8] Y. Liu, Y. Takiguchi, P. Davis, T. Aida, and S. Saito,Appl. Phys. Lett. , 4306 (2002).[9] M. Ciszak, C. R. Mirasso, R. Toral, and O. Calvo,Phys. Rev. E , 046203 (2009).[10] M. Ciszak, O. Calvo, C. Masoller, C. R. Mirasso, andR. Toral, Phys. Rev. Lett. , 204102 (May 2003).[11] M. Ciszak, F. Marino, R. Toral, and S. Balle,Phys. Rev. Lett. , 114102 (2004).[12] R. Toral, C. Masoller, C. R. Mirasso, M. Ciszak, andO. Calvo, Physica A , 192 (2003).[13] Essentials of Neural Science and Behavior , edited byE. R. Kandel, J. H. Schwartz, and T. M. Jessell (Ap-pleton & Lange, Norwalk, 1995).[14] A. Arvanitaki, J. Neurophysiol , 89 (1942).[15] The Synaptic Organization of the Brain , edited by G. M.Shepherd (Oxford University Press, New York, 1998).[16] U. Kim, M. V. Sanchez-Vives, and D. A. McCormick,Science , 130 (1997).[17] D. Debay, J. Wolfart, Y. Le Franc, G. Le Masson, andT. Bal, J. Physiol., 540(2004).[18] G. B. Mindlin and R. Laje,
The Physics of Birdsong (Springer, 2005).[19] A. L. Hodgkin and A. F. Huxley, J. Neurophysiol. ,500 (1952).[20] C. Koch,