Analysis of Large-Amplitude Pulses in Short Time Intervals: Application to Neuron Interactions

Abstract

This paper deals with the analysis of a nonlinear dynamical system which characterizes the axons interaction and is based on a generalization of FitzHugh-Nagumo system. The parametric domain of stability is investigated for both the linear and third-order approximation. A further generalization is studied in presence of high-amplitude (time-dependent) pulse. The corresponding numerical solution for some given values of parameters are analyzed through the wavelet coefficients, showing both the sensitivity to local jumps and some unexpected inertia of neuron's as response to the high-amplitude spike.

1. Introduction

Theclassical model of one excitatory neuron [1–11] is considered under a spike train (time-dependent) excitation. The dynamical system as well as its evolution are investigated under the dependence on some featuring parameters and in presence of a high-amplitude time-dependent pulse. It was already shown, in a previous paper [12], dealing with Hodgkin-Huxley model, that neuron's multiple firing does not reflects immediately on a direct response and there is some kind of time delay (inertia) before the pulse becomes effective. It should be noticed that there are different models of neurons (see, e.g., [1] and references therein) such as the Integrate-and-fire, FitzHugh-Nagumo [2, 9], Morris-Lecar, and the more general Hodgkin-Huxley model [6–8]. In the early sixties, FitzHugh [2] proposed, as a simplified model, a generalization of the Van der Pol equation thus showing the existence of a limit cycle and some periodicity. In all models, if we consider a simple system consisting of a neuron and a synapse, it is known that the activity of stimulated axons can be detected by an abrupt change in the electrical potential. These pulse in a short time are called spikes or axons firing. However, as shown later these abrupt changes do not appear immediately, thus showing, even in this simplified model of neurons interactions, the existence of some kind of neurons inertia, which can be physically explained by their cooperative responce to firing.

In the following we will propose a suitable generalization of the FitzHugh-Nagumo model (FHN) and we will study both on the linear and third-order approximation. Equilibrium points, parametric domain of stability will be outlined in absence of time-dependence. The presence of high-amplitude time-dependent source of pulses can be studied on the nonlinear time series of the numerical approximation. The pulse action of the firing process, based on FitzHugh suggestion of delta Dirac function, is proposed in the form of a high-amplitude localized function with compact support in a short interval, which replicates a few times before it disappears. According to FitzHugh any change in the electrical potential are localized in a short time interval and within this interval has a significant amplitude.

2. Simplified FitzHugh-Nagumo Model

The most successful and general model in neuroscience describes the neuron activity in terms of two conductances and electric potentials. The Hodgkin-Huxley model [6–8] considers the neuron activity as an electrical circuit. Cells membrane store charges, electrochemical forces arise because of the imbalanced ion concentration inside and outside the cell. It can be written, in adimensional form, as 𝑑𝑢𝑑𝑡=𝑎₁𝑛⁴(𝑎₂−𝑢)+𝑎₃𝑚³ℎ(𝑎₄−𝑢)+𝑎₅(𝑎₆−𝑢)−𝐼_ext(𝑡),𝑑𝑛𝑑𝑡=𝛼_𝑛(𝑢)(1−𝑛)−𝛽_𝑛(𝑢)𝑛,𝑑𝑚𝑑𝑡=𝛼_𝑚(𝑢)(1−𝑚)−𝛽_𝑚(𝑢)𝑚,𝑑ℎ𝑑𝑡=𝛼_ℎ(𝑢)(1−ℎ)−𝛽_ℎ(𝑢)ℎ, (2.1) being 𝛼_𝑛𝑎(𝑢)=𝜙⋅⁷(𝑢−𝑎₈)1−𝑒^{−(𝑢−𝑎8)/𝑎9},𝛽_𝑛(𝑢)=𝜙⋅𝑎₁₀𝑒^{−(𝑢−𝑎11)/𝑎12},𝛼_𝑚𝑎(𝑢)=𝜙⋅¹³(𝑢−𝑎₁₄)1−𝑒^{−(𝑢−𝑎14)/𝑎15},𝛽_𝑚(𝑢)=𝜙⋅𝑎₁₆𝑒^{−(𝑢−𝑎17)/𝑎18},𝛼_ℎ(𝑢)=𝜙⋅𝑎₁₉𝑒^{−(𝑢−𝑎20)/𝑎21},𝛽_ℎ𝑎(𝑢)=𝜙⋅²²1+𝑒^{−(𝑢−𝑎23)/𝑎24}. (2.2)

An alternative simplified 2-dimensional model (FHN) that has been proposed by FitzHugh-Nagumo [2, 9] can be substantially represented by the system𝑑𝑥𝑥𝑑𝑡=𝑎(𝑥−³3)−𝑏𝑦+𝛾𝑧,𝑑𝑦𝑑𝑡=𝛼𝑥−𝛽𝑦, (2.3) where 𝑧=𝑧(𝑡) is a pulse function defined in a very short-range interval (pulse function), having as limiting case the delta Dirac. In the following we will consider a spike train. Parameter 𝛾 defines the amplitude of the pulse function 𝑧(𝑡).

This dynamical system depends [2, 9] on 𝑎,𝑏,𝛼,𝛽,and𝛾 in the sense that the evolution would be completely different, starting from some critical values. In the following we will show that the solution tends asymptotically to a constant value.

A suitable generalization of system (2.3) is𝑑𝑥𝑑𝑡=𝑎sin𝑥−𝑏ℎ(𝑦)+𝛾𝑧(𝑡),𝑑𝑦𝑑𝑡=𝛼𝑓(𝑥)−𝛽𝑔(𝑦), (2.4) with nonnegative parameters 𝑎≥0,𝑏≥0,𝛽≥0, (2.5) and odd functions 𝑓(𝑥)=−𝑓(−𝑥),ℎ(𝑦)=ℎ(−𝑦),𝑔(𝑦)=𝑔(−𝑦),𝑓(𝑥₀)=𝑓^′′(𝑥₀)=𝑓^′′′(𝑥₀)=0,ℎ(𝑥₀)=ℎ^′′(𝑥₀)=ℎ^′′′(𝑥₀)=0,𝑔(𝑥₀)=𝑔^′′(𝑥₀)=𝑔^′′′(𝑥₀)=0. (2.6) Up to the first-order it is𝑑𝑥𝑑𝑡=𝑎𝑥−𝑏𝑦+𝛾𝑧(𝑡),𝑑𝑦𝑑𝑡=𝛼𝑥−𝛽𝑦, (2.7) and to the third-order it is𝑑𝑥𝑥𝑑𝑡=𝑎(𝑥−³6)−𝑏𝑦+𝛾𝑧(𝑡),𝑑𝑦𝑑𝑡=𝛼𝑥−𝛽𝑦, (2.8) so that the FHN system (2.3) might be considered as the third-order approximation of system (2.4).

2.1. Autonomous System 𝛾=0

2.1.1. Linear Case

When 𝛾=0, and assuming by analogy with system (2.3) that𝑎=1,𝑏=1, (2.9) from the linear system (2.7), we have𝑑𝑥𝑑𝑡=𝑥−𝑦,𝑑𝑦𝑑𝑡=𝛼𝑥−𝛽𝑦. (2.10a)(2.10b) There is only one equilibrium point at the origin: (𝑥,𝑦)=(0,0). (2.11) The eigenvalues are 𝜆_1,2=12(1−𝛽±√Δ),Δ≡1−4𝛼+2𝛽+𝛽². (2.12) The parametric domain of feasible values for 𝛼 and 𝛽 is defined according to (2.5) and by the parabolic curve as in Figure 1: 1𝛼=4(𝛽+1)². (2.13) Thus we have 𝜆₁=𝜆₂𝛽∈ℝ,𝛼=²4+𝛽2+1, (2.14) where 𝛽>1, stable equilibrium in (0,0), and 𝛽<1, unstable equilibrium in (0,0); 𝜆₁≠𝜆₂𝛽∈ℝ,𝛼>²4+𝛽2+1, (2.15) where 𝜆₁<0,𝜆₂<0, stable equilibrium in (0,0),𝜆₁<0,𝜆₂>0, unstable equilibrium in (0,0) (saddle), and 𝜆₁>0,𝜆₂>0, unstable equilibrium in (0,0); 𝜆₁=𝜆₂𝛽∈ℂ,𝛼<²4+𝛽2+1, (2.16) where 𝛽>1,stable equilibrium in (0,0), and 𝛽<1, unstable equilibrium in (0,0).

Figure 1 

Parametric domain for the linear FHN model (8).

As a special case, by assuming that 𝑡₀=0,𝑥(0)=𝑥₀,and𝑦(0)=𝑦₀ with 𝛼=0, 𝛽≠0, from (2.10b) we have 𝑦(𝑡)=𝑦₀𝑒^−𝛽𝑡, (2.17) and from (2.10a) we have 𝑥𝑥(𝑡)=⁰𝑒^𝑡+𝑦₀1+𝛽[𝑒^−𝛽𝑡−𝑒^𝑡𝑒],𝛽≠−1,_𝑡(𝑥₀−𝑦₀𝑡),𝛽=−1. (2.18) Analogously, by assuming that 𝑡₀=0,𝑥(0)=𝑥₀,and𝑦(0)=𝑦₀ with 𝛽=0, 𝛼≠0, from (2.10b) we have (𝛼≠1/4) 𝑒𝑥(𝑡)=^𝑡/2√1−4𝛼[𝑥₀√√1−4𝛼cosh(1−4𝛼2𝑡)+(𝑥₀−2𝑦₀√)sinh(1−4𝛼2𝑒𝑡)],𝑦(𝑡)=^𝑡/2√1−4𝛼[𝑦₀√√1−4𝛼cosh(1−4𝛼2𝑡)+(2𝛼𝑥₀−𝑦₀√)sinh(1−4𝛼2𝑒𝑡)],𝑥(𝑡)=^𝑡/22[𝑥₀(2+𝑡)−2𝑦₀𝑒𝑡],𝑦(𝑡)=^𝑡/24[𝑥₀𝑡−2𝑦₀(𝑡−2)], (2.19) when 𝛼=1/4.

2.1.2. Third-Order Approximation

Assuming that𝑎=1,𝑏=1, (2.20) from the third-order system (2.8), we have𝑑𝑥𝑥𝑑𝑡=(𝑥−³6)−𝑦,𝑑𝑦𝑑𝑡=𝛼𝑥−𝛽𝑦. (2.21) The null clines intersect at the points (𝑥₁,𝑦₁()=(0,0),𝑥₂,𝑦₂)=(−6(𝛽−𝛼)𝛽𝛼,−𝛽6(𝛽−𝛼)𝛽(),𝑥₃,𝑦₃)=(6(𝛽−𝛼)𝛽,𝛼𝛽6(𝛽−𝛼)𝛽), (2.22) so that, being 𝛽>0, there are 3 disjoint equilibrium points (Figure 2) when 𝛽>𝛼, (2.23) and only one when 𝛼≤−𝛽. (2.24) The Jacobian of (2.21) is 𝑥𝐽(𝑥,𝑦)=1−²2−1𝛼−𝛽 (2.25) so that 𝐽(𝑥₂,𝑦₂)=𝐽(𝑥₃,𝑦₃1)=𝛽(3𝛼−2𝛽)−1𝛼−𝛽. (2.26) The eigenvalues are𝜆_1,2=12𝛽(3𝛼−2𝛽−𝛽²±Δ₁),Δ₁≡(3𝛼−2𝛽)²+𝛽²(𝛽²−4𝛽+2𝛼). (2.27) The parametric domain of feasible values for 𝛼 and 𝛽 is defined by the curve (3𝛼−2𝛽)²+𝛽²(𝛽²−4𝛽+2𝛼)=0 (2.28) as in Figure 3.

Figure 2 

Null clines for (2.21).

Figure 3 

Parametric domain of (2.21).

2.1.3. The General Case (2.4)

When 𝛾=0, assuming that 𝑎=1,𝑏=1,𝑓(𝑥)=𝑥,𝑔(𝑦)=𝑦 (2.29) from the system (2.4), we have 𝑑𝑥𝑑𝑡=sin𝑥−𝑦,𝑑𝑦𝑑𝑡=𝛼𝑥−𝛽𝑦. (2.30) The null clines have at least two intersections (Figure 4) when 𝛼≤𝛽or𝛽≥−4.6𝛼 (2.31) andonly one when 𝛼≤𝛽. (2.32) otherwise, there are more than one intersection.

Figure 4 

Null clines for (2.8).

The Jacobian of (2.30) and eigenvalues coincide with the linear case.

3. FHN Model with a High-Amplitude Spike

Let us consider the general nonautonomous system (2.4) when 𝛾≠0:𝑑𝑥𝑑𝑡=sin𝑥−𝑦+𝛾𝑧(𝑡),𝑑𝑦𝑑𝑡=𝛼𝑥−𝛽𝑦. (3.1) The time-dependent function 𝑧(𝑡) is a high-amplitude function (Figure 5), based on Haar wavelets 𝑧(𝑡)≡³_𝑘=−310𝜓⁴₈(𝑡−𝑘)+5𝜓⁴₉(𝑡−𝑘)+𝜓⁴₁0(𝑡−𝑘) (3.2) which replicates itself in a finite interval (Figure 6).

Figure 5 

High-amplitude function $10{\psi }_8^4(t)+5{\psi }_9^4(t)+{\psi }_{10}^4(t)$ .

Figure 6 

$z(t)$ .

The basic wavelet function 𝜓^𝑛_𝑘(𝑡) is defined as 𝜓^𝑛_𝑘(𝑡)≡2^𝑛/2𝑘1,2^𝑛≤𝑡<𝑘+0.52^𝑛,−1,𝑘+0.52^𝑛≤𝑡<𝑘+12^𝑛,0,elsewhere. (3.3)

3.1. Numerical Simulation

Let us assume for the values of parameters the following: 𝛾=10², 𝛽=2 and two different values for 𝛼, that is, 𝛼=−0.01, 𝛼=2 in the interval 𝑡∈[0,8]. As initial conditions it is assumed that 𝑥(0)=0, 𝑦(0)=0. By using the Runge-Kutta 4th-order method, with the accuracy 10⁻⁶, we obtain in the interval (0<𝑡≤8), the solution in correspondence with different initial conditions.

From a direct inspection of the solution (see Figures 7 and 8) it can be seen that under a spike firing there is some delay effect so that the perturbation show its influence only after some time delay. This perturbation of the system makes the orbits nearby the origin unstable. The uniqueness in phase space is going to be lost at least in the initial time interval. The orbit in the phase plane is self-crossing. The uniqueness of motion is lost. When 𝛼=−0.01, the origin behaves as an attractor, while in the case 𝛼=2, 𝑥, 𝑦 diverge asymptotically.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7 

Numerical solution of (3.1) with $\gamma =10^2\text{,}$$\beta =2$, and $\alpha =-0.01$ (a, c, e), $\alpha =2$ (b, d, f) in the interval $t\in [\mathrm{0,8}]$; the initial conditions are $x_0=0\text{,}$$y_0=0$.

(a)

(b)

Figure 8 

Phase orbits of (3.1) with $\gamma =10^2\text{,}$$\beta =2$,and $\alpha =-0.01$ (a), $\alpha =2$ (b) in the interval $t\in [\mathrm{0,8}]$; the initial conditions are $x_0=0\text{,}$$y_0=0$.

4. Wavelet Analysis

The Haar scaling function 𝜑(𝑡) is the characteristic function on [0,1]. By translation and dilation we get the family of functions defined (in [0,1]) as𝜑^𝑛_𝑘(𝑡)≡2^𝑛/2𝜑(2^𝑛𝑡−𝑘),(0≤𝑛,0≤𝑘≤2^𝑛−1),𝜑(2^𝑛𝑡−𝑘)=1,𝑡∈Ω^𝑛_𝑘,0,𝑡∉Ω^𝑛_𝑘,Ω^𝑛_𝑘𝑘≡[2^𝑛,𝑘+12^𝑛). (4.1) The Haar wavelet family {𝜓^𝑛_𝑘(𝑡)} is the orthonormal basis [13]: 𝜓^𝑛_𝑘(𝑡)≡2^𝑛/2𝜓(2^𝑛𝑡−𝑘),‖𝜓^𝑛_𝑘(𝑡)‖_𝐿2=1,𝜓(2^𝑛𝑘𝑡−𝑘)≡−1,𝑡∈[2^𝑛,𝑘+1/22^𝑛),1,𝑡∈[𝑘+1/22^𝑛,𝑘+12^𝑛),0,elsewhere.(0≤𝑛,0≤𝑘≤2^𝑛−1), (4.2) Without loss of generality, we can restrict ourselves to 0≤𝑛, 0≤𝑘≤2^𝑛−1⇒Ω^𝑛_𝑘⊆[0,1]. Let 𝐘≡{𝑌_𝑖}, (𝑖=0,…,2^𝑀−1,2^𝑀=𝑁<∞,𝑀∈ℕ), be a finite energy time-series; 𝑡_𝑖=𝑖/(2^𝑀−1), is the regular equispaced grid of dyadic points.

The discrete Haar wavelet transform is the operator 𝒲^𝑁 which maps the vector 𝐘 into the vector of the wavelet coefficients {𝛼,𝛽^𝑛_𝑘}: 𝒲^𝑁𝐘={𝛼,𝛽⁰₀,…,𝛽₂^𝑀−1_𝑀−1−1},𝐘={𝑌₀,𝑌₁,…,𝑌_𝑁−1}, (4.3) However in order to reduce the computational complexity it has been proposed [14] the short wavelet transform as follows. Let the set 𝐘={𝑌_𝑖} of 𝑁 data be segmented into 𝜎 segments (in general) of different length. Each segment 𝐘^𝑠,𝑠=0,…,𝜎−1 is made of 𝑝_𝑠=2^𝑚𝑠, (∑_𝑠𝑝_𝑠=𝑁), data: 𝐘={𝑌_𝑖}_{𝑖=0,…,𝑁−1}=^𝜎−1_𝑠=0{𝐘^𝑠},𝐘^𝑠≡{𝑌_𝑠𝑝𝑠,𝑌_{𝑠𝑝𝑠+1},…,𝑌_{𝑠𝑝𝑠+𝑝𝑠−1}}, (4.4) being, in general, 𝑝_𝑠≠𝑝_𝑟. The short discrete Haar wavelet transform of 𝐘 is (see [14]) 𝒲^{𝑝𝑠,𝜎}𝐘, 𝒲^{𝑝𝑠,𝜎}≡^𝜎−1_𝑠=0𝒲^𝑝_𝑠,𝐘=^𝜎−1_𝑠=0𝐘^𝑠,𝒲^{𝑝𝑠,𝜎}𝐘=(^𝜎−1_𝑠=0𝒲^𝑝𝑠)𝐘=(^𝜎−1_𝑠=0𝒲^𝑝𝑠𝐘^𝑠𝒲),_2𝑚𝑠𝐘^𝑠={𝛼₀^0(𝑠),𝛽₀^0(𝑠),𝛽₀^1(𝑠),𝛽₁^1(𝑠),…,𝛽^𝑚𝑠₂^−1(𝑠)_{𝑚𝑠−1−1}}, (4.5) with 2^𝑚𝑠=𝑝_𝑠,∑^𝜎−1_𝑠=0𝑝_𝑠=𝑁.

There follows that, the matrix of the wavelet transform is expressed as a direct sum of lower-order matrices so that the short transform is a sparse matrix [14]. When the short wavelet transform maps short interval values into a few set of wavelet coefficients, it can be considered as a first-order approximation. However, since the wavelet transform maps the original signal into uncorrelated sequences, the short wavelet transform describes, for each sequence of detail coefficients, its local behavior. When 𝑝_𝑠=𝑝=𝑁,𝜎=1, the above coincides with the ordinary wavelet transform. We assume, in the following, that 𝑝_𝑠=𝑝=𝑁/𝜎,𝑠=0,…,𝜎−1,(𝜎>1).

4.1. Critical Analysis through the Wavelet Coefficients

Let us take the short Haar wavelet transform for the two time series, obtained by using the Runge-Kutta 4th order method, with the accuracy 10⁻⁶, of system (3.1) with the following values of parameters-initial conditions: 𝛾=10²,𝛽=2,𝑥₀=0,𝑦₀=0,𝑡∈[0,8]. (4.6) The two time series correspond to the two values of 𝛼∶𝛼=−0.01 and 𝛼=2. It can be seen from Figures 9 and 10 that the jump is more visible in some coefficients. For instance both in Figure 9 and in Figure 10 the highest value of the amplitude of wavelet coefficients for 𝑥(𝑡) is in 𝛽¹₀ while for 𝑦(𝑡) is in 𝛽⁰₀. Probably this difference is due to the fact that 𝑥(𝑡) is less regular than 𝑦(𝑡).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9 

Wavelet coefficients of (3.1) with $\gamma =10^2\text{,}$$\beta =2$, and $\alpha =-0.01$ in the interval $t\in [\mathrm{0,8}]$; the initial conditions are $x_0=0\text{,}$$y_0=0$.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10 

Wavelet coefficients of (3.1) with $\gamma =10^2\text{,}$$\beta =2$, and $\alpha =2$ in the interval $t\in [\mathrm{0,8}]$; the initial conditions are $x_0=0\text{,}$$y_0=0$.

5. Conclusion

In this paper the FitzHugh-Nagumo model has been considered with a high-amplitude pulse. The analysis was done by using wavelet coefficients and it has been shown that the dynamical system shows some kind of inertia against the rapid jumps. In fact, the jumps are detected with some delay with respect to the time they appear.