FITZHUGH-NAGUMO EQUATIONS WITH GENERALIZED DIFFUSIVE COUPLING

Abstract. The aim of this work is to investigate the dynamics of a neural network, in which neurons, individually described by the FitzHugh-Nagumo model, are coupled by a generalized diffusive term. The formulation we are going to exploit is based on the general framework of graph theory. With the aim of defining the connection structure among the excitable elements, the discrete Laplacian matrix plays a fundamental role. In fact, it allows us to model the instantaneous propagation of signals between neurons, which need not be physically close to each other. This approach enables us to address three fundamental issues. Firstly, each neuron is described using the well-known FitzHugh-Nagumo model which might allow to differentiate their individual behaviour. Furthermore, exploiting the Laplacian matrix, a well defined connection structure is formalized. Finally, random networks and an ensemble of excitatory and inhibitory synapses are considered. Several simulations are performed to graphically present how dynamics within a network evolve. Thanks to an appropriate initial stimulus a wave is created: it propagates in a self-sustained way through the whole set of neurons. A novel graphical representation of the dynamics is shown.

1. Introduction.Signal dynamics within neural populations have received much attention in the past fifty years.How to tackle the issue of describing biological neural networks involves two master steps.Firstly, the choice of the model which describes each individual neuron has to be carried out.Secondly, the definition of interactions among all neurons of the network has to be explained and discussed.
Concerning the first point, the FitzHugh-Nagumo model [4] began as a dimensional reduction of the well-known Hodgkin-Huxley model [5].It extracts the Hodgkin-Huxley fast-slow phase plane and presents it in a simplified form.The resulting model is more analytically and numerically tractable and it maintains a certain biophysical meaning.Thus, the model is constituted by two equations in two variables v and r.The first is the fast variable called excitatory: it represents the transmembrane voltage.The second variable is the slow recovery variable: it describes the time dependence of several physical quantities, such as the electrical conductance of the ion currents across the membrane.The FitzHugh-Nagumo 204 ANNA CATTANI equations, using the notation in [6], are shown below: where, a, b, g ∈ R + are parameters of the model.The model describes neurons as excitable elements which have two key properties.Firstly, they are characterized by their excitability behaviour: a sufficiently large stimulus provokes a very large response, that is, a small perturbation to the quiescent state of a neuron can provoke a large excursion of its potential.Secondly, they are characterized by their refractoriness: the elements cannot be excited during the period which follows the stimulus.
Once the choice of the model for each individual neuron has been made, the challenge is how to describe the dynamics of the action potential within a network.We sustain that the graph theory is the most suitable context to tackle this issue.According to this choice we are able to establish that links among neurons satisfy a specific rule in which also not near neighbour interactions are allowed.Not near neighbour interactions imply that we are allowing axons to reach neurons that are far from the ones that generated the spike, as is typical in the brain.
The aim of this work is to present a description of a neural network and the dynamics of the action potential within it, where neurons, individually described by the FitzHugh-Nagumo model (1), are coupled by a generalized diffusive term.The word "generalized" underlines the previously specified key property of the model, i.e., the exchange of signals may exist among neurons even if they are not adjacent.A thorough description of the model in the framework of the graph theory is the starting point of the present work.Furthermore, the study of the stability of the solution's equilibrium point is performed in the cases of one and two coupled neurons.However, the core of the paper concerns the visualizations of the solutions obtained by exploiting different connection rules among neurons.In particular, a novel representation is shown.

Excitable neural network.
To be able to present the FitzHugh-Nagumo model relating to a neural network, it is fundamental to introduce basic concepts proper to graph theory.
Let us consider a graph G = (V, E), where V = {1, • • • , N } ⊂ N is the set of vertices and E ⊂ V × V is the set of edges.The so-called adjacency matrix A G = [a ij ] is a N × N matrix whose entries are: where i, j = 1, • • • , N and weights w ij = ±1, as proposed in [1].Moreover, the diagonal degree matrix Furthermore, we introduce the concept of Laplacian matrix which is defined as follow: (2) Hence, the Laplacian matrix is an N × N matrix whose elements are Arrangement of the graph vertices in the plane.Neuron labels are disposed over the nodes of a one dimensional closed chain, i.e., a ring.Graph edges are not represented.An example of connection structure over the ring is described in Figure 2.
Exploiting the notions proper to graph theory, let us tackle the issue of describing a neural network.Neurons are identified with integer labels 1, • • • , N which are collected in the set of vertices V while the links among them constitute the set of edges.According to the definition of the Laplacian matrix (3), we are able to describe the connection structure among neurons: entry w ij represents the presence of a synapse between neurons i and j.The weights w ij = ±1 allow both excitatory (if w ij = 1) and inhibitory (if w ij = −1) synapses.Although Eq. ( 3) is the general definition of the Laplacian matrix, the connection structures we will exploit allow us to consider a more specific expression of the equation.Some further assumptions have to be established.Firstly, let us assume that each neuron is linked with a finite number of others.In particular, we will focus on the case of two connections per neuron.These connections are invariant under discrete neuron labels translations, i.e., we keep the same connection rule for each neuron.Secondly, let us consider that neuron labels are disposed over the nodes of a one dimensional chain with periodicity, i.e., a ring as shown in Figure 1.These ingredients mean that a sparse Laplacian matrix with the following entries is considered: where (5) Hence, in Eq. (4) the number of connections for all neuron i is shown on the diagonal while the linked neurons which take effect on it are neurons i − q and i + k with q, k ∈ N. The function φ in Eq. ( 5) describes the periodicity of the ring.Let us underlined that Eq. ( 4) is, in general, a non-symmetric matrix whenever q = k.The case of four connections per neuron is studied in [3].
In order to present the model that we will exploit in the simulations, let us make some biophysical assumptions.Firstly, we hypothesize that neurons are identical entities, i.e., they show the same dynamics when stimulated in the same way.Furthermore, we suppose infinite velocity in the exchange of communications among An example of link structure over the graph ring is represented.The connection structure is described by the Laplacian matrix (6).neurons.A future work will tackle the issue of delay in synapses.Furthermore, signals can reach even not neighbouring neurons without passing through the adjacent ones.This means that, in accordance with the connection structure which exists in a specific neural network, direct links between vertices arise in Figure 1.For example, if q = k = 2 in Eq. ( 4), the symmetric Laplacian matrix is described as follow: and the corresponding graph is shown in Figure 2. As a result of this treatment, we are able to formalise the model which describes the dynamics of signals within a general neural network.Hence, for each neuron i = 1, • • • , N within the network, the FitzHugh-Nagumo model with generalized diffusive coupling can be expressed as below: with parameters a, b, g, d ∈ R + .Exploiting the adjacency matrix instead of the Laplacian matrix, an equivalent formulation is: in which the presence of the well-known diffusive coupling in d j a ij (v j − v i ) becomes more evident.
3. Analytical results.Before showing the dynamics produced by Eq. ( 7), we present some analytical results.In this section we introduce two propositions which ensure the stability of the equilibrium point in the case of an uncoupled neuron and two coupled neurons whatever the choice of parameters a, b, g, d ∈ R + is.However, since we are interested in restricting our analysis to the systems which exhibit excitability with only one equilibrium point at the origin, let us state that: Proposition 1.The equilibrium point of the systems Eq.(1) and Eq. ( 11) is unique and it is the trivial one if and only if The proof follows from straightforward computations and is hence omitted.
For the sake of completeness, both systems ( 1)-(11) exhibit a periodic limit cycle by adding an external current to the v differential equation.
3.1.Single neuron.Consider first the case of a single uncoupled neuron whose dynamics is described by Eq. ( 1).Point is the trivial stationary solution and it is the unique equilibrium point since Eq. ( 9) holds.
Proof.The linearization of Eq. ( 1) on the trivial equilibrium point (10 The corresponding characteristic equation is which provides the eigenvalues of the linearized system: As ∀a, b, g ∈ R + is Re(λ 1, 2 ) < 0, then we conclude that (v, r) = (0, 0) is stable.
3.2.Two coupled neurons.Let us assume that neurons are identical.The model (8) with Let us note that there are two different diffusion coefficients d 1 , d 2 > 0. Although in the next section we will exploit a unique diffusion coefficient for each coupling, this choice ensures the possibility to consider an asymmetric connection structure even in the case of two coupled neuron.
Proof.The steps of the proof are the same as in the case of a single neuron.Here, the linearization of (11) with respect to the trivial equilibrium point is the following: As ∀a, b, g, d ∈ R + is Re(λ 1, 2, 3, 4 ) < 0, then we conclude that (v 1 , r1 , v2 , r2 ) = (0, 0, 0, 0) is stable.
As in [2], the analytical results concerning the stability of the trivial equilibrium point are performed when N = 1 and N = 2.When N > 2, we only check stability by numerical integrations.
4. Numerical integrations.The aim of this section is to present several dynamics of the action potential in a set of N excitable FitzHugh-Nagumo elements with the diffusive coupling described by Eq. (7).Integrations are done using Matlab.The differential equations are advanced in time by the Runge-Kutta method (ode45 solver in Matlab, with default parameters).The parameters used in the following simulations are a = 0.25, b = 0.001, g = 0.003, d = 0.05 and N = 128.Although in the following dynamics N > 2, it is possible to check stability by analysing the simulations at large time from a qualitative point of view.In fact, dealing with this specific choice of parameters, at the end of the dynamics all neuron action potentials will return to the quiescent state, i.e., (v f , r f ) = (0, 0).
In order to highlight the dynamics we are going to show, neurons are disposed in line with increasing integer labels.The periodic boundary conditions are still considered.Below, we show several dynamics in which different connection rule are taken into account.The next two subsections show, firstly, the dynamics produced by all excitatory synapses and, secondly, those produced by an ensemble of excitatory and inhibitory synapses.4.1.All excitatory synapses.In this subsection we will consider the case of all excitatory synapses between linked neurons.This means that w ij = 1 in Eq. ( 4).In all the following dynamics, an initial stimulus is applied to the central neuron of the line.The stimulus consists in imposing a non-null initial action potential on the central node and we set (v 0 , r 0 ) = (0.5, 0).As we will see, different connection structures, which means different Laplacian matrices, produce a variety of behaviours in the whole set of neurons.However, whatever the dynamics produced, all neurons return to the quiescent state at the end.Neurons are modelled as excitable units and then, after the excitation, they undergo a long refractory period.In this period they are blind to any stimulus.This is the reason why, as we will see in the dynamics frames, two travelling waves that collide depress their signals.Let us now go into details of the specific dynamics.7) with the Laplacian discretization (4) where q = k = 1.In this simulation N = 128 and an external stimulus is applied to the central neuron.Let us underline that the external stimulus is produced by a non-quiescent initial state, i.e., (v 0 , r 0 )(N/2) = (0.5, 0).The rest values are v 0 = 0 and r 0 = 0.At the end of the dynamics all neurons the whole set of neurons is at the quiescent state.
Let us first consider equations ( 4)-( 7) with q = k = 1.Naively, a neuron which receives a signal fires in turn to the two adjacent ones (to the two neurons identified by the adjacent integer labels).Then, after having applied an initial stimulus to the central neuron of the line, the result is the propagation of two impulses away from that point of stimulation.The boundary conditions are considered periodic as described in Eq. ( 5).Some selected frames of the dynamics are presented in Figure 3.An alternative representation is shown in Figure 5.To graphically stress how neurons reach the quiescent state at the end of integrations, let us consider a non-null initial datum on the first neuron of the line.Through the periodic boundary conditions the waves travel in opposite directions and they collide in the center of the domain.In Figure 4 the dynamics is reproduced.The same phenomenon can be observed by analysing Figure 5 for the upper bound values of t where the action potential of the neurons returns to zero.Let us underline that, in this simple case, the dynamics produced coincides with those proposed by [4] and recalled in [8].In these two references, the propagation along an axon is described.Roughly speaking, an axon can be assimilate to our set of neurons disposed in line.
Differently from the case q = k = 1, by exploiting equations ( 4)-( 7) with q = k = 2 we obtain a dynamics in which odd-label neurons remain at the quiescent state.It follows that, if we consider a set of N = 64 neurons constituted by even-label ones, the dynamics is equivalent to the one presented in the case of q = k = 1.
The case of equations ( 4)-( 7) with q = k = 3 is shown in Figure 6.The dynamics involves all neurons due to this choice.Differently from the case q = k = 1, before waves reach the boundary, two on three neurons remain at the quiescent state.Even if at different integration times, due to the periodic boundary conditions, all neurons  7) with the Laplacian discretization ( 4) where q = k = 1.In this simulation N = 128 and an external stimulus is applied to the central neuron.Let us underline that the external stimulus is produced by a non-quiescent initial state, i.e., (v 0 , r 0 )(1) = (0.5, 0).The rest values are v 0 = 0 and r 0 = 0. Due to the periodic boundary conditions, when the waves collide they disappear.In fact, neurons which are in the refractory period, i.e. the period that follows the instant when they reach the maximum action potential, cannot propagate the signals in turn and cause the decay of the action potential.At the end of the dynamics all neurons are at the quiescent state.
excite and then return to the stable equilibrium point.A correct permutation of neurons makes it possible to obtain the dynamics of q = k = 1.The dynamics proposed in Figure 7 and Figure 8 are obtained by exploiting equations ( 4)-( 7) with, respectively, (q = 1, k = 2) and (q = 5, k = 2).In contrast with the previous cases, two asymmetric dynamics are produced.Nevertheless, all neurons are involved in the dynamics.
A way to generalize the dynamics between neurons is to make the structure of the connections random, i.e., to randomly construct the Laplacian matrix.Specifically, let us firstly define for each neuron the connection number n i such that it is uniformly distributed on number two, three and four; in symbols: n i ∼ U(2, 4).Then, for each neuron i, the integer-labels of the n i links are determined exploiting the normal distribution.Specifically, where X i is a vector of n i components.Imposing µ = 8 and σ 2 = 25, an example of Laplacian matrix, having definition presented in Eq. ( 4), is shown in Figure 9.The resulting dynamics is presented in Figure 10.In accordance with the Laplacian matrix, the wave in the dynamics travels to the right.

4.2.
Ensemble of excitatory and inhibitory synapses.In this section, both excitatory and inhibitory synapses are considered.This ensemble is a fundamental ingredient to make the model meaningful from a biological perspective.In order to  7) with the Laplacian discretization (4) where q = k = 3.The parameters and the initial datum are the same used in Figure 3. Differently from the case of q = k = 2, all neurons will be excited.All neurons will return to the rest state (v, r) = (0, 0) at the end of the integration.describe a certain number of inhibitory connections, several sub-diagonal and superdiagonal entries w ij of the Laplacian matrix will be equal to −1.This approach is proposed, among others, in [7].Specifically, if a row i exists such that l i,i+k = 1, then the link between neurons i and i + k produces an inhibitory synapse.In Integration of Eq. ( 7) with the Laplacian discretization (4) where q = 1 and k = 2.The initial datum is (v 0 , r 0 )(N/2) = (0.5, 0).At the end of the integration all neurons will return to the quiescent state (v, r) = (0, 0).7) with the Laplacian discretization (4) where q = 5 and k = 2.The initial data are (v 0 , r 0 )(N/2) = (0.5, 0).At the end of the integration all neurons will return to the quiescent state (v, r) = (0, 0).from the left.Neurons 9, 24 receive inhibitory synapses from the right.In Figure 12 a comparison between the dynamics shown in Figures 3-11 is presented.It is important to underline that the action potential of neurons which receive an inhibitory synapse reaches lower values than excitatory synapses.Moreover, throughout the whole dynamics, inhibitory synapses produce a slower wave.Both phenomena can be observed in Figure 12.As underlined in the case of all excitatory synapses, all neurons return to the resting state at the end of the integration.Due to the periodicity of boundary conditions, the travelling waves collide.Neurons in the refractory  7) with the Laplacian discretization (4) built up as explained above in the text.The parameters and the initial datum are the same as in the previous dynamics.Accordingly to Eq. ( 12), the dynamics describes a wave which travel on the right.At the end of the integration, each neuron returns to the resting state.
period cannot propagate the signals in turn and cause the decay of the action potential in the whole set of neurons.Let us underline that, if we had assumed that four adjacent-labelled neurons receive an inhibitory synapse, the signal would abruptly disappear and neurons would suddenly return to their quiescent state.

5.
Conclusions.In summary, we have proposed the model (7) which is able to describe the dynamics of the action potential and the recovery variable within a  7) with the Laplacian discretization (4) with q = k = 1.In Eq. ( 4), admitted weights w ij are +1 and −1; this translates in considering both excitatory and inhibitory synapses.neural network (Sec II).A thoroughly explained mathematical structure allows us to formally describe several fundamental features of interactions in neural populations.
In fact, the description of not only near neighbor interactions and the presence of inhibitory synapses has been achieved.Then, the stability of the equilibrium point of the solution has been performed in two sample model cases (Sec III).To provide a graphical explanation of the model with different connection structures, several frames of dynamics have been shown (Sec IV).A novel representation offers a relevant understanding of the solution provided by the model.
Future work should explore some important aspects concerning the model.First of all, a way to make the model more biophysically relevant is to consider a nonnegligible time delay in coupling.After having introduced this ingredient we will be able to compare the results provided by the model with in-vitro experiments.Moreover, a currently ongoing work deals with the solutions obtained by exploiting the model with the number of neurons that tends to infinity in a bounded area.We expect this technique to open up a new way to study signal dynamics within large populations of neurons.

Figure 3 .
Figure 3. Integration of Eq. (7) with the Laplacian discretization (4) where q = k = 1.In this simulation N = 128 and an external stimulus is applied to the central neuron.Let us underline that the external stimulus is produced by a non-quiescent initial state, i.e., (v 0 , r 0 )(N/2) = (0.5, 0).The rest values are v 0 = 0 and r 0 = 0.At the end of the dynamics all neurons the whole set of neurons is at the quiescent state.

Figure 4 .
Figure 4. Integration of Eq. (7) with the Laplacian discretization (4) where q = k = 1.In this simulation N = 128 and an external stimulus is applied to the central neuron.Let us underline that the external stimulus is produced by a non-quiescent initial state, i.e., (v 0 , r 0 )(1) = (0.5, 0).The rest values are v 0 = 0 and r 0 = 0. Due to the periodic boundary conditions, when the waves collide they disappear.In fact, neurons which are in the refractory period, i.e. the period that follows the instant when they reach the maximum action potential, cannot propagate the signals in turn and cause the decay of the action potential.At the end of the dynamics all neurons are at the quiescent state.

Figure 5 .Figure 6 .
Figure 5.An alternative representation of the dynamics shown in Figure3is proposed.The action potential v(t, x) for each neuron, disposed over the y-axis, by time increasing, displayed on the x-axis, is described.

Figure 9 .Figure 10 .
Figure 9. Non-zero elements in the Laplacian matrix which produces the dynamics shown in Figure 10.

Figure 12 .
Figure 12.Comparison between dynamics shown in Figures 11-3.Black dots represent the dynamics with several inhibitory synapses.Red diamonds describe the the dynamics with only excitatory synapses.