Inhibition of spikes in an array of coupled FitzHugh – Nagumo oscillators by external periodic forcing

Elena Adomaitienė, Gytis Mykolaitis, Skaidra Bumelienė, Arūnas Tamaševičius Department of Electronics, Center for Physical Sciences and Technology, Saulėtekio ave. 3, LT-10223 Vilnius, Lithuania elena.tamaseviciute@ftmc.lt; skaidra.bumeliene@ftmc.lt; arunas.tamasevicius@ftmc.lt Department of Physics, Vilnius Gediminas Technical University, Saulėtekio ave. 11, LT-10223 Vilnius, Lithuania gytis.mykolaitis@vgtu.lt


Introduction
Undesirable instabilities in dynamical systems can be avoided by applying conventional proportional feedback techniques [8,11].An example is a simple second-order system, where the proportional feedback is given by a linear term with a control coefficient k: Here F and G are either linear or nonlinear functions, the x * is a reference point, e.g. a steady state coordinate of the system.However, in many real systems, especially in biology, chemistry, physiology, etc., the exact locations of these states are unknown.Moreover, their positions may vary with time because of unknown and unpredictable forces.Therefore, adaptive methods, automatically tracing and stabilizing the steady states, are required.A large number of adaptive control techniques have been developed so far, e.g. the tracking filter method [9,10,15], and applied to a variety of dynamical systems (see [1,19] and references therein).To implement the tracking filter technique, system (1) c Vilnius University, 2017 should be provided with an additional equation describing the dynamical variable z of the first-order filter: where ω f is the cut-off frequency of the filter (usually ω f 1).An alternative control method is a non-feedback technique based on applying to the system external periodic force: In system (3), the frequency of the external forcing ω should be high enough in comparison with the natural frequency of the uncontrolled dynamical system.A specific example is the stabilization of the unstable upside-down position of a mechanical pendulum by vibrating its pivot up and down at a relatively high frequency [20].Recently, this "mechanical" idea has been exploited in a seemingly unexpected field [12], namely, to get insight into the mechanism of the so-called deep brain stimulation (DBS) conventionally used for patients with the Parkinson's disease, essential tremor [2][3][4], and other brain malfunctions.
In this paper, we extend this research by demonstrating that external periodic forcing can inhibit spikes in an array of coupled neuronal oscillators.To be specific, we consider an array of the mean-field coupled electronic FitzHugh-Nagumo (FHN) oscillators, also known in literature as the Bonhoeffer-van der Pol oscillators [14].

Analogue circuits
The corresponding circuit diagrams are presented in Fig. 1.CN is a coupling node.It is assumed that the CN is not accessible directly from the outside, but via some passive resistance network represented here for simplicity by an equivalent resistance R g .DN is an accessible damping node.
In ).It is a slight modification of an oscillator described in [16] and essentially differs from the earlier asymmetric version of the FHN-type oscillator suggested in [5].In the experiments, we employed a hardware array with N = 30 described in details (without any external control) elsewhere [17].The external inhibitory current I inh (t) = I A sin(2πf t) was injected from an external sine wave generator via the damping node DN.

Electronic experiments
In real neuronal systems, the typical spike parameters are the following: spike height is about 100 mV, spike width is about 1 ms, interspike interval ranges from several hundreds to several tens of milliseconds [7].Consequently, the repetition rate f 0 of the spikes is from several Hz to several tens of Hz.Correspondingly, the practical DBS frequencies are chosen between 50 and 300 Hz [4].Though an analog electronic circuit exhibiting the above neuron parameters, e.g. with f 0 =10 Hz, can be implemented as a hardware [18], in the present modelling experiment, we employed an array of FHN-type oscillators described in [17].In this array, the frequencies of the individual units were intentionally set at 12 kHz, i.e. by about three orders higher than in real neurons.Such an increase of working frequency makes the analog experiment more convenient.All experimental data, in particular phase portraits, Poincaré sections, and power spectra can be taken on a real time scale [17].Also, time averaging procedure of the signals become extremely fast.
For the best performance, it is necessary to choose an appropriate drive frequency f and amplitude I A .The f should be much higher than the natural frequency f 0 of the spiking oscillators.The threshold value of the drive amplitude I * A depends on the frequency f .The I * A is relatively low in the interval from 50 to 500 kHz.Outside this range (below 50 and above 500 kHz), I * A increases rapidly.Similar amplitude-frequency dependence (a valley in a certain frequency range) was observed in real DBS experiments [3] with f from 100 until 5000 Hz, also for single electronic neuron model in low frequency experiment [18], where the optimal range of f was from 40 to about 400 Hz.In our present high frequency electronic experiments, we chose f = 150 kHz, i.e. in the middle between 50 and 500 kHz.The selected frequency provided the lowest threshold I * A = 50 mA.The experimental results are shown in Figs. 2 and 3 by the waveforms and the phase portraits, respectively.Here the V C is the mean-field voltage of the voltages V Ci from the individual oscillators (i = 1, 2, . . ., 30).
The time average of the high frequency non-spiking voltage V C (right-hand side of the bottom plot), taken over the period (T = 1/f ) of the external current, is ŪC ≈ −0.18 V.It is non-zero value because of the DC bias V 0 .The ŪC is noticeably different from the natural steady state V 0C = −0.27V measured in a non-oscillatory mode (when the all coils L are short-circuited).The self-sustained low frequency (f 0 ≈ 12 kHz) spikes of about 3 V height are totally suppressed when the inhibitory current I A I * A = 50 mA is injected.However, we have a finite (≈ 13%) higher frequency artefact.The voltage oscillates around the time average ŪC with the amplitude of about 0.4 V at the external drive frequency f .Moreover, the artefact voltage continues to change (Fig. 4) when the external drive amplitude I A is increased above the threshold value I * A (the amplitude I A should be higher than the threshold to guarantee robust inhibition).For example, at a double drive amplitude, I A /I * A = 2 the average voltage changes its sign.Similar behaviour was observed earlier, but not emphasized in the numerically simulated bifurcation diagram for the Hodgkin-Huxley (HH) single neuron model [12].

Mathematical model
Applying the Kirchhoff laws to the circuits in Fig. 1 with R 1 = R 2 and R 7 max{ L/C, R 3 , R 4 , R 5 , R 6 }, the following differential equations are derived: The nonlinear current-voltage (I-V ) characteristic 4) is approximated by three segments of linear functions Here V * is the breakpoint voltage of the forward I-V characteristic of the diodes (V * ≈ 0.6 V).In system (4), the individual oscillators are coupled via the mean-field voltage We introduce the following set of dimensionless variables and parameters: also, two additional dimensionless parameters for the external sine wave forcing: and we arrive to a set of 2N coupled non-autonomous differential equations convenient for numerical integration: The f (x i ) in system ( 8) is a nonlinear function presented by a piecewise linear function Note that, due to d 1 d 2 , the f (x i ) is an essentially asymmetric function [16] in contrast to the common FHN cubic parabola x 3 introduced by FitzHugh [6].The DC bias parameters c i are intentionally set different for each individual oscillator, thus making them non-identical units.

Numerical results
Integration of system (8) has been performed using the Wolfram Mathematica package.The numerical results are presented in Fig. 5.They are in a good agreement with the experimental plots in Fig. 2. The mean-field variable x does not converge to a constant steady state, but oscillates around it at the drive frequency.Strictly speaking, the nonautonomous (externally driven) dynamical systems, e.g.given by system (8), do not possess steady states at all.Only in the case of high frequency (f f 0 ) drive, we can introduce the average values taken over the external period.These averages are related to the real steady states.6 Mean-field approach Analysis of system ( 8) can be essentially simplified if we consider only the mean-field variables obtained by the direct averaging of x i , y i , f (x i ), k( x − x i ), and c i in the original equations: As a result, the mean of the coupling term k( x −x i ) = k( x − x ) = 0 in system (9) has been nullified independently on the value of k.Further, we assume that all |x i | 1.
According to definition of the nonlinear function f (x i ), this leads to f Eventually, we obtain a set of linear differential equations, which do not describe the full dynamics of the mean field, but provide its steady state.In the absence of the external drive (A = 0), it has the following coordinates (for ab < 1 and Stability analysis of system (9) shows that, for A = 0 and a > b, the steady state is unstable (the real parts of the both eigenvalues of the corresponding second-order characteristic equation are both positive).If (in addition to a > b) the sum a + b > 2, then the eigenvalues are real (no imaginary parts).Thus, the steady state is an unstable node.Whereas the external periodic forcing (A = 0), similarly to the mechanical pendulum [20] and the single HH neuron [12], stabilizes the originally unstable steady state.
For the parameter values employed in numerical simulations: a = 3.4, b = 0.16, and c i = −44/(24 + i), the steady-state coordinates have the following numerical values: x 0 = −0.41,y 0 = −2.57.Using the definitions of the dimensionless variables introduced in (6), we estimate the means of the steady-state coordinates of the original system: V 0C ≈ −0.25 V, I 0L ≈ −1 mA.The estimated steady-state voltage V 0C is close to its experimental value −0.27 V.

Conclusions
The recent papers [12,13] on suppression of neuronal spikes by means of periodic force consider mathematical models of single neurons.Our paper [16] on an adaptive feedback technique for damping neuronal activity also deals with a single oscillator only.The present work extends the above investigations to an array of coupled (synchronized) neuronal oscillators.In addition to numerical simulations, we have carried out a hardware experiment using an analog electronic network described in [17].This research can serve for better understanding the mechanism of the DBS technique.
The influence of strongly perturbed steady states of the neurons on the effectiveness of DBS technique has not been investigated yet.However, one can suppose that the high frequency artefact oscillations (Fig. 2) and especially its unnatural DC component (Fig. 4) probably can cause the undesirable side effects in the real neuronal systems.

Figure 2 .Figure 3 .
Figure 2. Experimental waveforms of the external periodic current I inh and the mean-field voltage of the array V C ; f =150 kHz.

Figure 4 .
Figure 4. Time average of the mean-field voltage ŪC , taken over the period (T = 1/f ) of the external inhibitory current I inh , as a function of the normalized amplitude I A /I * A of the external current.I * A = 50 mA.Extrapolation to zero control (I A = 0) provides a value of ŪC close to the natural steady state V 0C = −0.27V (dashed line in the plot).