title: Phase-bistability between anticipated and delayed synchronization in neuronal populations
authors: Julio Nunes Machado Fernanda Selingardi Matias
text: 
  
  
  
  
  
  
  
  
  
  
  Instituto de Física, Universidade Federal de Alagoas, Maceió, Alagoas 57072-970 Brazil.
  
  
  
  fernanda@fis.ufal.br
  Instituto de Física, Universidade Federal de Alagoas, Maceió, Alagoas 57072-970 Brazil.
  
  
  
  
  
  
  Two dynamical systems unidirectionally coupled in a sender-receiver configuration can synchronize with a nonzero phase-lag. In particular, the system
  can exhibit anticipated synchronization (AS), which is characterized by a negative phase-lag, if the receiver (R) also receives a delayed negative self-feedback. 
  Recently, AS was shown to occur between cortical-like neuronal populations in which the self-feedback is mediated by inhibitory synapses.
  In this biologically plausible scenario, a transition from the usual delayed synchronization (DS, with positive phase-lag) to AS
  can be mediated by the inhibitory conductances in the receiver population.
  Here we show that depending on the relation between excitatory and inhibitory synaptic conductances the system can also exhibit 
  phase-bistability between anticipated and delayed synchronization.
  Furthermore, we show that the amount of noise at the receiver and the synaptic conductances 
  can mediate the transition from stable phase-locking to a bistable regime and eventually to a phase-drift (PD).
  We suggest that our spiking neuronal populations model could be potentially useful to study phase-bistability in cortical regions related to bistable perception.
  
  
  
  Phase-bistability between anticipated and delayed synchronization in neuronal populations
      Fernanda Selingardi Matias
        ;  
  =========================================================================================
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Multistable perception is the brain's ability to alternate between two or more perceptual states that occur when sensory information is ambiguous <cit.>.
  It has been hypothesized that the perception of two different features of the same visual cue such as in the Necker cube <cit.>
  (see Fig. <ref>(a), which can be interpreted to have either the upper-left or the lower-right square as its front side)
  should be related to the nonlinear interaction between brain rhythms <cit.>. 
  This means that an almost fixed anatomical connectivity should allow flexible changes from one functional connectivity 
  pattern to another on timescales relevant to behavior. 
  Recently it has been reported that bistable phase-differences in magnetoencephalography (MEG) recordings appear when participants listening to bistable speech sequences that could be perceived as two distinct word sequences repeated over time <cit.>.
  This result suggests that phase-bistability in cortical regions could be related to bistable perception.
  
  
  Here we propose a model of two neuronal populations that can synchronize with a bistable phase-lag.
  Typically, when two dynamical systems are unidirectionally coupled they can synchronize with a positive phase-lag in which the sender is also the leader. 
  This regime is called delayed synchronization (DS).
  However, it has been shown that in a sender-receiver configuration if the receiver is subjected to a delayed self-feedback the system can synchronize with a negative phase-lag <cit.>. 
  This counterintuitive situation indicates that the receiver leads the sender and it is called anticipated synchronization (AS). 
  Formally, Voss has shown that if a system is described by the following equations:
  
      Ṡ    =     f( S(t)), 
      Ṙ    =     f( R(t)) +  K[ S(t)- R(t-t_d)],
  
  with arbitrary continuous f and coupling matrix K, it may has a stable solution R(t)= S(t-t_d), which characterizes AS i.e. the receiver predicts the sender. 
  
  In the last two decades AS has been verified in both theoretical <cit.> and
  experimental <cit.> studies. AS can also occur if the delayed self-feedback is replaced by different parameter mismatches at the receiver <cit.>, 
  including inhibitory loops mediated by chemical synapses <cit.>.
  Especially it has been shown that a faster internal dynamics of the receiver could promote AS between unidirectionally coupled oscillators <cit.>.
  
  Recently, AS has been verified between unidirectionally coupled cortical-like populations <cit.>. 
  The neuronal population model exhibits a smooth transition from DS to AS that can be mediated by synaptic conductances as well as by external noise.
  The model could explain electrophysiological results in non-human primates showing that
  unidirectional Granger causality can be accompanied by both positive or negative phase difference between cortical areas <cit.>.
  Furthermore, AS in the alpha band has been recently reported between synchronized electrodes in human EEG <cit.> as well as in human behavior <cit.>.
  
  
  Here we show that a simple biologically plausible motif with two unidirectionally coupled neuronal populations 
  can exhibit phase-bistability between delayed and anticipated synchronization regimes.
  Especially our populations present a bimodal distribution of phase-differences with a positive (DS) and a negative (AS) peak. 
  In Sec. <ref> we describe our motif as well as neuronal and synaptic models.
  In Sec. <ref>, we report our results, showing that the phase-bistability can be mediated by synaptic couplings and noise.
  Concluding remarks and a brief discussion of the significance of our findings for neuroscience are presented in Sec. <ref>.
  
  
  
  
  § NEURONAL POPULATIONS MODEL
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Our neuronal motif is composed of two unidirectionally coupled cortical-like 
  neuronal populations: a sender (S) and a receiver (R),
  see Fig. <ref>(b). Each one is composed of 400 excitatory and 100 inhibitory
  neurons <cit.> described by the Izhikevich
  model <cit.>:
  
      dv/dt   =    0.04v^2+5v+140-u +∑_x I_x, 
      du/dt   =    a(bv-u).
  
  In Eqs. <ref> and <ref> v is the membrane potential and
  u the recovery variable which accounts for activation
  of K^+ and inactivation of Na^+ ionic currents.  I_x are the synaptic currents provided by
  the interaction with other neurons and external inputs.  If
  v≥30 mV, v is reset to c and u to u+d.  To account for
  the natural heterogeneity of neuronal populations, which can exhibit a variety of neuronal dynamics (regular spiking, 
  bursting, chatering, fast spiking, etc. <cit.>), 
  the dimensionless parameters are randomly sampled as follows: (a,b)=(0.02,0.2) and
  (c,d)=(-65,8)+(15,-6)σ^2 for excitatory neurons and (a,b)=(0.02,0.25)+(0.08,-0.05)σ and
  (c,d)=(-65,2) for inhibitory neurons, where σ is a
  random variable uniformly distributed on the interval [0,1] <cit.>. 
  Equations were integrated with the Euler method and a time step of 0.05 ms.
  
  
  The connections between neurons in each population are assumed to be
  fast unidirectional excitatory and inhibitory chemical synapses
  mediated by AMPA and GABA_A. The synaptic currents are given
  by:
  
      I_x = g_xr_x(v-V_x),
  
  where x=E,I (excitatory and inhibitory mediated by AMPA and
  GABA_A, respectively), V_E=0 mV, V_I=-65 mV, g_x is the maximal
  synaptic conductance and r_x
  is the fraction of bound synaptic receptors whose dynamics is given
  by:
  
      τ_xdr_x/dt=-r_x + D ∑_k δ(t-t_k),
  
  where the summation over k stands for pre-synaptic spikes at times
  t_k. D is
  taken, without loss of generality, equal to 0.05. The time decays are τ_E=5.26 ms τ_I=5.6 ms.  Each
  neuron is subject to an independent noisy spike train described by a
  Poisson distribution with rate R. The input mimics excitatory synapses 
  (with conductances g_P) from n pre-synaptic neurons external to the population,
  each one spiking with a Poisson rate R/n which, together with a constant external current I_c,  
  determine the main frequency of mean membrane potential of each population. Unless otherwise stated, we have employed R=2400 Hz and I_c=0. 
  We have fixed the Poissonian synaptic conductance at the sender population g_P^S=0.5 nS and varied g_P only at the receiver population.  
  
  Connectivity within each
  population randomly targets 10% of the neurons, with excitatory
  conductances set at g^S_E =g^R_E = 0.5 nS. Inhibitory conductances are fixed at the sender population g^S_I = 4.0 nS 
  and g_I at the receiver population is varied throughout the study (see Fig. <ref>).
  Each neuron at the R population receives 
  20 fast synapses (with conductance g_E) from random excitatory
  neurons of the S population.
  The bistability studied in this paper only happens when the synaptic conductances g_E, g_I, and g_P have comparable values.
  The phase-locking regimes presented by the model for g_E≥ 0.5 nS and g_I≥ 1.0 nS have been studied in Ref. <cit.>.
  
  
  
  
  
   §.§ Characterizing time delay in the model
  
  The mean membrane potential V_x (x=S, R) is calculated as the total sum of the membrane potential 
  v of each neuron in the x population in a given time t, divided by the total number of neurons in that population.
  Fig. <ref>(a) and (b) show examples of V_S and  V_R in different regimes. 
  Since V_x, which we assume as a crude approximation of the measured local field potential (LFP), is
  noisy, we average within a sliding window of width 5-8 ms to obtain a smoothened signal, from which we can extract 
  the peak times {t^x_i} (where i indexes the peak). The
  period of a given population in each cycle is thus T^x_i ≡ t^x_i+1-t^x_i. 
  For a sufficiently long time series, we compute the mean period
  T_x and its variance.  
  
  In a similar way we can define the time delay between the sender and the receiver populations in each cycle
  τ_i=t^R_i-t^S_i (Fig. <ref>(a) and (b)). Then, if τ_i obeys a unimodal distribution,  we calculate τ as the mean
  value of τ_i and σ_τ as its variance.  
  
  If T_S  ≈  T_R
  and τ is independent of the initial conditions, the populations exhibit
  oscillatory synchronization with a phase-locking regime. In all those calculations we discard the
  transient time.  
  
  
  
  
  § RESULTS
  
  
  
  
  
  
   §.§ Phase-locking regimes: delayed and anticipated synzhronization
  
  
  
  
  
  
  In order to mimic the oscillatory activity in cortical regions 
  we simulated the sender population in such a way that the external noise 
  and the internal coupling are enough to allow the mean membrane potential to oscillate with f≃8 Hz (equivalent to T^S_i≃ 125 ms).
  Depending on the internal parameters of the receiver population, the sender-receiver coupling g_E can synchronize the activity of both areas or not. 
  The phase-locked regimes present non-zero phase-lags. 
  In Fig. <ref>(a) and (b) we show simulated time series of the S and R population. 
  Fig. <ref>(c) and (d) show the time delay τ_i in the i-th period as a function of the period index i 
  and Fig. <ref>(e) and (f) display their probability densities.
  
  The phase-lockings can be characterized by the mean time delay τ and its standard deviation.
  For sufficiently large g_E the mean time delay τ is positive which indicates that the sender population leads the receiver. 
  This is the usual delayed synchronization regime (DS).
  Left panels of Fig. <ref> show an example of DS for g_E=0.8 nS and g_I=0.02 nS. 
  For this set of parameters, the peak at the receiver population occurs on average ≈ 4.5 ms after the peak of the sender, 
  which is close to the magnitude of the synaptic time scales (as mentioned in Sec. <ref>: τ_E=5.26 ms)
  
  For larger values of inhibition g_I, the receiver population leads the sender on average, which is characterized by a negative mean time delay: τ<0. 
  This non-intuitive situation is the so-called anticipated synchronization <cit.>. 
  The right panels in Fig. <ref> exhibit an example of AS with g_E=0.5 nS and g_I=0.8 nS.
  For this set of parameters, the system exhibits τ≃-35.8 ms. Such value could not be inferred from any model parameter as, for example, the synaptic time scales. 
  This means that the anticipation time is a result of the nonlinear dynamics of the system.
  Furthermore, the anticipation does not occur every cycle, but the distribution of τ_i (see  Fig. <ref>(f)) clearly 
  shows that in the AS regime the majority of the peaks happens in the receiver-sender order. 
  
  
  
   §.§ Bistability between DS and AS
  
  
  
  
  
  
  
  A system that spends the vast majority of
  time near two well-separated regions exhibits a stationary density with two sharp peaks. 
  This phenomenon characterizes bistability.
  The Kramers problem is a standard example of a bistable system,
  which consists of diffusion in a double-well potential in the presence of noise.
  The distribution of the time spent in each region before a transition
  depends on the amount of noise.
  
  
  Here we show that depending on the relation between excitatory and inhibitory conductances the system can present a bistable regime between a DS and an AS regime (see Fig. <ref>). 
  The time delay τ_i is positive for a few cycles, with a well-defined mean value and standard deviation, which is similar to a DS regime for a certain amount of time. 
  Then, the system randomly switches for different dynamics in which τ_i is negative during a few other cycles. By analyzing the system only for these few periods
  of oscillation, one could wrongly characterize the system as in an AS regime. But, then, suddenly again, the system can jump back to the first attractor close to DS.
  Therefore, the histogram of time delays between the two populations is a bigaussian with one positive and one negative peak as shown in Fig. <ref>(c).
  In this regime, the system cannot be simply characterized by the mean time delay τ.
  
  
  If the system remains close to the DS-attractor (AS-attractor) for more than two cycles we define this as a DS-event (AS-event). 
  The DS-events are represented by the upper states in Fig. <ref>(b) whereas AS-events are the lower states in the figure.
  We define the event size as the number of cycles in which the system stays close to one of the two regions. 
  Since the mean period of oscillation of the sender is 125 ms we can use the size of the event as a proxy for the temporal dynamics of the bistable regime. 
  For example, an event that lasts for 8 cycles would have ≈1 s duration.
  This would be especially useful if one needs to compare the temporal dynamics of the model with behavioral data.
  
  In Fig. <ref> we show the size of the events and their histograms.
  The distribution of the number of occurrences of a specific size is different for DS-events and AS-events. 
  The probability to find very small events (up to 9 cycles) or very large events (larger than 300 cycles) is larger for DS than for AS-events. 
  On the other hand, events of intermediate sizes are more probable close to the AS region. 
  The distributions are qualitatively comparable to temporal dynamics of binocular rivalry during fMRI <cit.>.
  
  However, investigating in more detail the statistical properties of these distributions, to compare with cognitive data would require a significative computational effort to simulate even longer time series which is out of the scope of this study.
  
  
  
  
  
  
  
  
  
   §.§ Phase-drift regime
  
  
  
  
  
  
  
  For small values of the sender-receiver coupling g_E, the system can also exhibit a phase-drift (PD) regime in which the receiver is faster than the sender (T_S  > T_R).
  Fig. <ref> shows an example of such a regime:
  the time delay τ_i changes every cycle in a quasi-periodic configuration.
  
  The histogram of τ_i for g_E=0.3 nS seems flattened. 
  In the limit of g_E=0, which characterizes the totally uncoupled situation,
  every τ_i is equiprobable. 
  As in the bistable regime, in the PD, we do not use the mean time delay τ to characterize the regime.
  
  
  
  
  
  
  
  In Fig. <ref> we display a heatmap of the return map τ_i versus τ_i-1 which is useful to illustrate the different dynamical regimes.
  Several features in those curves are important to mention.
  Due to the noise, in every regime, there is a probability to find both positive and negative time delays per cycle τ_i. 
  However, in the DS regime, the probability density is clearly larger at the first quadrant, whereas in the AS the system remains for longer times in the third quadrant.
  Therefore, the bistable regime presents two denser regions: one in the first and the other in the third quadrant.
  In the PD regime there are continuous regions of the return map that are almost equally accessed by the system (see Fig. <ref>(d)). 
  This reflects the fact that the time delay varies from positive to negative in small steps during a few cycles (see Fig. <ref>(b)).
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§  Scanning parameter space
  
  
  
  The dependence of the system's dynamics on the synaptic conductances g_E and g_I is shown in
  Fig. <ref>. As one could expect, for small enough sender-receiver coupling there is no phase-locked regime. 
  As we decrease the g_E the system eventually reaches the PD.
  It is worth emphasizing that, similarly 
  to previous results on AS <cit.>, 
  the transition from DS to AS can be mediated by both internal inhibition and external noise at the receiver.
  However, differently from the previous studies <cit.>, here, the DS-AS transition
  does not occur via zero-lag synchronization,
  but through a bistable regime.
  
  To better understand the effects of external noise in the dynamical regimes, we vary the synaptic conductance g_P
  of the Poissonian input received by each neuron of the receiver population. This external input mimics synaptic currents received from  other 
  cortical regiones that are not included in our simple two-population model (see Sec. <ref> for more details).
  Fig. <ref> displays a two dimensional (g_I,g_P) projection of parameter space.
  The effect of increasing the noise at the receiver g_P is similar to the effect of decreasing the sender-receiver coupling g_E.
  In the PD regime, the time delay in each cycle varies almost continuously, 
  whereas in the bistable regime τ_i suddenly changes from positive to negative values in a non-predictable way. 
  
  
  Fig. <ref> shows illustrative examples of the time delay distribution as we change the external Poissonian input at the receiver population g_P 
  (along the vertical line  g_I=0.6 nS in Fig. <ref>. ). 
  By choosing conductances in such a way that the system presents AS when the amount of noise in both populations are the same (g_P=g_P^S=0.5 nS, g_E=0.5 nS and g_I=0.6 nS ), 
  as we increase the noise at the R population (g_P>0.5 nS) the system goes to a PD regime. 
  The AS-PD transition has been previously reported in neuronal microcircuits <cit.> but not in neuronal populations.
  On the other hand, by decreasing the noise at R (g_P>0.484 nS) the system undergoes an AS-DS transition via bistability. 
  
  We vary the Poissonian conductance g_P every 0.001 nS and the other synaptic conductances g_E and  g_I every 0.1 nS in order 
  to capture all the qualitative important features of the diagrams in Fig. <ref> and  <ref>
  The conditions to define the boundaries between the regimes in Figs <ref> and  <ref> have been determined 
  by analyzing the time delay distributions as in Fig. <ref>.
  For the DS regime the only required condition is a positive τ. 
  For the AS regime, besides a negative τ, the system should present an AS peak at least three times larger than the DS peak. 
  If this is not the case, the system could be in a bistable regime or a phase-drift.
  In the bistable regime the smallest peak should be at least 7 times higher than the probability to find intermediate values of τ_i between the peaks, otherwise it is a PD. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § CONCLUDING REMARKS
  
  
  To summarize, we have shown that a simple but biophysically plausible model 
  of two unidirectionally connected neuronal populations can present phase-bistability between anticipated and delayed synchronization.
  To the best of our knowledge, this the first verification of such a regime.
  Unlike previous studies on AS in neuronal models <cit.>, 
  here, the transition from DS to AS does not occur via zero-lag synchronization, but through the bistable regime.
  We have also shown that the interplay among the inhibitory conductance at the receiver, 
  the sender-receiver coupling and the external noise determines the dynamical regime.
  Moreover, for sufficiently large noise or small coupling the system can eventually reach a phase-drift regime in which the receiver population is faster than the sender. 
  
  
  Multi-stability of neuronal networks has been suggested as the mechanism underlying switches between
  different perceptions or behaviors <cit.>.
  Recently, multi-stability has also been associated with different oscillatory states of brain dynamics <cit.>.
  In particular, perceptual neuronal states models based on noise and adaptation have been used to qualitatively describe neurophysiological 
  experiments on human visual bistable perception <cit.>. This model alternates between two different active states 
  and reproduces probability distributions of dominance durations and their relation with the amount of noise. 
  However, these studies were not investigating phase relations during ambiguous perception.
  Therefore, we suggest that our results, using  populations of spiking neurons, could be potentially useful to study phase-bistability in cortical regions during bistable perception <cit.>.
  In fact, our model shows fixed structural connectivity that allows flexible dynamics which could change in timescales relevant for behavior. 
  As far as we know, this is the first spiking neuronal population model to present a bistable regime between two synchronized regimes with a positive and a negative phase-difference.
  
  
  
  Our results offer a number of possibilities for further investigation. 
  The DS-AS transition could possibly explain commonly reported short latency in visual systems <cit.>, 
  olfactory circuits <cit.>, songbirds brain <cit.> and human perception <cit.>.
  Future studies can be conducted in the light of understanding the functional significance of the diversity in the phase-relations between oscillatory brain regions <cit.>.
  Furthermore, including the effects of synaptic plasticity, especially spike-timing-dependent plasticity <cit.> 
  in the bistable regime is a natural next step which we are currently pursuing. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  The authors thank FAPEAL, UFAL, CNPq (grant 432429/2016-6) and CAPES (grant 88881.120309/2016-01) for finnancial support.
  
  
  
  
  
  
  
  