## Abstract

Synchronized cortical activities in the central nervous systems of mammals are crucial for sensory perception, coordination and locomotory function. The neuronal mechanisms that generate synchronous synaptic inputs in the neocortex are far from being fully understood. In this paper, we study the emergence of synchronization in networks of bursting neurons as a highly non-trivial, combined effect of electrical and inhibitory connections. We report a counterintuitive find that combined electrical and inhibitory coupling can synergistically induce robust synchronization in a range of parameters where electrical coupling alone promotes anti-phase spiking and inhibition induces anti-phase bursting. We reveal the underlying mechanism, which uses a balance between hidden properties of electrical and inhibitory coupling to act together to synchronize neuronal bursting. We show that this balance is controlled by the duty cycle of the self-coupled system which governs the synchronized bursting rhythm.

This article is part of the themed issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.

## 1. Introduction

Neuronal synchrony has been shown to be central to the function and dysfunction of human cognitive processing, memory and locomotion [1,2]. Synchrony plays a positive role in hippocampal networks and its disruption due to traumatic brain injury has been shown to severely impair cognitive processing and memory function for many years post-injury [3]. At the same time, synchronized neuronal firing is notoriously known to induce pathological brain states, especially during epilepsy and Parkinson’s tremors [4–7]. In particular, epilepsy is widely considered a dynamical network disease and is characterized by short bursts of synchronized neuronal activity and long events called seizures. This abnormal synchrony either is localized in a specific area of the brain, yielding a simple focal seizure, or spreads across the whole brain region, usually paralysing a patient and resulting in a complex generalized seizure [7]. Although there have been considerable advances in the treatment and understanding of the origin of epileptic seizures, the question of why the vast regions of brain become excitable and susceptible to synchronization remains open.

The emergence of synchronized rhythms in simple and complex networks of spiking and bursting neurons has been extensively studied in the literature [8–20]. Bursting occurs when neuron activities alternate between a period of quiescence and fast repetitive spiking [21–23]. There is experimental evidence that epileptic seizures are accompanied by changes in neuronal bursting activities [24,25] where individual spikes play an important role. In contrast with spiking neurons, whose synchronous behaviour is quite simple, coupled bursting neurons are capable of generating various forms of neuronal synchrony. These include synchronization of individual spikes, burst synchronization when only the envelopes of spikes synchronize while the spikes remain unlocked, and complete synchronization. The onset of neuronal synchrony is controlled by a non-trivial interplay between the intrinsic dynamics of individual neurons, the type of synaptic connections and network architecture.

Electrical and synaptic connections often play different roles in inducing synchronization or anti-phase spiking and bursting [26–32]. In contrast with slow or delayed inhibitory connections that favour neuronal synchrony [33–37], fast non-delayed inhibition is known to promote pairwise anti-phase bursting in a network with purely inhibitory synapses [8]. This is always the case in a pair of spiking neurons with fast non-delayed inhibitory connections, unless each neuron has more than one slow intrinsic variable [35]. It was also demonstrated that weak fast non-delayed reciprocal inhibition can favour the coexistence of in-phase and anti-phase bursting in networks of some bursting cell models; however, the in-phase rhythm is fragile and has a small basin of attraction [38,39]. In a recent work, we have shown that the addition of pairwise repulsive inhibition to excitatory networks of bursting neurons can induce synchrony, in contrast with one’s expectations [40,41]. This synergistic effect originates from the transition between different types of bursting, caused by strong excitatory–inhibitory synaptic coupling.

Many experimental findings indicate the presence of electrical coupling in GABAergic interneurons in the central nervous systems of mammals [42], particularly among neocortical neurons of the same class [43]. Networks of low-threshold-spiking neocortical interneurons with fast inhibitory synapses were found to be connected locally by synchronizing electrical coupling, a phenomenon which may be central to the coordination of strongly connected cortical subnetworks [44]. Indeed, GABAergic networks in the neocortex are known to control spike timing and influence rhythmogenesis throughout the entire neocortex despite a relatively small number of such inhibitory neurons [45]. Notably, it was shown that both GABA inhibitory currents and gap-junctional coupling are required for synchronized bursting in hippocampal interneurons of the rat [26]. The role of the duty cycle, the fraction of the period during which the neuron bursts, in promoting anti-phase bursting in networks with pairwise inhibitory and gap-junctional connections was previously discussed [46]. It was demonstrated that a short duty cycle can destabilize anti-phase bursting in an inhibitory network, but the addition of electrical coupling can restabilize the anti-phase pattern [46].

In this paper, we contribute to further understanding of cooperative dynamics in networks of bursting neurons with both gap-junctional (electrical) and inhibitory connections. We report a non-trivial synchronization mechanism of the combined coupling where electrical and inhibitory connections can synergistically induce synchronized bursting in a range of parameters where electrical coupling alone promotes anti-phase spiking and inhibition induces anti-phase bursting. The synchronization mechanism, where ‘two wrongs make a right’, is based on the properties of (i) weak electrical coupling to stabilize burst but destabilize spike synchronization and (ii) inhibition to generally promote anti-phase bursting but stabilize spike synchrony when initial conditions are close. The combined action of the two couplings uses the best of the two worlds to foster synchronized bursting, provided that a balance between the repulsive and attracting components of the combined coupling is preserved. Through analysis and numerics, we demonstrate that this balance is controlled by the duty cycle of the self-coupled system which governs the synchronized bursting rhythm.

This synergistic synchronization effect differs from the ones previously observed in networks with combined chemical and electrical synapses [30–32]. More specifically, it was shown [30] that a small amount of electrical coupling added to already significant inhibition of the network can increase synchronization more than a very large increase in the synchronizing inhibitory coupling. Notably, each kind of synapse in this network setting [30] alone fosters synchrony, but the resultant effect is much more pronounced. It was also demonstrated [31,32] that combining electrical synapses with inhibition in a network of spiking cells can enhance synchrony, whereas electrical synapses alone would impede synchronization. For this property to be true, the coupling strength of both electrical and chemical synapses should be sufficiently strong. In this setting, electrical and inhibitory synapses may both foster synchrony, or may compete, with one being an attractive force while the other repulses the cells. By contrast, the synergistic effect reported in this paper arises from nonlinear interaction of electrical and chemical synapses in a range of coupling strengths where both synapses alone impede complete synchrony. The discovered synergistic effect is due to nonlinear interactions of bursting cells at the level of bursts and spikes and is not observed in networks of spiking cells or reduced phase models of neurons. In this regard, our study along with our previous work [40,41] promotes the use of the detailed biophysical models that take into account neuronal spikes and bursts.

The layout of this paper is as follows. First, in §2, we describe the individual neuron and network models. In §3, we report the synchronization effect observed in the simplest two-cell network with electrical and inhibitory connections. We also use Poincaré maps to demonstrate how the number of coexisting phase-locked states changes as a function of the inhibitory coupling strength. Then, we employ a slow–fast decomposition of the networked system to isolate the impacts of the electrical and inhibitory coupling on the emergence of synchronized bursting. In §4, we introduce the variational equations for the stability of the synchronous bursting solution and explain the main synchronization mechanism via the calculations of averaged synaptic terms and their dependence on the duty cycle of synchronous bursting. In §5, we demonstrate that the synergistic effect is also present in larger networks and identify network topologies with the highest and lowest resilience of synchronized bursting. In §6, a brief discussion of the obtained results is given. Appendix A contains a rigorous proof of the stabilizing role of strong electrical coupling in synchronization of bursting Sherman cells. Finally, appendix B describes the numerical methods used in our study.

## 2. The network model

We consider a network of *N* Hodgkin–Huxley-type neuronal models [47] with electrical and inhibitory synapses:
2.1The intrinsic dynamics of the *i*th cell is represented by the membrane potential *V*_{i}, and the gating variables *n*_{i} and *S*_{i} are the opening probabilities of the fast and slow potassium currents, respectively. Function *F*(*V*_{i},*n*_{i},*S*_{i})=−[*I*_{Ca}(*V*_{i})+*I*_{K}(*V*_{i},*n*_{i})+*I*_{S}(*V*_{i},*S*_{i})] describes three intrinsic currents: fast calcium, *I*_{Ca}, persistent potassium, *I*_{K}, and slow potassium, *I*_{S}, currents such that
According to the Hodgkin–Huxley formalism, the steady-state values for the activation and inactivation of the fast and slow currents are represented by the Boltzmann equations as functions of *V*_{i},
Other intrinsic parameters are chosen and fixed as follows: *τ*=20 (ms), *τ*_{S}=10 000 (ms), (nS), *E*_{Ca}=25 (mV), (nS), *E*_{K}=−75 (mV) and (nS). The individual unit of the network (2.1), the Sherman cell model [47], was originally introduced to mimic the electrical activity of a pancreatic β-cell. This model is known to exhibit different types of bursting such as square-wave, plateau and pseudo-plateau bursting [48], and is often used as a generic Hodgkin–Huxley-type model to describe neuro-computational properties of bursting neurons and networks [23]. In the given set of parameters, the uncoupled cell generates square-wave bursting [21] (figure 1). The presence of the large parameter *τ*_{S}=10 000 (mV) makes the system (2.1) slow–fast such that the (*V*,*n*)-equations represent the two-dimensional fast ‘spiking’ subsystem; the *S*-equation corresponds to the slow one-dimensional ‘bursting’ system. The dynamics is centred around nullcline of the fast (*V*,*n*)-subsystem. The intersection between and the nullcline of the slow *S*-subsystem yields a saddle fixed point [40].

The cells are identical, and the coupling strength of the electrical (*g*_{el}) and inhibitory (*g*_{inh}) synapses is uniform for each type of coupling. The electrical coupling between cells *i* and *j* is modelled via the difference between the membrane potentials *V*_{i} and *V*_{j}. In order to make the chemical synapse inhibitory, the reversal potential is chosen at the level *E*_{inh}=−75 (mV), such that *E*_{inh}<*V*_{i}(*t*) for all permissible values of *V*_{i}. The inhibitory coupling is instantaneous and non-delayed; a smooth approximation of the Heaviside function is used to model the synaptic coupling function known as the fast threshold modulation [9]. The synaptic threshold *Θ*_{s}=−40 (mV) is chosen such that spikes in the single cell burst can cross the threshold (figure 1). Therefore, a spike in presynaptic cell *j* activates the synaptic current entering postsynaptic cell *i* (via *Γ*(*V*_{j}) switching from 0 to 1).

In (2.1), *N*×*N* connectivity matrices ** C**=(

*c*

_{ij}) and

*D*=(

*d*

_{ij}) describe the network structure of the electrical and inhibitory synapses, respectively. The electrical coupling matrix

**=(**

*C**c*

_{ij}) is symmetric as the electrical coupling is always undirected such that

*c*

_{ij}=

*c*

_{ij}and

*c*

_{ij}=1, if neuron

*i*receives an input from neuron

*j*. The nodes of the electrical network may have different in-degrees and receive a different number of inputs. The inhibitory coupling matrix

**can be asymmetric such that both mutual and unidirectional couplings are allowed. As in matrix**

*D***,**

*C**d*

_{ij}=1 if neuron

*i*receives an input from neuron

*j*; however,

*d*

_{ii}=0. We require the connectivity matrix

**to have all row-sums equal to**

*D**k*

_{inh}. This property implies that each cell on the inhibitory network receives

*k*

_{inh}inputs from other cells and this number is uniform for each cell. This requirement is a necessary condition for the existence of the synchronization subspace

*M*={

*V*

_{1}=⋯=

*V*

_{N}=

*V*(

*t*),

*n*

_{1}=⋯=

*n*

_{N}=

*n*(

*t*),

*S*

_{1}=⋯=

*S*

_{N}=

*S*(

*t*)} that defines complete synchronization between the cells. The dynamics of completely synchronized cells is governed by the following system: 2.2It is worth noticing that the synchronous behaviour differs from that of the uncoupled cell with

*g*

_{inh}=0 and

*g*

_{el}=0 due to the presence of the additional inhibitory synaptic term. As the electrical coupling disappears when

*V*

_{i}=

*V*

_{j}, the electrical coupling term is not present in system (2.2). As a result, changing the strength of inhibitory coupling can change the synchronous dynamics. In the following, we will show that these changes, induced by moderately weak inhibitory coupling, result in small variations of the duty cycle of synchronous bursting in system (2.2) and lead to stable synchronization.

## 3. Tug-of-war synchronization effect of combined coupling

We begin with the simplest network (2.1) where two cells are symmetrically coupled through electrical and inhibitory connections with *k*_{inh}=1. We will study this two-cell network to reveal the synergistic effect of combined coupling and describe its stability mechanism. We will then demonstrate that this effect is also present in larger networks and discuss the role of network structure.

### (a) Multistability and emergent synchronized bursting

Figure 2 demonstrates the onset of synchronized bursting in the two-cell network as a function of the electrical (*g*_{el}) and inhibitory (*g*_{inh}) coupling strengths. Figure 2*a*,*b* indicates that a strong electrical coupling, exceeding a threshold value *g*_{el}≈0.18, synchronizes the cells in the absence of inhibition (*g*_{inh}=0). The addition of inhibition to the strong electrical coupling impedes complete synchrony, as one would expect, and gradually increases the threshold value of electrical coupling *g*_{el}. Appendix A contains a rigorous derivation of an upper bound for the electrical coupling threshold, required for stable synchronization in the absence of inhibition. This analytical bound is very conservative and yields the synchronization threshold at *g*_{el}=3.925 compared with the actual threshold *g*_{el}≈0.18 (figure 2). However, it rigorously proves that the electrical coupling always promotes synchronization when it reaches the threshold value. Surprisingly, there is a range of much weaker electrical coupling *g*_{el}∈(0,0.02) (figure 2*c*,*d*) where the electrical coupling alone always impedes synchronization but the addition of inhibition can yield complete synchrony. The fact that increasing the electrical coupling within the range *g*_{el}∈(0,0.02) makes the synchronization solution more unstable is verified via the calculation of the largest transversal Lyapunov exponent of the synchronous solution which is positive and monotonically increases within the interval *g*_{el}∈(0,0.02) (figure 3) (the details on the calculation of the transversal Lyapunov exponent are given in appendix B). Note that, once the electrical coupling becomes stronger and lies beyond this coupling interval, its repulsive force becomes attractive and promotes synchrony, ultimately stabilizing complete synchrony at the threshold value *g*_{el}≈0.18.

Figure 2*c*,*d* also demonstrates that inhibition alone can foster or destabilize complete synchrony, depending on the coupling strength and initial conditions. When one cell is initially in the spiking phase, and the other is in the quiescent inactive state, inhibition also impedes synchrony and promotes anti-phase bursting in the absence of electrical coupling (see the dark grey (red) colour area adjacent to the *g*_{inh}-axis in figure 2*c*). When both cells start close to each other, inhibition can promote complete synchrony via the mechanism of nonlinear interaction between the spikes, described in [38,39]. For this property to be true, the inhibitory coupling must be weak such that *g*_{inh}∈(0,0.008) (see the black (blue) tongue-shaped region adjacent to the *g*_{inh}-axis in figure 2*d*). Further increase of *g*_{inh} makes the inhibition desynchronizing, independent from the initial conditions. Remarkably, the combination of the electrical and inhibitory coupling, where each synapse alone impedes complete synchrony, can promote synchrony, regardless of the initial conditions (see the black (blue) areas in figure 2*c*,*d*), even though the synchronization effect is much more pronounced when the cells start from close initial conditions. Figure 2*d* illustrates the synergistic effect, when ‘two wrongs make a right’, at point *C*, which corresponds to the combined attractive action of the repulsive electrical and inhibitory coupling. Note the instability of synchronization at points *B* and *D*, where the electrical and inhibitory coupling alone destabilizes complete synchrony. We are especially interested in the transition from point *B* via point *C* to point *E*. This transition along the horizontal line *g*_{el}=0.01 corresponds to a route where the repulsive electrical coupling first competes with the weak synchronizing inhibition (within the range *g*_{inh}∈(0,0.008)), then acts synergistically with the repulsive inhibition to promote synchrony (point *C*), and finally cooperates with the repulsive inhibition in a linear fashion to promote anti-phase bursting (point *E*). Similarly, the transition from *D* via point *C* to point *F* along the vertical line *g*_{inh}=0.01 is accompanied by the transition from anti-phase bursting at point *D* to complete synchrony at point *C*, and back to out-of-phase bursting at point *F*. Here, increasing *g*_{el} from 0 makes the electrical coupling alone more repulsive (figure 3); yet, it yields a region of stable synchrony when combined with the (repulsive) inhibition at *g*_{inh}=0.01. Thus, the same forces can switch their stabilizing and destabilizing roles, similar to playing tug of war. In addition to the coexistence of complete synchrony and anti-phase bursting, the combined coupling can also induce multiple coexisting phase-locked states, as shown in figure 4.

Our goal is to explain this counterintuitive synergistic effect and reveal the properties of the coupled system (2.1) which make the combined coupling attractive. Towards this goal, we will first use Poincaré maps for the phase differences between two interacting cells to reveal the existence of multistable phase-locked states and their dependence on the strength of electrical and inhibitory coupling (see §8 for details of how the phases are introduced and calculated). Figure 4*a* illustrates how the phase differences between two cells stabilize after 40 bursts. The number of bursts to skip (*k*=40) is chosen large enough to avoid transient stages. Note that the *B*−*C*−*E* transition (see figure 2) originates from the phase-locked state where the electrical coupling alone tends to establish burst synchrony; at the same time, it promotes anti-phase spiking (figure 4*a*, (*B*)). Increasing the inhibition up to point *C* (see figure 2) helps to synchronize the spikes within the bursts (figure 4*a*, (*C*)). Further increase in *g*_{inh} up to point *E* destabilizes complete synchrony and establishes a phase-locked state close to anti-phase bursting (figure 4*a*, (*E*)). Remarkably, while being destabilizing for complete synchronization when acting alone, the impact of electrical and inhibitory coupling is different. The electrical coupling promotes (destabilizes) burst (spike) synchrony, whereas inhibition does the reverse and fosters spike synchrony and impedes burst synchrony. This suggests how the combined action of both types of coupling can stabilize complete synchrony. To further validate this observation and isolate the impact of the electrical and inhibitory coupling, we will use a slow–fast decomposition of the two-cell coupled system (2.1).

### (b) Insight from the slow–fast decomposition

To better understand the action of electrical and inhibitory coupling on the dynamics of two cells during two distinct phases of active spiking and quiescence, we employ the slow–fast property of square-wave bursting and dissect the network dynamics into the fast and slow components (figure 5). We choose and fix the slow variable *S* at some level *S*=0.18, which corresponds to the middle of the spiking phase (see figure 1). This yields the fast (*V*_{i},*n*_{i}) systems (*i*=1,2), which are coupled via the electrical and inhibitory coupling and mimic the interaction between the cells during the spiking phase of both cells. Similarly, decreasing the time constant *τ* for the second variable *n*_{i}, we effectively get rid of all the spikes and turn the coupled system into a two-relaxation-oscillator network. This network aims at reproducing the cooperative dynamics of the full system (2.1) during the stage when one cell is active while the other is inactive. In this setting, the active cell keeps the inhibition on such that the inactive (inhibited) cell is kept at the inactive state as long as the active cell is in the spiking phase, causing anti-phase bursting.

This slow–fast decomposition reveals striking differences between the impacts of the electrical and inhibitory synapses on synchronization in the networks of fast and slow subsystems (figure 6). In the given range, the electrical synapses always repel the spikes and attract their envelopes (bursts) (compare the circle from figure 6*a*,*c* with the triangle from figure 6*b*,*d*). At the same time, the inhibitory connections bring the spikes together but push the bursts apart (compare the square with the diamond from the diagrams). While the heatmaps of figure 6 can slightly differ, depending on the value of *S* (not shown), they remain qualitatively the same and indicate the same effect. Combined together, the two seemingly counter-actions of electrical and inhibitory synapses make up a rich multistable pattern in the full system and induce the synchronization mechanism that we have called a ‘tug of war’.

## 4. Stability mechanism: why does the duty cycle matter?

To better quantify the stability mechanism and reveal the property of the coupled system that controls the stability of complete synchronization, we use the stability equations for the infinitesimal transverse perturbations Δ*V* =*V*_{1}−*V*_{2}, Δ*n*=*n*_{1}−*n*_{2}, Δ*S*=*s*_{1}−*s*_{2} [13]:
4.1where *Ω*(*V*)=*S*_{1}+*S*_{2} with *S*_{1}=*Γ*(*V*) and *S*_{2}=(*E*_{inh}−*V*)*Γ*_{V}(*V*) is due to the contribution of the inhibitory coupling. Here, {*V* (*t*),*n*(*t*),*s*(*t*)} denotes the synchronous solution which corresponds to the self-coupled system (2.2), *Γ*_{V}(*V*) is the partial derivative of *Γ*(*V*) with respect to *V*. The stability of the completely synchronous solution corresponds to the zero fixed point {Δ*V* =0,Δ*n*=0,Δ*S*=0} of the variational equations (4.1). The function *Ω*(*x*) promotes the stability of synchronization when it becomes positive and has a destabilizing impact when it is negative [13].

The two terms *S*_{1} and *S*_{2}, composing *Ω*(*V*), play opposite roles in stabilizing synchronization. The first (stabilizing) term *S*_{1}≥0 remains turned on when the voltage *V* (*t*) is above the synaptic threshold *Θ*_{s}. The second (destabilizing) term *S*_{2}≤0 contains the derivative *Γ*_{V}(*V*), which has a negative peak around *Θ*_{s} (in the case of the Heaviside function, *Γ*_{V}(*V*) turns into the negative delta function). Hence, the term *S*_{2} switches and remains on for the values of *V* close to the threshold *Θ*_{s} when the spikes cross the threshold (figure 7). Thus, the terms *S*_{1} and *S*_{2} compete with each other to stabilize and destabilize the completely synchronous rhythm.

The contribution of the electrical coupling to the stability of the variational equations is always favourable due to the negative term −2*g*_{el}Δ*V*. As we seek to quantify the *B*−*C*−*E* transition (see figure 2) where the electrical coupling is fixed, we study the changes in the overall dynamics of the variational equations (4.1) as a function of the averaged contribution of the inhibitory synaptic terms *S*_{1} and *S*_{2} (figure 8). However, it is important to emphasize that, once the phase difference between the cells is no longer infinitesimal, the variational equations lose their credibility. As a result, the role of the electrical coupling for non-infinitesimal voltage differences cannot be assessed from the variational equations. As the phase map and slow–fast decomposition analysis suggest, this role is destabilizing for spike synchrony.

It is also important to stress that increasing the inhibitory coupling strength *g*_{inh} from 0 along the *B*−*C*−*E* route changes the dynamics of the self-coupled system (2.2) and alters the duty cycle of synchronous bursting in a nonlinear fashion (figure 8). This change turns out to be the critical quantity which shifts the balance between the competing terms *S*_{1} and *S*_{2}. As a result, a shorter duty cycle maximizes the averaged contribution of the resultant force *S*_{1}+*S*_{2} and induces complete synchronization.

## 5. Larger networks

The combined effect of electrical and inhibitory coupling is also present in larger networks (2.1). Figure 9 presents stability diagrams for synchronization in four-cell networks with different network structures of electrical and inhibitory connections. Our previous results on synchronization in excitatory–inhibitory networks [41] suggest that the structure of the added inhibitory connections is not important and only the number of inhibitory inputs controls the onset of synchronization, independent from all other details of their network topology. However, this is only true if the synchronizing excitatory connections are strong enough and form a connected graph which involves all the cells [41].

As figure 9 indicates, the sparse network topology with both local electrical and inhibitory connections (figure 9*d*) has the maximal horizontal and vertical size of the stability region. Notice that each cell in this network receives two inhibitory inputs such that *k*_{inh}=2. Similarly to the above-mentioned scaling law in excitatory–inhibitory networks [41], the horizontal size of the stability region in network configurations, where the combined effect is observed (the three networks in figure 9*a*–*c*), is inversely proportional to the number of incoming inhibitory connections *k*_{inh}. For example, the network with both global electrical and inhibitory connections (figure 9*a*) has the stability region whose horizontal size scales down by a factor to offset the effect of increasing the number of inhibitory inputs, *k*_{inh}, from two (as in the locally connected network with the maximal stability region) to three (as in the fully connected four-cell network). This scaling law originates from the self-coupled system (2.2) which governs the synchronous rhythm via the term *k*_{inh}*g*_{inh}(*E*_{inh}−*V*)*Γ*(*V*), whose impact remains the same as long as the quantity *k*_{inh}*g*_{inh} is preserved.

At the same time, the interplay between network structures of electrical and inhibitory connections and its impact on the stability of synchronization are highly non-trivial. Figure 9 demonstrates that global electrical connections should be compensated for by global inhibitory connections to enlarge the stability region (see figure 9*a*,*c*, network configurations). Note that the figure 9*b*,*e* networks do not exhibit the combined synchronizing effect in the region where inhibition is repulsive. Indeed, only two electrical connections in the figure 9*b*,*e* network can not burst synchronize all four cells, such that repulsive inhibition induces anti-phase bursting between the left and right sides of the network. The global structure of electrical connections in the figure 9*d* network does induce burst synchrony (not shown); however, the sparse directed inhibitory coupling is insufficient to overcome the impact of electrical coupling and synchronize the spikes. As a result, the combined effect cannot be achieved.

To show that the combined effect of electrical and inhibitory coupling appears in neural networks with complex topologies, we have simulated a 30-cell random network where each cell receives eight inhibitory connections (*k*_{inh}=8), whereas the number of electrical connections varies from one cell to another (figure 10). The uniform node in-degree of the inhibitory connections (*k*_{inh}=8) is preserved to guarantee the existence of the completely synchronized solution. Node degrees of the electrical connections do not affect this condition, and, therefore, were chosen freely. Figure 10*b* demonstrates the emergence of stable synchrony as a result of the synergistic interaction between the electrical and inhibitory coupling. Note that the size of the stability zone (black (blue) region) has shrunk, compared with the four-cell networks of figure 9. This is due to the increased node-degree of the inhibitory network and the above-mentioned scaling law.

The detailed analysis of the interplay between synchronization and the network structure of electrical and inhibitory connections is beyond the scope of this paper and will be reported elsewhere. This analysis can be based on the variational equations, similar to (4.1), and the application of the Connection Graph method [49–51], which uses Lyapunov functions and graph theoretical reasoning.

## 6. Conclusion

We have discovered a highly nonlinear effect of combined electrical and chemical synapses in promoting synchronization of bursting cells in a parameter region where each type of synapse alone destabilizes synchronization. This unexpected effect where ‘two wrongs make a right’ is caused by a sudden decay in the duty cycle of synchronous bursting. This change can induce stable synchronization as a result of the separable and counter-balancing effect of both coupling types on the slow (bursting) and fast (spiking) subsystems, corresponding to potassium and calcium ion channels. More precisely, fairly weak electrical coupling stabilizes burst synchronization but repels the individual spikes, whereas the inhibition does the opposite, promoting spike synchrony when the phases of the cells are close to each other and destabilizing burst synchronization when one of the cells is in an inactive state. The duty cycle controls a fragile balance between the two opposite forces, such that shorter duty cycles with a longer quiescent period increase the stabilizing impact of the electrical coupling in establishing burst synchronization. By the same token, these short duty cycles maximize the impact of inhibition in stabilizing spike synchrony. This dependence is non-trivial and increasing the inhibitory strength changes the duty cycle of synchronous bursting via the self-coupled system in a nonlinear fashion. The observed combined effect is not limited to networks of bursting Sherman cells but is also present in, for example, coupled Purkinje neuron models [52], capable of generating square-wave bursting. Our preliminary studies also indicate that synchronized bursting induced by the combined coupling persists under small parameter mismatch, including the intrinsic parameters and coupling strength.

Our study reinforces previous work [41], in which it was shown that the addition of strong pairwise repulsive inhibition to excitatory networks of bursting neurons can induce synchrony due to the transition between different types of bursting. Remarkably, the addition of the inhibitory coupling can promote synchronization much more significantly than strengthening the present excitatory connections. In contrast with excitatory–inhibitory networks studied in [41] where the excitatory connections are synchronizing, the combined effect of electrical and inhibitory coupling reported in this paper originates from two types of connections which are both repulsive. Our studies of neuronal synchronization form a basis for understanding the counterintuitive dynamics of bursting networks, which may yield meaningful insight into the phenomenon of pathological synchrony in epileptic networks. Epileptic seizures are strongly associated with a synchronized state of certain brain networks. Our results together with [41] suggest that promoting normally repulsive inhibition in an attempt to prevent seizures can have an unintended effect of inducing pathological synchrony.

## Authors' contributions

R.R. contributed to stability analysis and performed simulations. K.D. performed simulations. I.B. conceived the project and provided guidance. All authors contributed to writing the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the US Army Research Office Network Science Division under grant no. W911NF-15-1-0267 and the National Science Foundation (USA) under grant no. DMS-1616345. R.R. acknowledges support as a Georgia State University Brains & Behavior Program Fellow.

## Appendix A. The synchronizing role of strong electrical coupling

In this appendix, we rigorously prove that strong electrical coupling never changes its synchronizing role as long as it exceeds a synchronization threshold. We derive a rigorous upper bound on the strength of electrical coupling sufficient to induce globally stable synchronization in the two-cell network (2.1) in the absence of inhibitory connections.

### Theorem A.1.

*Complete synchronization in the network (*2.1*) of two mutually coupled cells with only electrical synapses is globally asymptotically stable if g*_{el}*≥g***, where* .

### Proof.

It can be seen from system (2.1) that each single cell with or without connections has an absorbing domain: 0≤*n*,*S*≤1,*E*_{K} ≤*V* ≤*E*_{Ca} such that any trajectory will eventually converge to this domain. Henceforth, we can assume that all the values of our system variables are inside this absorbing domain.

Similar to (4.1), we introduce the differences Δ*V* =*V*_{1}−*V*_{2}, Δ*n*=*n*_{1}−*n*_{2}, Δ*S*=*s*_{1}−*s*_{2}. As we target the global stability of synchronization, these differences do not have to be infinitesimal as in (4.1). Therefore, we obtain the following difference equation system from system (2.1) with *g*_{inh}=0:
A 1To have the explicit presence of Δ*V*, Δ*n* and Δ*S*, we apply the mean value theorem such that
A 2where and Strictly speaking, the values of and in the partial derivatives of functions *F*,*G* and *H* are not the same. However, we will later bound them by the same conservative quantity, so we keep this abused notation.

Calculating the partial derivatives and regrouping terms in (A 1) yields
A 3In order to prove the global asymptotic stability of synchronization, it suffices to show that the origin of system (A 3) Δ*V* =Δ*n*=Δ*S*=0 becomes globally stable when *g*_{el} exceeds some critical value. To this end, we construct a Lyapunov function in the following form:
where parameter *E*_{Ca} is always greater than *E*_{K}. The exponent 2.5 is chosen to minimize the bound on the synchronization threshold.

We need to show that the derivative of this quadratic form with respect to the trajectories of system (A 3) is negative everywhere except of the origin. Thus,
This quadratic form simplifies as follows:
where
To prove that the quadratic form is positive definite, we use the Sylvester criterion:
A 4All the conditions are satisfied if the last one is true: *A*−*B*^{2}−*C*^{2}>0. From which, it follows that
where we have used the absorbing domain bounds and ▪

The theoretical estimate for the synchronization threshold bound is conservative (*g*_{el}>3.925) compared with the numerically computed bound *g*_{el}≈0.18 shown in figure 2. However, it guarantees that the electrical coupling remains synchronizing as long as it exceeds the synchronization threshold.

We point the reader to the connection graph method [49–51], which allows us to use the bound *g** for the two-cell network to calculate the critical value of electrical coupling sufficient for globally stable synchronization in large *N*-cell networks (2.1) with arbitrary topologies of electrical connections in the absence of inhibition. The following proposition is a direct application of the connection graph method [49] to the network (2.1) with only electrical connections.

### Theorem A.2.

*Complete synchronization in the network (*2.1*) of N electrically coupled Sherman models without inhibitory connections (g*_{inh}*=0) is globally asymptotically stable if for every edge k on the connection graph associated with the connectivity matrix* ** C**
A 5

*where g**

*is the bound given in theorem A.1 for the two-cell network and the quantity*

*is the sum of the lengths of all chosen paths P*

_{ij}

*which pass through a given edge k.*

### Proof.

The proof directly follows from that of the main theorem of the connection graph method [49]. ▪

More details on the calculation of graph quantity *b*_{k} for a given network topology are given in [49,53].

## Appendix B. Numerical methods

**(a) Phase difference**

The phase difference Δ*ϕ*_{n} used for the Poincaré maps in figure 4 was introduced via a time delay between the *n*th onsets of bursting in the two cells. The time delay was normalized over the full period of bursting oscillations such that Δ*ϕ*_{n}=0 corresponds to complete synchrony and Δ*ϕ*_{n}=0.5 indicates anti-phase bursting. The phase of bursting in either cell is reset to zero when the voltage of the cell increases from its quiescent state to reach an auxiliary threshold *Θ*_{aux}=−50 (mV). The auxiliary threshold *Θ*_{aux}=−50 (mV) is chosen such that it lies between the minimum values of spikes and the quiescent state. Therefore, the time when the voltage of the reference cell increases from its quiescent state and crosses this threshold is the time for the onset of bursting. More details on the calculations of the phase differences and the Poincaré maps using this procedure can be found in [39].

**(b) Transversal Lyapunov exponent**

Transversal Lyapunov exponents for the stability of synchronization correspond to eigenvectors transversal to the synchronization subspace *M*. When all *N*−1 transversal Lyapunov exponents are negative, an initial synchronization error converges to zero, yielding stable synchronization. The largest transversal Lyapunov exponent λ_{⊥} shown in figures 3 and 8 was calculated from simulated time-series data of the variational equations (4.1) via the orbit separation algorithm [54] and the standard fourth-order explicit Runge–Kutta method of numerical integration.

## Footnotes

One contribution of 14 to a theme issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.

- Accepted March 8, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.