Suppression of Phase Synchronization in Scale-Free Neural Networks Using External Pulsed Current Protocols

Abstract

The synchronization of neurons is fundamental for the functioning of the brain since its lack or excess may be related to neurological disorders, such as autism, Parkinson’s and neuropathies such as epilepsy. In this way, the study of synchronization, as well as its suppression in coupled neurons systems, consists of an important multidisciplinary research field where there are still questions to be answered. Here, through mathematical modeling and numerical approach, we simulated a neural network composed of 5000 bursting neurons in a scale-free connection scheme where non-trivial synchronization phenomenon is observed. We proposed two different protocols to the suppression of phase synchronization, which is related to deep brain stimulation and delayed feedback control. Through an optimization process, it is possible to suppression the abnormal synchronization in the neural network.

1. Introduction

For years, the study of synchronization has been important since this phenomenon has been observed in many different biological systems, e.g., firefly communities, pacemaker cells of the heart, and crickets that chirp in unison [1,2,3,4,5]. In addition, the complexity seen in the brain is directly related to the distinct activation patterns of the neurons, which can be understood as a synchronization phenomenon. Particularly, the synchronization of neurons is important since anomalous synchronization can disrupt the brain functioning, generating disorders, such as Parkinson’s disease (PD) and autism [6,7,8,9,10].

A possible neurosurgical treatment for PD is called deep brain stimulation (DBS), which consists of the insertion of an electric probe that emits electromagnetic signals in a target brain area [9,11,12]. A more recently developed technique is the noninvasive DBS, which consists of temporally interfering electric fields [12] generated outside of the cranium. Despite its long history of use, it is still unclear how DBS works [9,10]. Some studies indicate that high-frequency DBS replaces pathological low-frequency network oscillations in the rat model of Parkinson’s disease with a regularized pattern of neuronal firing [13], and there is evidence that the DBS releases the activity patterns of groups of cells in the subthalamic nucleus that present abnormal synchronization due to PD, which destroys neurons in basal ganglia [14]. Depending on the frequency of the signal, it allows suppressing the symptoms of Parkinson’s disease [13,15].

A healthy human brain consists of $\sim {10}^{11}$ interconnected neurons through $\sim {10}^{15}$ synapses [16,17]. In the theoretical point of view, a possible way to study coupled neurons is given by the computational simulation of complex networks, where each site of the network corresponds to a neuron and its connections are represented by the edges of the network [18]. In this scenario, distinct topologies or connection architectures have been successfully used to simulate the interconnections of the human brain, such as small-world, scale-free and random topologies [18,19,20,21,22,23].

In this study, we simulated a neural network composed of $N=5000$ neurons in a scale-free topology, where each neuron was modeled by a Hodgkin–Huxley-type model [24,25,26]. This model is characterized by the insertion of two temperature sensitive parameters, and two additional slow ionic currents to the original ideas of Hodgkin and Huxley [27], which can be understood as the contribution of calcium ion channels [28].

It is observed in the literature that a neural network under a small-world topology can display abnormal phase synchronization for weak coupling regime since the phase synchronization regime in this region is characterized by a non-monotonic evolution of synchronization levels as a function of the coupling between neurons [29,30,31]. In fact, this kind of behavior has also been observed in non-identical coupled Rösller oscillators [32]. Recently, the mechanism behind the abnormal synchronization in a neural network composed of bursting neurons is explored and the relationship between the individual neuron behavior and the network synchronization helps to understand the phenomenon [33,34]. In [33], it is shown that the occurring of abnormal synchronization is related to the periodic inter-burst interval of the uncoupled neuron. Besides that, in [34], it is observed that the abnormal synchronization occurs due to an interplay between the periodic individual behavior and the influence of coupling strength, which is strong enough to induce the network to phase synchronization without destroying the influence of individual periodic behavior.

Here, we extend the study of abnormal synchronization to the scale-free connection architecture, since the topology plays an important role in the dynamics of systems [35,36,37]. Scale-free topology is different from the small-world one since the scale-free scheme presents a high degree of heterogeneity where neurons with a high connectivity degree are connected with neurons with low connectivity degree [37,38]. Thus, we study the existence of abnormal synchronization in a scale-free neural network and its suppression by the application of a disturbance in the network neurons. This perturbation is characterized by the application of a pulsed external current, which can be described as a theoretical interpretation of the DBS treatment.

We considered another suppression strategy that consists of the reapplication of a fraction of the signal generated by the neurons, which is called delayed-feedback-control (DFC). It was first applied experimentally in vitro in a spontaneously bursting neural network [39,40] and is frequently used in neural stimulation treatments [41,42].

To quantify synchronization of the network, we used the order parameter proposed by Kuramoto [43], which is able to capture information about phase synchronization of the system using data of each neuron. In this sense, we show that the suppression methods proposed are able to suppress the anomalous synchronization without affecting the regular synchronized states, which occur for large values of the coupling.

The paper is organized as follows. In Section 2, we introduce the details of the connection scheme and the used neuronal model. In Section 3, we introduce the quantification of PS by using of the Kuramoto order parameter. In Section 4, details of the perturbation methods imposed into the system, as well as the results obtained with each perturbation are discussed. Our conclusions are in the last section.

2. Neural Model and Connection Scheme

We studied the dynamical behavior of a neural network composed of 5000 neurons connected in a scale-free topology. In this case, the number of connections per node presents a statistical power-law dependence $P(n)\sim n^{-\kappa }$ [44,45]. The values of the scaling exponent are within $2\le \kappa \le 2.2$ and the average connection $\langle n\rangle \approx 4$ [21,41,46], where $P(n)dn$ is the probability to find a node with a degree in the interval from n to $n+dn$.

Scale-free networks can be obtained by the Barabasi–Albert procedure through a sequence of steps starting from an initial lattice with a small number $N_0$ of nodes randomly connected [44,45]. At each step, a new node is inserted in the network, which is randomly connected to $n\ge 2$ nodes. The process is repeated until the network reaches the desired number of nodes. In this work, we used $N=5000$ nodes. To generate a scale-free network, we used the Python library NetworkX [47], which have us $\kappa \approx 2.2$.

To simulate the individual neuron dynamics, we used a Hodgkin–Huxley-type model [25,26], where the adaptation takes into account the addition of two slow ionic fluxes. Mathematically, the neuronal model used in this work describes the temporal dynamics of the neuron membrane potential as a function of the ionic fluxes. This adaptation also includes temperature sensitive parameters. The temporal evolution of the membrane potential $V_i$ is described by

$$\mathcal{C}\frac{d{V}_i}{dt}=-J_{i,\mathrm{Na}}-J_{i,\mathrm{K}}-J_{i,\mathrm{sd}}-J_{i,\mathrm{sa}}-J_{i,\mathrm{L}}+J_{i,\mathrm{coup}},$$

(1)

where $\mathcal{C}$ is the specific membrane capacitance of neurons measured in µF/cm²; $V_i$ is measured in mV; $J_{i,\mathrm{Na}}$, $J_{i,\mathrm{K}}$, and $J_{i,\mathrm{L}}$ are the sodium, potassium and non-gated channels fluxes of the original Hodgkin–Huxley model [27], respectively, which are measured in µA/cm²; and $J_{i,\mathrm{sd}}$ and $J_{i,\mathrm{sa}}$ are the two slow ionic fluxes added by Braun et al. to this model and the are associate to calcium flux [28].

The electrical fluxes related to the ion and leak channels are given by conductance-based expressions [25]

$$\begin{array}{ccc}\hfill J_{i,\mathrm{Na}}& =& \rho {\overline{g}}_{\mathrm{Na}}{\alpha }_{i,\mathrm{Na}}(V_i-E_{\mathrm{Na}}),\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill J_{i,\mathrm{K}}& =& \rho {\overline{g}}_{\mathrm{K}}{\alpha }_{i,\mathrm{K}}(V_i-E_{\mathrm{K}}),\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill J_{i,\mathrm{sd}}& =& \rho {\overline{g}}_{\mathrm{sd}}{\alpha }_{i,\mathrm{sd}}(V_i-E_{\mathrm{sd}}),\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill J_{i,\mathrm{sa}}& =& \rho {\overline{g}}_{\mathrm{sa}}{\alpha }_{i,\mathrm{sa}}(V_i-E_{\mathrm{sa}}),\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill J_{i,\mathrm{L}}& =& {\overline{g}}_{\mathrm{L}}(V_i-E_{\mathrm{L}}),\hfill \end{array}$$

(6)

where ${\overline{g}}_{\mathrm{Na}}$, ${\overline{g}}_{\mathrm{K}}$, ${\overline{g}}_{\mathrm{sd}}$, ${\overline{g}}_{\mathrm{sa}}$, and ${\overline{g}}_{\mathrm{L}}$ are the maximum (specific) conductances measured in mS/cm², and $E_{\mathrm{Na}}$, $E_{\mathrm{K}}$, $E_{\mathrm{sd}}$, $E_{\mathrm{sa}}$, and $E_{\mathrm{L}}$ denote the reversal Nernst potentials for each ionic current measured in mV. The term $\rho$ refers to a temperature dependence of the model and it is described by

$$\rho ={\rho }_0^{(T-T_0)/{\tau }_0},$$

(7)

where T, $T_0$ and ${\tau }_0$ and ${\rho }_0$ are constants of the model.

The temporal evolution of the activation functions ${\alpha }_{i,\mathrm{Na}}$, ${\alpha }_{i,\mathrm{K}}$,${\alpha }_{i,\mathrm{sd}}$, and ${\alpha }_{i,\mathrm{sa}}$ are described by

$$\begin{array}{ccc}\hfill \frac{d{\alpha }_{i,\mathrm{Na}}}{dt}& =& \frac{\varphi }{{\tau }_{\mathrm{Na}}}({\alpha }_{i,\mathrm{Na},\infty }-{\alpha }_{i,\mathrm{Na}}),\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill \frac{d{\alpha }_{i,\mathrm{K}}}{dt}& =& \frac{\varphi }{{\tau }_{\mathrm{K}}}({\alpha }_{i,\mathrm{K},\infty }-{\alpha }_{i,\mathrm{K}}),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill \frac{d{\alpha }_{i,\mathrm{sd}}}{dt}& =& \frac{\varphi }{{\tau }_{\mathrm{sd}}}({\alpha }_{i,\mathrm{sd},\infty }-{\alpha }_{i,\mathrm{sd}}),\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \frac{d{\alpha }_{i,\mathrm{sa}}}{dt}& =& \frac{\varphi }{{\tau }_{\mathrm{sa}}}(-\eta J_{i,\mathrm{sd}}-\gamma {\alpha }_{\mathrm{i},\mathrm{sa}}),\hfill \end{array}$$

(11)

where ${\tau }_{\mathrm{Na}}$, ${\tau }_{\mathrm{K}}$, ${\tau }_{\mathrm{sd}}$, and ${\tau }_{\mathrm{sa}}$ are constants [25]. The parameter $\eta$ works to increase calcium ion concentration following $J_{i,\mathrm{sa}}$, while $\gamma$ accounts for active the elimination of intracellular calcium [28]. $\varphi$ is another temperature dependence of the model, given by

$$\varphi ={\varphi }_0^{(T-T_0)/{\tau }_0}.$$

(12)

The functions ${\alpha }_{i,\mathrm{N}a,\infty }$, ${\alpha }_{i,\mathrm{K},\infty }$, and ${\alpha }_{i,\mathrm{s}d,\infty }$ are described by

$${\alpha }_{i,\mathrm{Na},\infty }=\frac{1}{1+exp[-s_{\mathrm{Na}}(V_i-V_{0\mathrm{Na}})]},$$

(13)

$${\alpha }_{i,\mathrm{K},\infty }=\frac{1}{1+exp[-s_{\mathrm{K}}(V_i-V_{0\mathrm{K}})]},$$

(14)

$${\alpha }_{i,\mathrm{sd},\infty }=\frac{1}{1+exp[-s_{\mathrm{sd}}(V_i-V_{0\mathrm{sd}})]},$$

(15)

where $s_{\mathrm{Na}}$, $s_{\mathrm{K}}$, $s_{\mathrm{sd}}$, $V_{0\mathrm{Na}}$, $V_{0\mathrm{K}}$, and $V_{0\mathrm{sd}}$ are constants whose values are given in

Table 1. Parameter values of the neuronal model according to Braun et al. [25].

Membrane Capacitance	$\mathcal{C}$ = 1.0 µF/cm2
Characteristic Times (ms)	${\tau }_{\mathrm{Na}}=0.05$	${\tau }_{\mathrm{K}}=2.0$	${\tau }_{\mathrm{sd}}=10$	${\tau }_{\mathrm{sa}}=20$
Maximum Conductances (mS/cm2)	${\overline{g}}_{\mathrm{Na}}=1.5$	${\overline{g}}_{\mathrm{K}}=2.0$	${\overline{g}}_{\mathrm{sd}}=0.25$	${\overline{g}}_{\mathrm{sa}}=0.4$
	${\overline{g}}_{\mathrm{L}}=0.1$
Reversal Potentials (mV)	$E_{\mathrm{Na}}=50$	$E_{\mathrm{K}}=-90$	$E_{\mathrm{sd}}=50$	$E_{\mathrm{sa}}=-90$
	$E_{\mathrm{L}}=-60$	$V_{0\mathrm{Na}}=-25$	$V_{0\mathrm{K}}=-25$	$V_{0\mathrm{sd}}=-40$
Other Parameters	${\rho }_0=1.3$	${\varphi }_0=3.0$	$T_0=25$ °C	${\tau }_0=10$ °C
	$s_{\mathrm{Na}}=0.25$ (1/mV)	$\eta =0.012$ µA	$\gamma =0.17$	$s_{\mathrm{K}}=0.25$ (1/mV)
	$s_{\mathrm{sd}}=0.09$ (1/mV)	$T=13$ °C

following Braun et al. [25].

The coupling term $J_{i,\mathrm{coup}}$, in Equation (1), is an excitatory chemical synapse, since the synapse does not occur directly. In this way, the ith postsynaptic neuron receives signals from presynaptic ones [48]

$$J_{i,\mathrm{coup}}=\frac{\epsilon }{\langle n\rangle }\sum _{j=1}^N{e}_{i,j}{r}_j(V_{\mathrm{syn}}-V_i),$$

(16)

where $\epsilon$ is the coupling parameter that controls the coupling intensity. $\langle n\rangle$ is the normalization factor, defined as the average of connections number, which is $\langle n\rangle \approx 4$. $V_{syn}$ is the synaptic reversal potential, set here as 20 mV, which assures that the contribution coming from the coupling is positive for all instant of time, characterizing an excitatory synapse. $e_{i,j}$ represents the elements of the adjacency matrix, which is a scale-free type. In this case, if the ith and jth neurons are connected, $e_{i,j}=1$; otherwise, $e_{i,j}=0$.

Added to the kinetic variable of the model, $r_i$ refers to the fraction of bound receptors available to receive a connection. We used the equation of $r_i$ proposed by Destexhe et al. [48],

$$\frac{d{r}_i}{dt}=\left(\frac{1}{{\tau }_r}-\frac{1}{{\tau }_d}\right)\frac{1-r_i}{1+exp[-s_0(V_i-V_0)]}-\frac{r_i}{{\tau }_d},$$

(17)

where $s_0$ is a unitary constant measured in (1/mV), $V_0=-20$ mV, and ${\tau }_r=0.5$ ms and ${\tau }_d=8$ ms are constants associated to the rises and decays of the synaptic transition, respectively.

To integrate the set of coupled equations composing the model, we used Adams’ predictor-corrector method [49] with an absolute tolerance less than ${10}^{-8}$. Figure 1a depicts the typical membrane potential for a neuron, using the fixed set of parameter values shown in Table 1. As observed, the neuron depicts bursting dynamics characterized by a sequence of spikes followed by a resting time [50]. We refer to this dynamics as bursting regime and, throughout the coupling interval used here, the suppression process ensured that this regime is not lost, which makes possible the phase association and synchronization evaluation for all interval of coupling and suppression strength.

3. Phase Synchronization Quantifier

To quantify phase synchronization in the bursting regime, we associated a geometric phase to the sequence of bursts for each neuron. Figure 1b shows the auxiliary variable $U_i\equiv 1/{\alpha }_{i,\mathrm{sa}}$ computed using Equation (11), where each maximum of $U_i$ corresponds to the beginning of a burst of the ith neuron [51]. If $t_{k,i}$ is the beginning time of the kth burst of the ith neuron, the duration of the burst would be $t_{k+1,i}-t_{k,i}$, with $k=0,1,2\cdots$ and $i=1,2,\cdots ,N$, consequently the phase would vary from $2\pi k$ to $2\pi (k+1)$, and it is defined for specific time t as [52]

$${\theta }_i(t)=2\pi k_i+2\pi \frac{t-t_{k,i}}{t_{k+1,i}-t_{k,i}},t_{k,i}\le t<t_{k+1,i}.$$

(18)

Considering the geometric phase variable ${\theta }_i$ as defined in Equation (18), to quantify PS, we used the modulus of the Kuramoto order parameter $R(t)$ [53]

$$R(t)=\left|\frac{1}{N}\sum _{i=1}^N{e}^{i{\theta }_i(t)}\right|,$$

(19)

where $R(t)$ gives us a number between 0 (completely unsynchronized) and 1 (completely phase synchronized).

The order parameter oscillates in time for not fully synchronized neurons [43] and its temporal mean value is defined as

$$\langle R\rangle =\frac{1}{M}\sum _{j=1}^{M}R(t_j^{\prime }),$$

(20)

being $t_1^{\prime }=t_i,t_2^{\prime }=t_i+h,\cdots ,t_M^{\prime }=t_f$, where $h=0.01$, and $t_i$ and $t_f$ are initial and final times of the computation of $R(t)$.

To show the synchronization behavior of the network, Figure 2 depicts $\langle R\rangle$ as a function of the coupling parameter $\epsilon$ of a neural network given by Equations (1)–(17). To avoid any trend in the result, we used random initial conditions in the following intervals: $V_i\in [-65.0;0.0]$; ${\alpha }_{\mathrm{Na},i},{\alpha }_{\mathrm{K},i},{\alpha }_{\mathrm{sd},i},{\alpha }_{\mathrm{sa},i},r_i\in [0.1;1.0]$. Observe that, for $\epsilon >{\epsilon }^{*}=0.007$ (mS/cm²), the network followed a route of globally stable phase synchronized state, since $\langle R\rangle$ approached 1 as the coupling strength was increased. For the interval of coupling strength $0.002<\epsilon <0.007$, the network exhibited a local maximum of phase synchronization ($\epsilon \approx 0.004$). This behavior was also observed in small-world network [30,31,51], which characterized a non-monotonic evolution of the synchronization level as a function of the coupling that could be understood as an abnormal synchronization since PS occurred for a coupling $\epsilon <{\epsilon }^{*}$. In this way, it is known that several brain disorders, such as Parkinson’s disease and autism, are related to abnormal neuronal synchronization [6,7,8,9], thus it is expected that the application of synchronization suppression methods may be useful to vanish the anomalous synchronization, as observed in Figure 2.

4. Results and Discussions

Considering a scale-free neural network composed of $N=5000$ identical neurons, we used the mean value of the Kuramoto order parameter $\langle R\rangle$ to evaluate the PS for different suppression synchronization protocols. Here, we used a transient time given by $t_i=150$ s and the total simulation time was set to $t_f=250$ s.

Motivated by experimental results [12,15], we made a perturbation in the network by applying an external pulsed current $\lambda (t)$ in Equation (1) of the neuronal model. Mathematically, the suppression method can be described as

$$\lambda (t)=\frac{{\lambda }_0}{2}+\sum _{m=1,3,\cdots }^{\infty }\frac{2{\lambda }_0}{m\pi }sin\left(\frac{2m\pi t}{\tau }\right),$$

(21)

where ${\lambda }_0$ is the amplitude of each pulse measured in µA/cm², and $\tau$ is the period for which the current is successively turned on and off and measured in seconds. Figure 3 shows the evolution of $\lambda (t)/{\lambda }_0$ as a function of $t/\tau$.

Figure 4 depicts $\langle R\rangle$ as a function of the coupling $\epsilon$ for different values of amplitude ${\lambda }_0$. For amplitudes lower than ${\lambda }_0<0.05$ (µA/cm²), the network still presented PS for small values of $\epsilon$. However, for ${\lambda }_0\ge 0.05$, the anomalous PS was suppressed without altering the globally stable state of synchronization for coupling value higher than the critical value ${\epsilon }^{*}$ (which remains constant in ${\epsilon }^{*}\approx 0.007$ for ${\lambda }_0<0.2$). The frequency $\nu =1/\tau$ was fixed at 140 Hz (which means $\tau \approx 0.0071$ s), since experimental results show that only a high frequency currents ($\nu >100$ Hz [9,10,54]) could restore normal neural behavior in Parkinson’s Disease (PD) [13].

The next step consisted of the study of the heterogeneity of the scale-free network because one of the characteristics of the scale-free topology is the existence of hubs, which are characterized by neurons with high connectivity in the network [41,44,45]. Intuitively, it is believed that the perturbation presents greater influence when applied in the hubs since they have high connectivity in the network. We made a change in the applied current to apply the current in select groups of neurons

$$\lambda (t,{\lambda }_0,\tau )\to {\lambda }_i(t,{\lambda }_{0,i},\tau )$$

where ${\lambda }_{0,i}={\lambda }_0$ if $i\subset \mathcal{G}$ or ${\lambda }_{0,i}=0$, otherwise $\mathcal{G}$ is a subset of neurons in the network. Here, it was studied how the PS varied for three different subsets $\mathcal{G}$. Firstly, we applied the current in the neurons with higher connectivity degree of the network, in this case, the order of $\mathcal{G}$ was $\left|\mathcal{G}\right|=N_{\mathrm{hubs}}$. In the second case, the pulsed current was applied in random neurons of the network, and then $\left|\mathcal{G}\right|=N_{\mathrm{rand}}$. In the latter case, a neuron in the network was randomly chosen and the current was applied in that neuron and its neighbors, which formed a package of neurons that received the application with $\left|\mathcal{G}\right|=N_{\mathrm{package}}$. Figure 5 depicts an example of a subset $\mathcal{G}$ with ${\lambda }_{0,i}/{\lambda }_0=1$, in this case $\left|\mathcal{G}\right|=1000$. The Figure 5a shows a subset of random neurons and Figure 5b a subset of a package of neurons.

Figure 6 depicts the mean value of Kuramoto order parameter $\langle R\rangle \times \epsilon$ as a function of the order of the subset $\left|\mathcal{G}\right|$ with ${\lambda }_0=0.1$. In Figure 6a a subset of the neurons with higher connectivity in the network with $\left|\mathcal{G}\right|=N_{\mathrm{hubs}}$ is chosen, Figure 6b a subset of random neurons with $\left|\mathcal{G}\right|=N_{\mathrm{rand}}$ and Figure 6b a package of neurons with $\left|\mathcal{G}\right|=N_{\mathrm{package}}$. Note that the three surfaces have the same shape. In this case, when the order of $\left|G\right|\gtrsim 2000$, the PS was suppressed, that is, now the network depicted a monotonic evolution of the synchronization as $\epsilon$ increased.

Another strategy consisted in the application of a delayed mean field signal $\overline{V}$ over the network

$$\xi (t)={\xi }_0\overline{V}(t-t_{\mathrm{delay}}),$$

(22)

where ${\xi }_0$ (similarly to ${\lambda }_0$) is the current amplitude given by ${10}^{-4}$ µA/cm², and $t_{\mathrm{delay}}$ (ms) is the delay time between the generation of the mean field and the re-application of the signal. $\overline{V}$ is the mean field potential of the network, which is defined by

$$\overline{V}(t)=\frac{1}{N}\sum _{i=1}^N{V}_i(t).$$

(23)

Figure 7 depicts the $\overline{V}(t)$ for three different values of coupling $\epsilon$. In Figure 7a, $\epsilon =0.001$, the mean field display a small amplitude variation since the network presents an unsynchronized behavior. In Figure 7b, $\epsilon =0.007$ and, in Figure 7c, $\epsilon =0.020$, the mean field display an oscillatory behavior, since in this regime the neurons of the network presents a signal of partial synchronization.

In this approach, $\overline{V}<0$, ${\xi }_0>0$ (${\xi }_0<0$) characterizes an inhibitory (excitatory) current which decreases (increases) the membrane potential $V_i$. The natural period of the Hodgkin–Huxley-type neuron with the parameters in Table 1 is $t_0\approx 1,250$ ms. In Figure 8, we show how the PS varies with the application of $\xi (t)$ in all the neurons as a function of the amplitude ${\xi }_0$ and the coupling parameter $\epsilon$ for three different delay time $t_{\mathrm{delay}}$, in panel Figure 8a, $t_{\mathrm{delay}}=0$, that is, the mean field affects the network instantaneously; when $-15<{\xi }_0<-5$, the network is characterized by a monotonic evolution of the synchronization the anomalous PS is suppressed; when ${\xi }_0>-5$, the anomalous regime still occurs; and, particularly, when ${\xi }_0>5$ the synchronization for coupling $\epsilon <{\epsilon }^{*}$ is amplified with $\langle R\rangle \approx 0.95$. In Figure 8b, $t_{\mathrm{delay}}=500$, which is $t_{\mathrm{delay}}\approx t_0/2$, and the mean field with delay $\overline{V}(t-t_0/2)$ is in anti-phase with the $\overline{V}(t)$, for ${\xi }_0<-7.5$ the anomalous PS is suppressed, otherwise the network still depicts a non-monotonic evolution of the $\langle R\rangle$ which characterizes abnormal synchronization. In Figure 8c, $t_{\mathrm{delay}}=1000\approx t_0$, as expected; the result is similar to Figure 8a because of the oscillatory behavior of the mean field, that is $\overline{V}(t-t_0)\approx \overline{V}(t)$, the numerical value of $\xi (t)$ is the same in both cases.

5. Conclusions

In this paper, we model a neural network composed of $N=5000$ Hodgkin–Huxley-type neurons to study the synchronization phenomena. A similar approach have been used to analyze small-world neural networks [31,33,55]. However, the influence of topology plays an important role in the synchronization characteristics [22,35,38]. In this way, we simulated a scale-free network since there are great differences regarding the heterogeneity of the network in comparison to the small-world one [37,38]. It was shown that the scale-free network displays a non-monotonic evolution of the phase synchronization as the coupling between neurons increases. A similar scenario has been observed in small-world networks, which is called “anomalous phase synchronization” [30,31,55], since the traditional behavior should monotonically transition to PS [33]. Especially, Parkinson’s disease and some episodes of seizure behavior generated by epilepsy may be associated to anomalous synchronization.

We have proposed two methods of suppression of the anomalous synchronization behavior, both based on treatment for neurological disorders, which consist in the application of an external current in the neurons of the network [9,10].

The first method consists of electrical pulses imposed all over the network. It was shown that, for an amplitude higher than a critical value ${\lambda }_0>0.05$, the anomalous synchronization was suppressed. As a second approach, we studied how the heterogeneity of the scale-free network affects the anomalous PS. We used three different protocols applying the pulsed current in subsets of hubs,; random neurons, and a selected package of neurons. We showed the existence of a threshold, 2000 neurons (40% of the network), which must be disturbed to reach the suppression of anomalous PS in all cases. Such a conclusion implies that the synchronization is related to the individual dynamics of each neuron rather than the network topology [33].

In the second method, only a small fraction ${\xi }_0>-0.0005$ µA/cm² (with $t_{\mathrm{delay}}=0$) of the delayed signal of the mean field was applied to all neurons and the abnormal synchronization was suppressed.

Finally, we showed that the delayed signal of the mean field potential had a greater contribution in the region where the suppression was not reached. When $t_{\mathrm{delay}}\approx t_0/2$, it was observed that the anomalous synchronization still persisted but with a lower intensity ($\langle R\rangle \approx 0.6$) compared to the non-delayed scenario ($t_{\mathrm{delay}}=0$).