Dynamics of Time Delay-Induced Multiple Synchronous Behaviors in Inhibitory Coupled Neurons

The inhibitory synapse can induce synchronous behaviors different from the anti-phase synchronous behaviors, which have been reported in recent studies. In the present paper, synchronous behaviors are investigated in the motif model composed of reciprocal inhibitory coupled neurons with endogenous bursting and time delay. When coupling strength is weak, synchronous behavior appears at a single interval of time delay within a bursting period. When coupling strength is strong, multiple synchronous behaviors appear at different intervals of time delay within a bursting period. The different bursting patterns of synchronous behaviors, and time delays and coupling strengths that can induce the synchronous bursting patterns can be well interpreted by the dynamics of the endogenous bursting pattern of isolated neuron, which is acquired by the fast-slow dissection method, combined with the inhibitory coupling current. For an isolated neuron, when a negative impulsive current with suitable strength is applied at different phases of the bursting, multiple different bursting patterns can be induced. For a neuron in the motif, the inhibitory coupling current, of which the application time and strength is modulated by time delay and coupling strength, can cause single or multiple synchronous firing patterns like the negative impulsive current when time delay and coupling strength is suitable. The difference compared to the previously reported multiple synchronous behaviors that appear at time delays wider than a period of the endogenous firing is discussed. The results present novel examples of synchronous behaviors in the neuronal network with inhibitory synapses and provide a reasonable explanation.


Introduction
Synchronization is an emerging phenomenon of a population of interacting units in nature. Synchronization has attracted great interest in many different fields such as mechanics, physics, chemistry, biology, and ecology [1,2]. Synchronization in biological systems such as the nervous system [3][4][5] and cardiac myocytes [6][7][8] has been extensively explored. Synchronous behaviors have been observed in many central nervous systems such as the visual cortex, investigated. Period-K (K = 6 chosen as representative) bursting is investigated and the bursting period is labeled as T. Different from the results of Refs. [27,29] wherein only the first of the multiple synchronous behaviors appears at time delays within 0 < τ < T and each of the remaining behaviors appears at time delay within (n−1)T < τ < nT (n = 2, 3), multiple synchronous behaviors in the present paper appear at time delay within 0 < τ < T. The multiple synchronous behaviors are period-(K − 2), period-(K − 1), period-K and period-(K+1) bursting patterns, which appear at different intervals of τ when coupling strength is suitable. Using bursting dynamics acquired by fast-slow variable dissection method and the response of bursting pattern to a negative impulsive current at different phase of fast spike or quiescent state, the multiple synchronous behaviors are well interpreted. An isolated neuron manifests burst with K − 2, K − 1, K, and K+1 spikes when a negative impulsive current with suitable strength is applied suitable phase within (K − 2) th , (K − 1) th , K th spikes, and quiescent state of a burst, respectively. For a neuron in the motif, the inhibitory coupling current received by a neuron in the motif closely resembles the negative impulsive current applied to an isolated single neuron. The time delay and the coupling strength can modulate the application time and the strength of the inhibitory coupling current, respectively. If the time delay and coupling strength is suitable, different synchronous behaviors can be induced. These results present novel synchronous behaviors induced by inhibitory synapses with time delays.

Leech heart interneuron model of single neuron
In fact, the multiple synchronous behaviors can be simulated in the reduced model of Leech heart interneuron, Chay model and Rulkov map model. Considering that the reduced model of Leech heart interneuron has been widely used to investigate complex dynamics of neuronal firing patterns and the spatiotemporal dynamics of neuronal motifs or networks [25,26,29,34], in the present paper, Leech heart interneuron model is used as representative and is described as follows, C _ V ¼ À½ g Na f ðÀ150; 0:0305; VÞ 3 h Na ðV À E Na Þ þ g K2 m 2 K2 ðV À E K Þ þ g L ðV À E L Þ þ I pol ð1Þ _ h Na ¼ ½f ð500; 0:0325; VÞ À h Na =t Na ð2Þ Here, V is the membrane potential, h Na and m K2 are the gating variables describing the activation of potassium current and inactivation of sodium current, respectively. C is the membrane capacitance; g K2 and g Na are the maximum conductance of potassium current and sodium current, respectively; E K2 and E Na are the reversal potentials of K + and Na + , respectively; g L and E L are the conductance and reversal potential of leak current, respectively; I pol is the polarization current; τ Na and τ K2 are the time constants of activation of potassium current and inactivation of sodium current, respectively; f(x,y,z) = 1 / (1 + exp{x(y + z)}) is a Boltzmann function describing kinetics of the currents. V shift K2 is the shift of the membrane potential of half-activation of potassium current from its canonical value.
The parameter τ K2 (0.9 s) is much larger than τ Na (0.0405 s), therefore the model composed of Eqs (1-3) is a fast-slow system. The Eqs (1-3) and the third equation is the full system and slow subsystem (SS), respectively. The fast subsystem (FS) is the first two equations when m K2 is chosen as the bifurcation parameter.

Model of the motif with time delay
The motif is composed of three identical neurons described by the reduced model of the Leech heart interneuron. The connections are reciprocal inhibition coupling and bidirectional. The model of the motif composed of three Leech heart interneurons is described as follow: Here, h i and m i corresponds to h Na and m K2 , respectively, and I syn i is synaptic current received by neuron i (i = 1, 2, 3) and described as follows [29]: Here, i, j = 1, 2, 3, and j 6 ¼ i to avoid the self-coupling. E syn is the reversal potential of synaptic current. g ij is the coupling strength from neuron j to neuron i. E syn = −0.0625 V, which is set below the minimal value of V i (t) to ensure that the synapse is inhibitory. θ syn is set as −0.03 V to ensure that every spike within a burst can cross the threshold.

The synchronization index
A similarity function S ij , which has been used to investigate synchronization between neuron i and neuron j in previous studies [42], is employed in this paper and described as follows, Here, < > represent the time average, i = 1, 2, 3, j = 1, 2, 3, j i. S ij = S ji . If the signals V i (t) and V j (t) are independent, the difference between them is the same order as the signals themselves. If V i (t) = V j (t), as in the case of complete synchronization, S ij reaches its minimum value being as 0. In the present paper, S max = max(S 12 , S 23 , S 31 ) and the three neurons are regarded to achieve synchronization if S max < 10 −4 .

Parameter configurations
In the present paper, the coupling conductances of counterclockwise and clockwise connections are set as g 12 = g 23 = g 31 = g 1 , and g 21 = g 13 = g 32 = g 2 . Because of the coexistence of bursting patterns, the probability of S max < 10 −4 for 100 different initial values for each combination of τ, g 1 , and g 2 , which is labeled as R, is also used to characterize the synchronous dynamics.
burst are induced, and the relationships between the number of spikes per burst of burst pattern and application phase and strength of the negative impulsive current are built. by the dynamics of the endogenous bursting pattern of isolated neuron, which is acquired by the fastslow dissection method, combined with the inhibitory coupling current.
Synchronous behaviors with different bursting patterns can be simulated in inhibitory coupled neurons with time delay when different values are assigned to time delay and/or coupling strength. For a neuron in the motif, the inhibitory coupling current, of which the application time and strength is modulated by time delay and coupling strength, can play roles like the negative impulsive current. When time delay and coupling strength are suitable, single or multiple synchronous firing patterns can be induced.  remaining five spikes revolve around the UB from left to right and run along V max and V min . After the sixth spike, the trajectory runs across the MB through the neighborhood of SH and to first and then along the LB from right to left to form the quiescent state of the period-6 bursting. As the trajectory runs farther to the left and away from the SN bifurcation point, the first spike appears via a fast increase of V. The behaviors of period-6 bursting transit from quiescent state to spikes first via the SN (fold) bifurcation and then back to quiescent state via the homoclinic bifurcation. The bursting should be fold/homoclinic bursting [37,38].

Dynamics of period-6 bursting pattern of isolated neuron
The response and burst pattern of the a single neuron to a negative impulsive current For such a period-6 bursting, if a square "negative" pulse current with a fixed time width (labeled with ΔT) and strength (labeled by A) is introduced at a suitable time (Δt), which is called action time in this paper, is measured by the time beginning from the circle shown in Fig  1(A), and lies between the k th peak and k th trough (k = 3, 4, 5), a burst with k spikes is induced when A is larger than a threshold, and the firing pattern is still a burst with 6 spikes when A is smaller than the threshold. For example, the threshold of a square If action time Δt is changed from the k th peak to (k + 1) th trough, the threshold of negative current with ΔT = 0.133 s decreases first and then increases drastically, as shown by the lines with triangle (k = 3), circle (k = 4), and dot (k = 5) in Fig 3. The threshold increases as k is decreased from 5 to 3, as shown in Fig 3. If Δt is between k th trough and (k + 1) th peak (k = 3, 4, 5), i.e., the ascending branch, it is much more difficult to induce a burst different from the burst with 6 spikes than Δt between k th peak and k th trough, i.e. the descending branch.
When Δt of a square "negative" current with ΔT = 0.133 s is between the 6 th peak and 6 th trough, the firing is still bursting with 6 spikes per burst when A is chosen as a value from 0 to a large enough value. The threshold of the negative current to induce a burst with 6 spikes is zero. Examples with A = −0. When Δt is within the middle (1.70 s Δt 1.93 s) and end (1.94 s Δt 2.88 s) of the quiescent state, a square "negative" current with ΔT = 0.93 s can induce a burst with 7 spikes (solid line) when A is larger than a threshold, and the firing is still burst with 6 spikes when A is smaller than the threshold. The example of Δt = 1.70 s is shown in Fig 2(E1) and 2(E2) and the threshold is −0.079 mA. After the negative impulsive current, the membrane potential decreases to a large extent firstly, and then the trajectory of the first spike moves to left to form a burst with 7 spikes.
With increasing Δt, the change of the threshold of the negative current to induce a burst with 7 spikes decreases, as shown in Fig 3. The threshold within the middle sub-region of the quiescent state is relatively large and decreases. The threshold within the end sub-region of the quiescent state remains at a roughly fixed level, which is less than that to induce a burst with 5 spikes and is larger than that to induce a burst with 6 spikes when Δt is within a descending branch of the 6 th spike and within the beginning sub-region of the quiescent state.
With exception of the single "negative" impulse, a pair of "negative" impulses is also considered. If Δt of a subthreshold square "negative" impulse current is within k th spike and Δt of a suprathreshold current is within (k + 1) th spike (k = 3, 4, 5), a burst with k + 1 spikes can be evoked. For example, the behavior when k = 4 is shown in Fig 4. The doublet impulses, the burst with 5 spikes, and the period-6 burst are shown by dash-and-dot, dotted, and solid lines, respectively.

Multiple synchronous behaviors and the coexisting firing patterns of the motif
When values of parameter g 1 , g 2 , and τ are suitable, the motif can exhibit multiple synchronous behaviors and the coexisting firing patterns. For example, when g 1 = 0.6 nS, the distribution of S max values in the plane (τ, g 2 ) for two different initial values are illustrated in Fig 5(A) and 5 (B), respectively. The synchronous behaviors (S max < 10 -4 ) can be found in multiple ranges of 0 < τ < T, as shown by the arrows in Fig 5(A) and 5(B). The positions of arrows from left to right correspond to those of the lines shown in Fig 3. The multiple synchronous behaviors mainly appear when 0 < τ < T, which is different from those in Refs. [23,38], wherein each of multiple synchronous behaviors appear within the range of (n−1)T < τ < nT (n = 1, 2, 3,), and the only one synchronous behavior appears within the range of 0 < τ < T.  The difference can be found between Fig 5(A) and Fig 5(B). For example, the value of S max shown by the bold arrows in the upper and lower panels is approximately zero and 1, respectively, showing the coexistence of different behaviors at some τ values. When g 1 = 0.6 nS and g 2 = 0.9 nS, the changes of S max values with respect to τ for the initial values corresponding to Fig 5(A) and 5(B) are shown by the solid and dotted lines in Fig 5(C). Distinction between the two lines, i.e., the coexisting behaviors, can be found, as shown by two arrows in Fig 5(C). When τ = 0.65 s (left arrow), the dotted line presents a non-synchronous bursting pattern with S max = 0.645 (Fig 6), and the solid line presents the coexisting synchronous period-4 bursting with S max = 0 (Fig 7(A)).
The spike trains of the period-4 bursting patterns of three neurons illustrated in Fig 7(A) are shown by solid lines, which show synchronous behaviors. The coupling currents of the three neurons are shown by dotted lines in Fig 7(A). There was a negative impulse of the coupling current that is located near the 4 th spike and is suprathreshold to induce a burst with 4 spikes per burst, which is similar to that shown in Fig 2(A).
Except for the synchronous period-4 bursting, many synchronous bursting patterns can be simulated with different combinations of τ, g 1 , and g 2 values when 0 < τ < T, 0 < g 1 < 1, and 0 < g 2 < 1, as shown in Fig 7( Fig 7(B). This condition is similar to that of the double impulses shown in Fig 4. The former impulse is subthreshold and the latter is suprathreshold. When τ is lengthened to 2.18 s, g 1 = 0.2 nS, and g 2 = 0.16 nS, bursts with 7 spikes are induced by the suprathreshold "negative" current of the coupling current, as shown in Fig 7(C), which is similar to Fig 2(E). If the coupling strength is decreased to g 1 = 0.01 nS and g 2 = 0.2 nS, the coupling current becomes subthreshold and the synchronous bursting pattern with 6 spikes appears, as shown in Fig 7(D).

Definition of different synchronous bursting patterns
Based on the relationships between the synchronous firing pattern and the corresponding coupling current, P k l or P k qc is used to discriminate different synchronous patterns. For example, the synchronous bursting shown in Fig 7(A), 7(B), 7(C), and 7(D) are P 4 4 , P 5 4 , P 7 qe , and P 6 qe , respectively. For both P k l and P k qc , "k" presents the number of spikes per burst of the synchronous firing pattern (k = 4, 5, 6, 7). For P k l , "l" is an integer number, which shows that the delayed coupling "inhibitory" current is at l th spike of the burst (l = 3, 4, 5, 6). For P k qc , "q" means that the action time of delayed coupling current is within the quiescent state, and "c" is a character of Fig 5. Multiple synchronous behaviors and coexistence of firing patterns. The numbers 3, 4, 5, 6, and 7 represent burst with 3, 4, 5, 6, and 7 spikes, respectively. BQ, MQ, and EQ lying between two dashed lines represent the beginning, middle, and end of the quiescent state of period-6 bursting (bursting pattern of isolated neuron), respectively. (a) The dependence of S max on both g 2 and τ; (b) the dependence of S max on both g 2 and τ with initial values different from those of Fig (a); (c) the dependence of S max on τ when g 1 = 0.6 nS and g 2 = 0.9 nS. The dotted (solid) line corresponds to initial values of Fig 5(a) (Fig 5(b)). doi:10.1371/journal.pone.0138593.g005 Time Delay-Induced Multiple Synchronizations "b", or "m", or "e", corresponding to the beginning, middle and end sub-regions of the quiescent states. Therefore, the synchronous bursting patterns shown in Fig 6E-6J are P 5 5 , P 6 5 , P 6 6 , P 6 qb , P 6 qm , and P 4 3 , respectively, which closely match the examples analyzed in Fig 2. The dependence of the synchronous behaviors on coupling strength and time delay Multiple synchronous behaviors in detail and the relationships between different synchronous bursting patterns can be further found in the distributions of R in the two-dimensional parameter space (τ, g 2 ) at four different levels of g 1 = 0.001 nS, 0.01 nS, 0.2 nS, and 0.6 nS, as shown in  Fig 8(B), the synchronous behavior appears only once within 0 < τ < T, and the synchronous bursting pattern is P 6 qb . The results also show that single synchronous behavior is changed to multiple synchronous behaviors within 0 < τ < T with increasing coupling strength, as shown in Fig 8. The relationships between various synchronous patterns can be clearly found. For example, the coupling strength of P 4 4 , P 5 5 , and P 7 qe is larger than that of P 5 4 , P 6 5 , and P 6 qe . From left to right, the synchronous bursting patterns with 4, 5, 6, and 7 spikes per burst appear. The τ ranges of bursting patterns with 6 spikes and with 7 spikes per burst are much wider than bursting patterns with 4 and 5 spikes per burst, which is consistent with the width of the curves of the threshold values shown in Fig 3. The synchronous bursting patterns with 6, 7, 5, 4 and 3 spikes per burst appear in sequence with increasing g 2 , which is consistent with the values of the threshold curve shown in Fig 3. With increasing coupling strengths g 1 and g 2 , the probabilities R of the same synchronous patterns increase or remain nearly 1.0, as shown in Fig 8. The changes of R with respect to τ along the five dotted lines in Fig 8   . For P k l (l = 3, 4, 5 or 6, k = 4, 5 or 6) and P k qc (k = 6 or 7, "c" represents ''b'', "m" or "e", respectively), the numbers 3, 4, 5, 6, and 7 of superscript represent synchronous burst with 3, 4, 5, 6, and 7 spikes, respectively. The numbers 3, 4, 5, 6 of subscript show that action positions of coupling current are within 3th, 4th, 5th, and 6 th spikes, respectively. The characters ''b'', "m" or "e" of subscript show that action positions of coupling current are within beginning, middle, and ending part of the quiescent state of period-6 bursting, respectively.  The distribution of R on plane (τ, g 2 ) at different g 1 levels. The numbers 3, 4, 5, 6, and 7 represent bursting pattern with 3, 4, 5, 6, and 7 spikes, respectively. P k l (l = 3, 4, 5 or 6, k = 4, 5 or 6) and P k qc (k = 6 or 7, "c" represents ''b'', "m" or "e", respectively) represent synchronous bursting patterns shown in Fig 7.  The dependence of R on g 1 and g 2 at different τ levels The distribution of R values in the two-dimensional parameter space (g 1 , g 2 ) at different τ levels is helpful to understand the influence of the cooperation between g 1 and g 2 on the synchronization dynamics, as shown in Fig 10. Five levels of τ, τ = 0 s (Fig 10(A)), τ = 0.65 s (Fig 10(B)), τ = 0.83 s (Fig 10(C)), τ = 2.3 s (Fig 10(D)), τ = 1.45 s (Fig 10(E) and 10(F)), are investigated as representatives.
When τ = 0, most R values equal 0 and the remaining equal 0.01, indicating that the motif is difficult to achieve synchronization, as shown in Fig 10(A). When τ 6 ¼ 0, the R values exhibit symmetrical characteristics with g 1 = g 2 in the (g 1 , g 2 ) plane, showing that g 1 and g 2 play equivalent roles in influencing the synchronization dynamics, as shown in Fig 10B-10F. When τ = 0.65 s, τ = 0.83 s, and τ = 2.3 s, two separated regions are shown by arrows in Fig 10(B), 10(C), and 10(D), respectively. The lower left region is P 5 4 , P 6 5 , and P 6 qe , and right upper region is P 4 4 , P 5 5 , and P 7 qe , which is consistent with Fig 8. It should be noted that the motif exhibits synchronization when the sum of g 1 and g 2 are strong. When τ 1.45 s, there is only one synchronous bursting pattern (P 6 qb ), as shown in Fig 10(E) and 10(F) (Fig 10(F) is the enlargement of the downleft corner).
When comparing Fig 10(B)-10(D) with Fig 10(F), it can be found that the threshold of coupling strengths that can induce synchronous behaviors of P 6 q , P 7 qe , P 5 5 and P 4 4 increases in sequence, which is consistent with that of Fig 8. The multiple synchronous behaviors of the motif composed of neurons with other bursting patterns When the behavior of the neuron in the motif is chosen as period-5 bursting with V shift K2 = −0.008 V, the multiple synchronous behaviors appear within a period (2.778 s) of period-5 bursting of time delay at different coupling strength, as shown in Fig 11(A) (g 1 = 0.2 nS and g 2 = 0.16 nS), 11(B) (g 1 = 0.3 nS and g 2 = 0.18 nS), and 11(C) (g 1 = 0.6 nS and g 2 = 0.6 nS). If the bursting pattern of the motif is chosen as a bursting pattern within a bifurcation scenario with respect to V shift K2 (Fig 12, similar firing pattern has been investigated in Refs. [43,44]), the motif can exhibit multiple synchronous behaviors similar to those of the present paper (not shown here).

Discussion
Time delay induced multiple synchronous behaviors are simulated in inhibitory coupled neurons with endogenous bursting, which appear within a period of the endogenous bursting. This is different from the previous studies in which only one synchronous behavior appeared within a period of the endogenous firing, and the multiple synchronous behaviors appear when time delay is lengthened to a range of integer multiples of periods of the endogenous firing [27,29]. In the present paper, when coupling strength is relatively small, only synchronous period-6 bursting that corresponds to the first of the multiple synchronous behaviors in the previous The changes of R values with respect to τ at different coupling strengths (g 1 ,g 2 ). The numbers 3, 4, 5, 6, and 7 represent burst with 3, 4, 5, 6, and 7 spikes, respectively. BQ, MQ, and EQ represent the beginning, middle, and end of the quiescent state of period-6 bursting, respectively. P k l (l = 3, 4, 5 or 6, k = 4, 5 or 6) and P k qc (k = 6 or 7, "c" represents ''b'', "m" or "e", respectively) represent synchronous bursting patterns shown in Fig 7.  studies is induced at approximately half of the period of the endogenous bursting (Fig 8A and  8B). When coupling strength is relatively large, multiple synchronous behaviors are simulated when time delay is shorter than a period of endogenous bursting in inhibitory coupled neurons, which is a novel example of multiple synchronous behaviors. We also examined other busting patterns with the same dynamics as the endogenous periodic bursting, such as period-5 bursting, and found that time-delay induced multiple synchronous behaviors can be found within a period of endogenous bursting. It shows that the multiple synchronous behaviors are independent to number of spikes per burst, but dependent to bursting. We also investigated multiple synchronous behaviors when time delay is lengthened to (n−1)T < τ < nT (n = 2, 3) (not shown here). Although not shown in the present paper, we examine the multiple synchronous behaviors simulated in network composed of Chay model with fold/homoclinic bursting or Rulkov map model with bursting. However, the multiple synchronous behaviors appear at time delay within a period of endogenous spiking has not been simulated in network composed of Hodgkin-Huxley model with spiking pattern. The results show that multiple synchronous behaviors are not dependent to the kind of neuronal model, but dependent to the kind of bursting. All these show that the multiple synchronous behaviors investigated in this paper are common in network composed of bursting neurons.
Furthermore, the bursting patterns of the synchronous behaviors, which are also dependent on time delay and coupling strength, show a close relationship with the endogenous period-6 bursting pattern. For example, the inhibitory coupled current with strong strength can terminate the behavior of period-6 bursting at k th spike to form the period-k (k = 4, 5), which can be well interpreted by the dynamics of endogenous period-6 bursting acquired with the fast-slow variable dissection method. With respect to the increase of time delay, the synchronous bursting patterns investigated in the present paper show a period-adding bifurcation from period-4 to period-7, which is similar to the previous studies wherein synchronous dynamics of the inhibitory coupled chaotic Hindmarsh-Rose neurons with distributed time delay are investigated [33]. However, the cause of the synchronous period-k (k = 1, 2, 3, 4, 5, 6) bursting patterns simulated at different time delays were not addressed in the investigation [33]. We suggest that the generation of these period-k (k = 1, 2, 3, 4, 5, 6) bursting patterns investigated in Ref [33] may obey the same rule as the synchronous behaviors of the present investigation. The present paper presents an example that different bursting patterns of the multiple synchronous behaviors of motif can be well explained by the dynamics of the endogenous bursting pattern of isolated neuron.
In the previous studies, synchronous behavior can be induced by inhibitory synapses with zero time delay in case of fast synapse at special conditions as investigated in Refs. [25,26] and slow synapses [20][21][22][23][24] and inhibitory synapses with time delay [27][28][29][30][31][32][33]. The results of the present paper are different from these cases, presenting novel synchronous behaviors induced by inhibitory synapses with time delay. Our findings suggest that time delay and coupling strength are important for controlling or adjusting the novel synchronous bursting patterns. The results Time Delay-Induced Multiple Synchronizations provide an inhibitory synaptic mechanism responsible for the generation of multiple synchronous behaviors in the nervous system. The period-6 bursting is identified as fold/homoclinic, which is a very common bursting pattern observed in different realistic nervous tissues or simulated neuronal models [39][40][41]. In addition, the motif composed of two coupled neurons with slow inhibitory synapses can achieve synchronous behavior, which has been demonstrated in the biological experiment [24]. The synchronous behavior is similar to the single synchronous behavior investigated in the present paper. If the coupling strength is increased, it is expected that the two coupled bursting neurons can achieve multiple synchronous behaviors, which should be demonstrated in biological experiment in future.