The derivation of continuum limits of neuronal networks with gap-junction couplings

We consider an idealized network, formed by N neurons individually described by the FitzHugh-Nagumo equations and connected by electrical synapses. The limit for N to infinity of the resulting discrete model is thoroughly investigated, with the aim of identifying a model for a continuum of neurons having an equivalent behaviour. Two strategies for passing to the limit are analysed: i) a more conventional approach, based on a fixed nearest-neighbour connection topology accompanied by a suitable scaling of the diffusion coefficients; ii) a new approach, in which the number of connections to any given neuron varies with N according to a precise law, which simultaneously guarantees the non-triviality of the limit and the locality of neuronal interactions. Both approaches yield in the limit a pde-based model, in which the distribution of action potential obeys a nonlinear reaction-convection-diffusion equation; convection accounts for the possible lack of symmetry in the connection topology. Several convergence issues are discussed, both theoretically and numerically.


Introduction
The computer simulation of the behaviour of complex networks with a huge number of nodes, such as networks of neurons in some portion of the brain, is a formidable challenge. The intrinsic difficulties of such a task may be alleviated to some extent by identifying one or more multiscale structures within the networks; this allows one to describe and simulate different scales by different models, while posing the problem of the interaction among scales.
Within a multiscale framework, the co-existence of discrete and continuous models is a natural option, which may lead to significant savings. Higher-level nodes, or interactions, may be affordably given an individual description (e.g., by a system of coupled ordinary differential equations), whenever their number is small to moderate. On the contrary, this approach would be computationally prohibitive for the description of lower-level nodes or interactions, if their number is exceedingly large. In this case, a possible alternative may consist in modelling the huge population of individuals by a continuum, confined in some spatial region, and describe its behaviour by a limited number of variables, e.g., submitted to satisfy partial differential equations. One recognizes here a process underlying the mathematical description of many physical phenomena, e.g., in Fluid Dynamics.
The derivation of a continuum model may be accomplished by "passing to the limit" in a discrete model, assuming that the number of individuals tends to infinity. The present paper aims at investigating one such limit process. Specifically, the discrete model we start from is inspired by the connections of electrical (rather than chemical) nature within a neuronal network. In order to better motivate our model, we provide now some biological background.
Despite in the past years almost all networks have been represented as constituted by neurons that are interconnected by chemical synapses, electrical synapses are largely present in the nervous system. In the sequel, we will use indifferently the terms electrical synapses and gap junctions. However, for the sake of completeness, gap junctions are the morphological equivalent of electrical synapses. In particular, as specified in [8], gap junctions exist between near-neighbour neurons and they allow low-resistance electrical transmissions. Indeed, at an electrical synapse a current I gap is generated which is proportional to the difference between the action potentials v of the post-synaptic and pre-synaptic neurons (see, e.g., [3] eq. (7.12)); explicitly, we have for some d > 0 This establishes a diffusive coupling between neighbouring neurons. Unfortunately, the analysis of electrical synapses in situ presents severe technical difficulties and therefore their specific roles are still largely unexplored. Nevertheless, in the past ten years, the topic concerning gap junction networks has been object of several investigations, sometimes leading to paradoxical results (see e.g. [6,9,15]).
In order to build up the sample network we will consider, several ingredients are taken into account. First of all, we model each single cell as an excitable element by exploiting the FitzHugh-Nagumo model (see [4]). The excitable feature means that neurons may not fire intrinsically without any synaptic inputs. Furthermore, each cell belongs to the same functional class, avoiding the presence of heterogeneity. This agrees with authors in [6] who stress that electrical synapses exist exclusively between neurons of a specific class. In particular, despite many works underline the presence of electrical synapses between inhibitory neurons (see e.g. [6]), the existence of electrical connections between excitatory neurons is demonstrated in the early postnatal stages (see [15]). Finally, as we will specify, we consider the presence of both bidirectional (non-rectifying) and unidirectional (rectifying) synapses as claimed in [9].
We now describe the content of this paper. After setting our mathematical model of an idealized neuronal network with electrical-type coupling between neurons, we carefully investigate the "passage to the limit" as the number of neurons tends to infinity, while they remain confined in a fixed and bounded spatial region. We identify two different manners of increasing the population of the network so that a non-trivial continuum limit is obtained. The first one assumes a fixed topology of the network (nearest-neighbour connections) but makes the proportionality coefficient in (1.1) to depend upon the total number of neurons according to a specific law; conversely, the second manner keeps this coefficient fixed but suitably increases the number of connections per neuron. Both methods lead to equivalent continuous models, in which the action potential is the solution of a reaction-diffusion partial differential equation (or a reaction-convection-diffusion equation if connections are not symmetric, i.e., if rectifying synapses are allowed). Our arguments apply in any spatial dimension, although we detail them in 1D and we sketch their extension to 2D. Clear numerical evidence confirms all theoretical results. At last, an example of random connections is also presented.

The FitzHugh-Nagumo model of a single neuron
The FitzHugh-Nagumo model [4] was introduced as a dimensional reduction of the well-known Hodgkin-Huxley model [7]. 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. The model is constituted by two equations in two variables v and r. The first one 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 equations, using the notation in [10], are given by: where, a, b, c ∈ 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. Throughout the paper, following [14], we will adopt the FitzHugh-Nagumo model with a = 0.25, b = 0.001, c = 0.003, to describe the behaviour of each neuron in our network.

Diffusive coupling within the network
We suppose that our network contains N neurons, identified by integer labels i = 1, · · · , N ; labels may refer to the physical position of the neurons, but other ways to index neurons could be used in a more convenient way. Electricaltype connections in the neuronal network are easily described by basic concepts from graph theory (see, e.g., [1]). 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 an N × N matrix whose entries are: where i, j = 1, · · · , N and the weights are strictly positive. Exploiting the adjacency matrix, and assuming the gap-junction law (1.1) for the interaction between adjacent neurons, we define the FitzHugh-Nagumo model with diffusive coupling as follows: Specifically, the summation describes the influence on the i−th neuron of all neurons linked to it; it produces a diffusion effect within the network. The simplest example is given by the expression in (4.2), which models nearestneighbour interactions in a chain of neurons.
Introducing the diagonal degree matrix D G := diag(d i ) with d i = j =i a ij , and the Laplacian matrix L G := D G − A G = (l ij ), the previous system can be written asv In the sequel, we assume that all weights w ij are equal and precisely w ij = d for some d > 0, which we will call the diffusion coefficient. Let us introduce the set Q(i) of all indices q such that neuron i + q is linked to neuron i, i.e., a i,i+q = 0. Then, the model (3.1) can be written aṡ In most cases, we shall consider Q(i) = Q independent of i, thus assuming a homogeneous network topology.
We are interested in describing the behaviour of the network as the number of neurons increases, identifying conditions on the model which lead to nontrivial asymptotic patterns in the limit N → ∞. We assume that the network is contained in a bounded region B (independent of N ) of the Euclidean space R m , for some 1 ≤ m ≤ 3; let us denote by x i ∈ B the physical position of the i−th neuron. Then, we assume that the distance of any pointx ∈ B from the network tends to zero as N → ∞, and the distance of each neuron from its neighbours in the network has a similar behaviour.
If interactions between neurons are local, we can give an expression of the diffusive term in (3.2) which is based on the Taylor expansion of the differences ∆v i,q = v i+q − v i . Precisely, let us assume that at each time there exists a sufficiently smooth function where ∇v denotes the gradient vector of v, whereas Hv denotes the Hessian matrix of v. Substituting this expression into (3.2), we obtain a representation of the diffusive term by which we can find the conditions on the coefficient d and/or the sets Q(i) (depending on the network) yielding a non-trivial limit as N → ∞. We will detail our analysis assuming a specific distribution of neurons in the one-dimensional case first, and then we will consider the multidimensional extension.

One-dimensional dynamics
We consider neurons disposed over a closed chain, i.e., a ring. Each neuron occupies a specific physical position x i in the interval B = [0, 1] given by where N is the number of elements equally distributed along the chain and, consequently, ∆x = 1/N is the distance between any two adjacent ones. Since the chain is closed, we assume periodic boundary conditions, i.e., we set v 0 = v N and v i+kN = v i for any k ∈ Z.

Nearest-neighbour interactions
Let us first consider two symmetric nearest-neighbour interactions for each neuron. This translates in considering the set of connections per neuron Q(i) = Q = {±1}. In this case, the diffusive coupling assumes the following form: (4.2) An interesting dynamics produced by (3.2), which will represent a test case for the subsequent discussion, is obtained by applying an initial stimulus to the central neuron (i = N/2, assuming N even) of the line. Specifically, its action potential is initially set to the value 2, whereas all the other variables are set to 0. Considering the diffusion coefficient d = 0.05 (see [14]), the resulting dynamics is constituted by two pulses that travel in opposite directions in the whole set of neurons (see e.g. [8] for the analysis of travelling pulses). A sample dynamics is shown in Figure 4.1. We observe for further reference that a similar dynamics is obtained starting from an initial stimulus of the action potential given by a Gaussian function concentrated around the central neuron. In all cases, at the end of dynamics, neurons return to the quiescent state. In fact, 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 two travelling pulses that collide depress their signals.
We now focus on how the dynamics produced by our model depends upon N . The first observation is that, if the diffusion parameter d is kept fixed, then the diffusive effect tends to vanish as N → ∞. This can be seen in two ways. On the one hand, if neurons get close to each other and the action potential varies in a smooth manner, then the differences on the right-hand side of (4.2) tend to zero, implying the vanishing of the diffusion term L G v in each node. On the other hand, considering the test case introduced above, it is easily seen that the effect of, say, doubling N is equivalent to have a chain of neurons with the same spacing but with double length; this means that on the original chain, waves have half the length and propagate with half the speed. In order to obtain non-trivial diffusion effects in the limit, one possibility -that we call Approach I -consists in letting the parameter d grow with N , i.e., d = d N . The precise dependence can be found by exploiting the Taylor expansion (3.3), which in the present setting becomes where the prime indicates differentiation with respect to the spatial variable x. Therefore, the following expression holds for the diffusive term: for a fixed constant d * > 0. Hence, we obtain i.e., d N is proportional to the square of the number of neurons. The fact that d N is proportional to N 2 is not surprising: the spectral gap of the Laplacian matrix has the same behaviour as 1/(N 2 ). As N → ∞, the discrete model "converges" to a continuous model. To support this statement, we observe that the quantity h.o.t. in (4.4) is given by It follows that if we fix any pointx ∈ [0, 1] and, for each N , we consider a neuron of index i = i(N ) such that We conclude that a continuum of neurons is the results of the limit process of letting N → ∞, and is the system of nonlinear partial differential equations of incomplete parabolic type which describes the action potential and the recovery variable in the whole set of neurons. Note that the first equation is similar to the so-called cable equation, which describes the distribution of the potential along the axon of a single neuron (see, e.g. [3,12]). Reaction-diffusion models like (4.9) are studied e.g. in [5,13]. We observe that the discrete model (4.6)-(4.7) can be viewed as a numerical semi-discretization (in space) of the PDE system (4.9), obtained by using a second-order centered finite difference method on the equally-spaced (4.1). Thus, if the solution of (4.9) is sufficiently smooth as in the case of an initial Gaussian stimulus, we expect to have quadratic convergence in ∆x of the discrete solutions, at any fixed time t > 0, as it can be deduced from the fact that the error term on the right-hand side of (4.8) is proportional to ∆x 2 .
We now give an example. Following the choice of parameters presented in [14], we set d = 0.05 and we consider the case N = 128 as a reference one, i.e., we impose d N = d for N = 128, which yields d * = 0.05 128 2 = 3.0518 · 10 −6 . (4.10) A comparison of several discrete solutions is presented in Figure 4.2. The (b) plots clearly document the convergence of the discrete dynamics towards a limit one. Note that these dynamics are generated by applying an initial stimulus v i|t=0 = 2 to a number of neurons proportional to N around the center of the chain; in the limit, the initial action potential takes the value 2 in an interval of positive length symmetrically placed around the point x = 1/2, and vanishes elsewhere.
2) also depending on i and q, i.e., the diffusive coupling law is replaced by and d * is a smooth function. In the limit, the diffusion term in (4.9) is replaced by ∂ ∂x d * ∂v ∂x . For simplicity, we confine ourselves to the constant-coefficient case.

Extended range interactions
We now introduce a second approach to reproduce the same limit dynamics emerged above, which avoids rescaling the diffusion coefficient with the square of the number of neurons. This alternative way -which we call Approach II -consists of increasing the number of connections per neuron according to a specific law (and just slightly adjust the diffusion coefficient).
Since the core idea is to consider a number of connections per neuron that varies as a function of N , let us define the following set: where Q N is a positive integer to be determined. Thus, neurons linked to the i-th one belong to the interval Using again Taylor expansions, the sum in the diffusive coupling becomes Introducing the function ϕ : and invoking the identity We would like to choose Q N in such a way that for a fixed constant d * > 0. This equation admits a unique solution, say Q r N , which however need not be an integer. Therefore, we choose Q N as the nearest integer to Q r N . Proposition 4.2. The number of neurons linked to any given one grows proportionally to the power N 2/3 of the total number of neurons. Indeed 1 , 1 For any two non-negative sequences A N and B N , we will use the symbols Proof. By definition, Q r N satisfies (4.17) The result follows recalling that ϕ(x) ∼ x 3 3 for x → ∞. Let us underline that, although the number of neurons linked to any given one grows with N , interactions remain local, i.e., these neurons belong to a neighbourhood whose size decays with N . Indeed, considering the i−th neuron and recalling (4.13), we have Thus, we expect that the limit model, as N → ∞, be again described by partial differential equations. As specified above, the slight shift from Q r N to Q N provokes the necessity of slightly modifying the diffusion coefficient. Precisely, we define d N so that the identity is satisfied. An alternative possibility, which will be explored later on and which leads to similar effects, would be to define d * N so that The coefficient d N is really a small perturbation of d, as the next proposition indicates. Proof. From (4.17) and (4.19), we obtain the following equality: (4.21) Since Q N is defined as the nearest integer to Q r N , and then, Q N = Q r N + ε N with a proper choice of ε N , such that |ε N | ≤ 1/2. Writing ϕ(Q r N ) = ϕ(Q N ) + ϕ(Q r N ) − ϕ(Q N ) and substituting in (4.21), we get Using (4.22) and omitting computations, we conclude that |ϕ(Q r N ) − ϕ(Q N )| N 4/3 while ϕ(Q N ) ≃ N 2 . This gives the desired estimate.
In order to obtain the continuous model as a limit of the discrete model for N → ∞, we observe that, if the fourth derivative of v is continuous in [0, 1], the diffusion term (4.15) takes the form This means that Approach II yields in the limit the same system (4.9) of partial differential equations, that we got from Approach I. We now illustrate the asymptotic behaviour of the quantities defined above, for the same test case considered in the previous subsection. We choose again d = 0.05, and we enforce that for N = N 0 = 128 we have Q N0 = Q r N0 = 1, which corresponds to the nearest-neighbour interaction previously considered; we also enforce d N0 = d, and consequently we get which is precisely (4.10). Increasing N by powers of 2, i.e., setting N = N 0 2 p with p ≥ 1, the algorithm presented above produces the values of Q N and d N shown in Table 1. The last column of this table, as well as  Figure 4.5. While the shapes of the pulses are already well captured, their speed of propagation is less accurately reproduced; this should be related to the fourth-order error term on the right-hand side of (4.24), whose decay is slower than in Approach I as indicated by (4.25) compared to (4.8).

Non-symmetric interactions
A more general configuration of the network admits non-symmetric links for each neuron, which correspond to unidirectional connections (the so-called rectifying synapses). A natural extension of the symmetric case consists in choosing    where would be an obvious alternative.) We will prove that a suitable choice of Q C N depending on N leads to modify the limit model (4.9), by adding a first order term to the action potential equation. With our definitions, the sum in the diffusive coupling becomes Exploiting the Taylor expansion (4.3), we obtain

(4.28)
Recalling the definition (4.14) of the function ϕ, and introducing the function ψ : R + → R + defined as and such that n q=1 q = ψ(n), it is easily seen that the diffusive coupling takes the form Ideally, we would like to find integers Q D N and Q C N > Q D N satisfying the system   We observe that, given any arbitrary d * and c * , there always exists an integer N * such that condition (4.32) is satisfied for all N ≥ N * . Definition 4.5. Under the assumption (4.32), we define Q D N and Q C N , resp., to be the nearest integers to Q D,r N and Q C,r N , resp., which are the unique solutions of the system Proposition 4.6. The following asymptotic behaviour of the integers Q D N and Q C N holds: Proof. It is enough to estimatex = Q D,r N andŷ = Q C,r N . We recall that they satisfy (4.33). With the ansatzŷ ≃ N α , we have ϕ(ŷ) ≃ N 3α . On the other hand, the inequalityx <ŷ and the monotonicity of ϕ yield ϕ(ŷ) ≤ ϕ(x)+ϕ(ŷ) ≤ 2ϕ(ŷ). Since A N ≃ N 2 , we deduce that ϕ(ŷ) ≃ N 2 , whence α = 2/3. On the other hand, ψ(ŷ) ≃ N for somex <ẑ <ŷ; since ψ ′ (ẑ) =ẑ + 1/2 ≃ N 2 3 , we conclude thatŷ −x ≃ N Even for the present model, interactions are local. Indeed, all neurons linked to the i−th one belong to the interval . In order to accommodate the effect of the slight shift from (Q D,r N , Q C,r N ) to (Q D N , Q C N ), we introduce perturbations (d * N , c * N ) of (d * , c * ). They are defined in such a way that (Q D N , Q C N ) is the solution of the system (4.36) The size of the perturbation can be estimated as follows.
Proposition 4.7. The perturbed coefficients d * N and c * N introduced above satisfy Proof. Using (4.35) and (4.36), we get As in the proof of Proposition 4.3, we have |ϕ( , which gives the result. Finally, we study the limit behaviour of our model as N → ∞. To this end, we make use of the following expression for the higher order terms in (4.30): which holds under the assumption that the fourth derivative of v is continuous Thus, we obtain the following result.
Therefore, the discrete model (3.2) with Q given by (4.26) and Q D N , Q C N defined in Definition 4.5 leads for N → ∞ to the continuous model  i) Observe that having a larger number of neurons influencing a given neuron from its right rather than from its left results in a convective term, whose coefficient c * is positive; this corresponds to a negative speed of convective propagation, i.e., waves moving from right to left, as documented by Fig. 4.6. Obviously, choosing c * = 0 yields Q C N = ∅, so one is back to the symmetric case considered in Sect. 4.2.
ii) The same limit model can be obtained with a nearest-neighbour interaction that extends the one considered in Sect. 4.1, i.e.,

39)
with d N = d * N 2 and c N = c * N 2 . iii) A generalization to variable coefficients d * and c * similar to the one discussed in Remark 4.1 is also possible, yielding the two last terms on the right-hand side of (4.38) being replaced by the conservation form ∂ ∂x d * ∂v ∂x + ∂ ∂x (c * v). We now provide some quantitative insights for our model. Extending the test case considered in the previous subsection, we choose d = 0.05 and we enforce that for N = N 0 = 128, we have Q D N0 = Q D,r N0 = 1 and Q C N0 = Q C,r N0 = 2, i.e., each neuron is influenced by its first neighbour on the left and by the two first neighbours on the right. Using

Multi-dimensional dynamics
In this section, we extend the previous one-dimensional treatment, and in particular the material of Section 4.3, to describe the dynamics of a multi-dimensional agglomeration of neurons. We will focus on the main aspects of the analysis, leaving to the reader those details that are straightforward extensions of the one-dimensional results.
We assume that neurons form a periodic lattice contained in B = [0, 1] m , m = 2 or m = 3. Precisely, given any integer n ≥ 2 and setting h = 1/n, each neuron is associated to a multi-index l ∈ {0, · · · , n − 1} m , which identifies its physical position x = hl ∈ B. Thus we have N = n m distinct neurons in   B, which are labelled by indices i ∈ {1, · · · , N } according to some rule; the i−th neurons has position x i = hl i , action potential v i and recovery variable r i . Periodicity means that we replicate the situation at x = hl in any y = h(l + nk) with k ∈ Z m . We adopt again the diffusion model (3.2), with Q given by (4.26). The definition of Q D N and Q C N is as follows: • given a radius R D N := hQ D N with Q D N > 0 (to be determined later on), we set • given a radius R C N := hQ C N with Q C N ≥ Q D N (to be determined later on), and a unit vector ν ∈ R m , we set 2) i.e., Q C N identifies neurons sitting on semi-balls of suitable radii centered at x i ; these semi-balls are obtained by cutting the corresponding balls by the hyperplane containing x i and perpendicular to ν, and retaining the halves oriented in the direction of ν (see Figure 5.1 for a pictorial representation of the sets Q D N and Q C N in two dimensions).

The effect of
Now, it is easily seen that by the form of K D N , the quantity since vectors in K D N can be arranged in couples that are symmetric with respect to each coordinate hyperplane. Thus, where ∆v = m α=1 D 2 αα v is the Laplacian of the function v. We observe for further reference that for any Q > 0, denoting by B(0, Q) the ball of center 0 and radius Q in R m , one has for any given α = 1, · · · , m At this point, we assume that ν = e 1 , the first element of the canonical basis in R m ; this choice is not at all restrictive, but simplifies the following arguments. Indeed, referring to the analogue of (5.3) in which Q D N , K D N resp., are replaced by Q C N , K C N resp., we have On the other hand, But now, We conclude that, going back to the case of an arbitrary ν, The global effect of Q C N Summing up (5.5) and (5.7), we obtain At this point, given two constants d * > 0 and c * ≥ 0, we would like to find This system is similar to (4.31) and we can discuss its solvability as done in Section 4.3. The conclusion is that for N large enough, the solution exists and satisfies whereas the number of neurons that should be connected to a given neuron scales like N 2 m+2 . We summarize our conclusions as follows.
Theorem 5.1. The discrete model (3.2), with Q given by (4.26)-(5.1)-(5.2) in which Q D N and Q C N are the solution of (5.8), tends for N → ∞ to the following continuous model of reaction-convection-diffusion type ∂v ∂t = f (v, r) + d * ∆v +ĉ * · ∇v , ∂r ∂t = g(v, r) , where the convective velocity is given by the vectorĉ * = c * ν.
The well-posedness of this model, as well as its numerical discretization, can be studied by adapting the arguments given in [2] and [11].
An example of a two-dimensional dynamics produced by the model described above is given in Figure 5.2. We fix d = 0.05 as for the one-dimensional models; then, we choose d * and c * in such a way that (5.8) is satisfied for n = 256 by Q D N = √ 2 and Q C N = 2. This gives d * = 3.8147 · 10 −6 and c * = 3.9063 · 10 −4 .
The vector ν is chosen to be e 1 . Figure 5.2 shows the evolution of the action potential in the periodic box B = [0, 1] 2 for n = 256, starting from an initial stimulus v |t=0 = 1 applied to the neurons lying in the circle of radius 1/32 around the center of the box. The stimulus propagates in all directions, but since c * > 0 the speed of propagation is faster in the direction of −ν.

Pseudo-random connections
While the models considered so far are fully deterministic, it is interesting to introduce some form of randomness and monitor its effects. In the simplest form, this can be accomplished by perturbing the model considered above via a (pseudo-)random removal of a fixed percentage of links among neurons. Connections to each neuron are turned-off with uniform distribution in the given percentage, independently of the other neurons; thus, the set Q(i) in (3.2) does depend upon i, in a (pseudo-)random manner.
As an example, we keep the same parameters d = 0.05, Q D N = √ 2, Q C N = 2 and n = 256, as well as the same initial datum as above, and we choose to turn 30% of connections off. In Figure 5.3, the resulting dynamics at the same time instants as in Figure 5.2 is shown. The random effects on the patterns are apparent. The reduction of active connections is reflected by a weaker propagation strength. The excitation front travels leftward only, with a lower speed than in the deterministic case. Furthermore, contours are irregular and, in some realizations not shown here, even disconnected.

Conclusions
We have considered an idealized network, formed by N neurons individually described by the FitzHugh-Nagumo equations and connected by electrical synapses. The limit for N → ∞ of the resulting discrete model has been thoroughly investigated, with the aim of identifying a model for a continuum of neurons having an equivalent behaviour. Two strategies for passing to the limit have been analyzed. A more conventional approach is based on a fixed nearest-neighbour connection topology accompanied by a suitable scaling of the diffusion coefficients. We have devised a new approach, in which the number of connections to any given neuron varies with N according to a precise law, which simultaneously guarantees the non-triviality of the limit and the locality of neuronal interactions. Both approaches yield in the limit a pde-based model, in which the distribution of action potential obeys a nonlinear reaction-convection-diffusion equation; convection accounts for the possible lack of symmetry in the connection topology. Several convergence issues are discussed, both theoretically and numerically.
Based on the present study, in a forthcoming work we will consider more realistic models describing both electrical and chemical synapses. The discrete models here analyzed will be coupled to models of chemical interactions within a population of excitatory/inhibitory neurons, such as those given in [3], eq.(9.6). Again, the focus will be on the limit process leading to coupled continuous models.