## Abstract

Neural continuum networks are an important aspect of the modeling of macroscopic parts of the cortex. Two classes of such networks are considered: voltage and activity based. In both cases, our networks contain an arbitrary number, *n*, of interacting neuron populations. Spatial nonsymmetric connectivity functions represent cortico-cortical, local connections, and external inputs represent nonlocal connections. Sigmoidal nonlinearities model the relationship between (average) membrane potential and activity. Departing from most of the previous work in this area, we do not assume the nonlinearity to be singular, that is, represented by the discontinuous Heaviside function. Another important difference from previous work is that we relax the assumption that the domain of definition where we study these networks is infinite, that is, equal to or . We explicitly consider the biologically more relevant case of a bounded subset Ω of , a better model of a piece of cortex. The time behavior of these networks is described by systems of integro-differential equations. Using methods of functional analysis, we study the existence and uniqueness of a stationary (i.e., time-independent) solution of these equations in the case of a stationary input. These solutions can be seen as ‘persistent’; they are also sometimes called *bumps*. We show that under very mild assumptions on the connectivity functions and because we do not use the Heaviside function for the nonlinearities, such solutions always exist. We also give sufficient conditions on the connectivity functions for the solution to be absolutely stable, that is, independent of the initial state of the network. We then study the sensitivity of the solutions to variations of such parameters as the connectivity functions, the sigmoids, the external inputs, and, last but not least, the shape of the domain of existence *Ω* of the neural continuum networks. These theoretical results are illustrated and corroborated by a large number of numerical experiments in most of the cases 2 ⩽ *n* ⩽ 3, 2 ⩽ q ⩽ 3.

## 1. Introduction

We analyze the ability of neuronal continuum networks to display localized persistent activity, or *bumps*. This type of activity is related, for example, to working memory, which involves holding and processing information on the timescale of seconds. Experiments in primates have shown that there exist neurons in the prefrontal cortex that have high firing rates during the period the animal is “remembering” the spatial location of an event before using the information being remembered (Colby, Duhamel, & Goldberg, 1995; Funahashi, Bruce, & Goldman-Rakic, 1989; Miller, Erickson, & Desimone, 1996). Realistic models for this type of activity have involved spatially extended networks of coupled neural elements or neural masses and the study of spatially localized areas of high activity in these systems. A neuronal continuum network is first built from a “local” description of the dynamics of a number of interacting neuron populations where the spatial structure of the connections is neglected. This local description can be thought of as representing such a structure as a cortical column (Mountcastle, 1957; Mountcastle, 1997; Buxhoeveden & Casanova, 2002). We call it a neural mass (Freeman, 1975). Probably the best-known neural mass model is that of Jansen and Rit (Jansen, Zouridakis, & Brandt, 1993), based on the original work of Lopes da Silva and colleagues (Lopes da Silva, Hoeks, & Zetterberg, 1974, 1976) and of Van Rotterdam and colleagues (van Rotterdam, Lopes da Silva, van den Ende, Viergever, & Hermans, 1982). A complete analysis of the bifurcation diagram of this model can be found in Grimbert and Faugeras (2006). The model has been used to simulate evoked potentials: EEG activities in normal (Jansen & Rit, 1995) and epileptic patients (Wendling, Bartolomei, Bellanger, & Chauvel, 2001, 2000). In a similar vein, David and Friston (2003) have used an extension of this model to simulate a large variety of cerebral rhythms (α, β, γ, δ, and γ) in MEG/EEG simulations. Another important class of such models is the one introduced by Wilson and Cowan (1973); Hoppenstaedt & Izhikevich (1997).

These local descriptions are then assembled spatially to form the neuronal continuum network. This continuum network is meant to represent a macroscopic part of the neocortex, for example, a visual area such as V1. The spatial connections are models of cortico-cortical connections. Other, nonlocal connections with, for example, such visual areas as the lateral geniculate nucleus, or V2, are also considered. Other researchers have used several interconnected neural masses to simulate epileptogenic zones (Wendling et al., 2000, 2001; Lopes da Silva et al., 2003) or to study the connectivity between cortical areas (David, Cosmelli, & Friston, 2004). In this letter we consider a continuum of neural masses.

## 2. The Models

We briefly discuss local and spatial models.

### 2.1. The Local Models.

We consider *n* interacting populations of neurons such as those shown in Figure 1. The figure is inspired by the work of Alex Thomson (Thomson & Bannister, 2003) and Wolfgang Maass (Haeusler & Maass, 2007). It shows six interacting populations of neurons. Black indicates excitation and gray inhibition. The thickness of the arrows pertains to the strength of the interaction. In this example, the six populations are located in layers 2/3, 4, and 5 of the neocortex. Each population being described by its state (defined below), we derive from first principles the equations that describe the variations over time of the states of the different interacting populations.

Our derivation follows closely that of Ermentrout (1998). We consider that each neural population *i* is described by its average membrane potential *V _{i}*(

*t*) or its average instantaneous firing rate ν

_{i}(

*t*), the relation between the two quantities being of the form ν

_{i}(

*t*) =

*S*(

_{i}*V*(

_{i}*t*)) (Gerstner and Kistler, 2002; Dayan & Abbott, 2001), where

*S*is sigmoidal and smooth. The functions

_{i}*S*,

_{i}*i*= 1, …,

*n*, satisfy the following properties introduced in the following definition:

*For all i = 1, …, n, |S_{i}| ⩽ S_{im} (boundedness). We note S_{m} = max_{i}S_{im}. For all i = 1, …, n, the derivative S′_{i} of S_{i} is positive and bounded by S′_{im}>0 (boundedness of the derivatives). We note DS_{m} = max_{i}S′_{im} and DS_{m} the diagonal matrix diag(S′_{im}).*

*S*is given in equation 2.1: This function is symmetric with respect to the “threshold” potential θ

_{i}_{i}and varies between 0 and 1. The positive parameter

*s*controls the slope of the

_{i}*i*th sigmoid at

*v*= θ

_{i}.

*s*= 0.5, 1, 10. We have

*S*= 1 and

_{im}*S*′

_{im}=

*s*. When

*s*→ ∞,

*S*converges to the Heaviside function

*Y*defined by

Neurons in population *j* are connected to neurons in population *i*. A single action potential from neurons in population *j* is seen as a postsynaptic potential *PSP _{ij}*(

*t*−

*s*) by neurons in population

*i*, where

*s*is the time of the spike hitting the terminal and

*t*the time after the spike. We neglect the delays due to the distance traveled down the axon by the spikes.

*i*due to action potentials of population

*j*is where the sum is taken over the arrival times of the spikes produced by the neurons in population

*j*. The number of spikes arriving between

*t*and

*t*+

*dt*is ν

_{j}(

*t*)

*dt*. Therefore, we have or, equivalently, The

*PSP*can depend on several variables in order to account for adaptation and learning, among other examples.

_{ij}There are two main simplifying assumptions that appear in the literature (Ermentrout, 1998) and produce two different models.

#### 2.1.1. The Voltage-Based Model.

*w*>0, the population

_{ij}*j*excites population

*i*, whereas it inhibits it when

*w*< 0.

_{ij}*A*in the weights

_{i}*w*and added an external current

_{ij}*I*

_{ext}(

*t*) to model the nonlocal connections of population

*i*. We introduce the

*n*×

*n*matrixes

**W**such that

*W*=

_{ij}*w*/τ

_{ij}_{i}and the function

**S**, such that

**S**(

**x**) is the vector of coordinates

*S*(

_{i}*x*), if

_{i}*x*= (

*x*

_{1}, …,

*x*). We rewrite equation 2.4 in vector form and obtain the following system of

_{n}*n*ordinary differential equations, where

**L**is the diagonal matrix

**L**=diag(1/τ

_{i}).

In terms of units, the left- and right-hand sides of this equations are in units of, say, mV × ms^{−1}. Therefore, **I**_{ext}, despite its name, is not a current. Note that since **S**(**V**) is an activity, its unit is ms^{−1}, and hence **W** is in mV.

#### 2.1.2. The Activity-Based Model.

*PSP*(

_{i}*t*) satisfies the differential equation 2.3 and define the time-averaged firing rate to be A similar derivation yields the following set of

*n*ordinary differential equations: We include the τ

_{i}s in the sigmoids

*S*and rewrite this in vector form: The units are ms

_{i}^{−2}for both sides of the equation.

**W**is expressed in mV × ms, and

**I**

_{ext}is in mV.

### 2.2. The Continuum Models.

We now combine these local models to form a continuum of neural masses, for example, in the case of a model of a significant part Ω of the cortex. We consider a subset Ω of , *q* = 1, 2, 3, which we assume to be connected and compact (i.e., closed and bounded). We note |Ω| its Lebesgue measure (length, area, volume). This encompasses several cases of interest.

When *q* = 1, we deal with one-dimensional sets of neural masses. Although this appears to be of limited biological interest, this is one of the most widely studied cases because of its relative mathematical simplicity and because of the insights one can gain of the more realistic situations.

When *q* = 2, we discuss properties of two-dimensional sets of neural masses. This is perhaps more interesting from a biological point of view since Ω can be viewed as a piece of cortex where the third dimension, its thickness, is neglected. This case has received by far less attention than the previous one, probably because of the increased mathematical difficulty. Note that we could also take into account the curvature of the cortical sheet at the cost of an increase in mathematical difficulty. This is outside the scope of this letter.

Finally *q* = 3 allows us to discuss properties of volumes of neural masses, such as cortical sheets where their thickness is taken into account (Kandel, Schwartz, & Jessel, 2000; Chalupa & Werner, 2004).

The theoretical results presented in this letter are independent of the value of *q*.

**V**(

**r**,

*t*) (resp.,

**A**(

**r**,

*t*)) the

*n*-dimensional state vector at the point

**r**of the continuum and at time

*t*. We introduce the

*n*×

*n*matrix function

**W**(

**r**,

**r**′), which describes how the neural mass at point

**r**′ influences that at point

**r**at time

*t*. We call

**W**the connectivity matrix function. In particular,

**W**(

**r**,

**r**) =

**W**, the matrix that appears in equations 2.5 and 2.6. More precisely,

*W*(

_{ij}**r**,

**r**′) describes how population

*j*at point

**r**′ influences population

*i*at point

**r**at time

*t*. Equation 2.5 can now be extended to and equation 2.6 to

It is important to discuss again the units of the quantities involved in these equations. For equation 2.7, as for equation 2.4, the unit is mV × ms^{−1} for both sides. Because of the spatial integration, **W** is in mV × ms^{−1} × mm^{−q}, and *q* is the dimension of the continuum. To obtain a dimensionless equation, we normalize (i.e., divide both sides of the equation) by the Frobenius norm ‖**W**‖_{F} of the connectivity matrix function **W** (see section A.1 for a definition). Equivalently, we assume that ‖**W**‖_{F} = 1.

We have given elsewhere (Faugeras, Grimbert, & Slotine, 2007)—but see proposition 2 below for completeness—sufficient conditions on **W** and **I**_{ext} for equations 2.7 and 2.8 to be well defined and studied the existence and stability of their solutions for general external currents. In this letter, we analyze in detail the case of stationary external currents (independent of the time variable) and investigate the existence and stability of the corresponding stationary solutions of 2.7 and 2.8.

A significant amount of work has been devoted to this or closely related problems, starting perhaps with the pioneering work of Wilson and Cowan (1973). A fairly recent review of this work, and much more, can be found in Coombes (2005). Amari (1977) investigated the problem in the case *n* = *q* = 1 when the sigmoid function is approximated by a Heaviside function and the connectivity function has a Mexican hat shape. He proved the existence of stable bumps in this case. His work has been extended to different firing rate and connectivity functions (Gutkin, Ermentrout, & O'Sullivan, 2000; Laing, Troy, Gutkin, & Ermentrout, 2002; Laing & Troy, 2003; Rubin & Troy, 2004; Guo & Chow, 2005a, 2005b).

The case *n* = 1, *q* = 2 has been considered by several authors, including Pinto and Ermentrout (2001a, 2001b) for general firing rate functions and gaussian-like connectivity functions (Blomquist, Wyller, & Einevoll, 2005) when the firing rate functions are approximated by Heaviside functions, and more recently for translation-invariant kernels including some discussion of sigmoids in Owen, Laing, and Coombes (2007).

Extending these analyses to a two- or three-dimensional continuum is difficult because of the increase in the degrees of freedom in the choice of the connectivity function. The case *n* = 2, *q* = 1 has been studied in Werner and Richter (2001) and Bressloff (2005) when the firing rate functions are approximated by Heaviside functions and the connectivity function is circularly symmetric, while the case *n* = 2, *q* = 2 is mentioned as difficult in Doubrovinski (2005).

In all of these contributions, the proof of the existence of a bump solution is based on Amari's original argument (1977), which works only when *q* = 1 and the firing rate function is approximated by a Heaviside function. Kubota and Aihara (2005); Kishimoto and Amari (1979) consider sigmoids in the case *q* = 1 for translation-invariant symmetric connectivity functions. Solutions are usually constructed using a variant of the method of the singular perturbation construction (Pinto & Ermentrout, 2001b). Sufficient conditions for their stability are obtained by a linear stability analysis, which in general requires the use of Heaviside functions instead of sigmoids.

The approach that we describe in this letter is a significant departure from previous ones. By using simple ideas of functional analysis, we are able to:

- •
Prove the existence and uniqueness of a stationary solution to equations 2.7 and 2.8 for any dimensions

*n*and*q*, arbitrary connectivity functions, and general firing rate functions. - •
Obtain very simple conditions for the absolute stability of the solution in terms of the spectrum of the differential of the nonlinear operator that appears on the right-hand side of equations 2.7 and 2.8.

- •
Construct a numerical approximation as accurately as needed of the solution, when it exists, for any stationary input.

- •
Characterize the sensitivity of the solutions to variations of the parameters, including the shape of the domain Ω.

To be complete, let us point out that equations of the type 2.7 and 2.8 have been studied in pure mathematics (see e.g., Hazewinkel, 2001). They are of the Hammerstein type (Hammerstein, 1930; Tricomi, 1985). This type of equation has received some recent attention (see Appell & Chen, 2006), and progress has been made toward a better understanding of their solutions. Our contributions are the articulation of the models of networks of neural masses with this type of equation, the characterization of persistent activity in these networks as fixed points of Hammerstein equations, the proof of the existence of solutions, the characterization of their stability, and the analysis of their sensitivity to variations of the parameters involved in the equations.

The rest of the letter is organized as follows. In section 3, we prove that under some mild hypotheses about the connectivity matrix, because Ω is bounded, there always exist stationary solutions to the neural field equations, 2.7 and 2.8, for stationary inputs. In section 4, we provide sufficient conditions for these stationary solutions to be stable. Section 5 provides numerical support for the theorems of the previous two sections: because we know that there exist stable solutions, we can derive a robust and accurate numerical scheme to compute them for various values of the parameters. In section 6, we probe the sensitivity of the solutions to variations of these parameters. This is a first step in the direction of the study of their bifurcations, which will be the topic of a forthcoming paper. In section 7, we summarize our results and open some perspectives.

## 3. Existence of Stationary Solutions

In this section we deal with the problem of the existence of stationary solutions to equations 2.7 and 2.8 for a given stationary external current **I**_{ext}.

**L**

^{2}

_{n}(Ω) of square integrable functions from Ω to . This is a Hilbert, hence a Banach, space for the usual inner product, where is the complex conjuguate of the vector

**V**. This inner product induces the norm (see section A.1). is the state space. Another important space is

**L**

^{2}

_{n×n}(Ω × Ω), the space of square integrable

*n*×

*n*matrices (see section A.1 for a precise definition). We assume that the connectivity matrix functions

**W**(·, ·) are in this space (see propositions 1 and 2 below).

We first recall some results on the existence of a solution to equations 2.7 and 2.8 that will be used in the sequel.

**X**(

**r**,

*t*) is a mapping , and equations 2.7 and 2.8 are formally recast as an initial value problem: where

**X**

_{0}is an element of and the function

*f*from is defined by the right-hand side of equation 2.7, in which case we call it

*f*, or equation 2.8, in which case we call it

_{v}*f*. In other words, equations 2.7 and 2.8 become differential equations defined on the Hilbert space .

_{a}We need the following two propositions, which we quote without proof (Faugeras et al., 2007):

*If the following two hypotheses are satisfied,*

*The connectivity function*is in*W**L*^{2}_{n×n}(Ω × Ω) (see section A.1),*At each time instant**t*∈*J*the external currentis in*I**, the set of continuous functions from**J*to*,*

*then the mappings*

*f*and_{v}*f*are from_{a}*to*

*, continuous, and Lipschitz continuous with respect to their second argument, uniformly with respect to the first.*

*If the following two hypotheses are satisfied,*

*The connectivity function*is in*W**L*^{2}_{n×n}(Ω × Ω),*The external current**I*_{ext}is in*, the set of continuous functions from*,*J*to

*then for any function*

*X*_{0}in*, there is a unique solution*

**, defined on***X**(and not only on*.

*J*) and continuously differentiable, of the abstract initial value problem, equation 3.1, for*f*=*f*and_{v}*f*=*f*_{a}This proposition says that given the two hypotheses and the initial condition, there exists a unique solution to equation 2.7 or 2.8 and that this solution is in , the set of continuously differentiable functions from to .

We now turn our attention to a special type of solution of equations 2.7 and 2.8, corresponding to stationary external currents. We call these solutions, when they exist, *stationary solutions*. The currents **I**_{ext} are simply in .

We now recall:

*A continuous mapping**(linear or nonlinear) is called compact provided that for each bounded subset**of**, the set**is relatively compact; that is, its closure is compact.*

*If W ∈ L^{2}_{n×n}(Ω × Ω), g^{L}_{v} and g^{L}_{a} are compact operators of*.

We know from proposition 1 that *g ^{L}_{v}* is continuous and prove that for each bounded sequence {

**V**

_{n}}

^{∞}

_{n=1}of , there exists a subsequence such that is convergent in .

Because of definition 1 of **S**, the sequence {**A**_{n} = **S**(**V**_{n})}^{∞}_{n=1} is bounded in by . We prove that there exists a subsequence such that converges in .

Since is separable, its unit ball is weakly compact, and because {**A**_{n}}^{∞}_{n=1} is bounded, there exists a a subsequence of {**A**_{n}}^{∞}_{n=1} that converges weakly in toward **A**. Because of Fubini's theorem, for almost all **r** ∈ Ω (noted almost surely), the function **r**′ → **W**(**r**, **r**′) is in . Therefore, almost surely,

Since , **A** is also bounded by *C* in . It is easy to show that , and we can apply Lebesgue's dominated convergence theorem to the sequence and conclude that , that is, is convergent in .

A small variation of the proof shows that *g ^{L}_{a}* is compact.

From proposition 3 follows:

*Under the hypotheses of proposition 3, if**, f^{L}_{v} and f^{L}_{a}, are compact operators of*.

The operators **X** → **I**^{L}_{ext} and **X** → **I**_{ext} are clearly compact under the hypothesis ; therefore, *f ^{L}_{v}* is the sum of two compact operators, hence compact. For the same reason,

*g*+

_{a}**I**

_{ext}is also compact, and so is

*f*=

^{L}_{a}**S**

^{L}(

*g*+

_{a}**I**

_{ext}) because

**S**

^{L}is smooth and bounded.

We can now prove:

Note that in the proofs of proposition 4 and theorem 1, the assumption that Ω is bounded (i.e., that |Ω| is finite) is essential.

## 4. Stability of the Stationary Solutions

### 4.1. The Voltage-Based Model.

**W**

_{cm}(

**r**,

**r**′) by

**W**

_{cm}=

**W**

*D*

**S**

_{m}, where

*D*

**S**

_{m}is defined in definition 1. We also define the corresponding linear operator , which is compact according to proposition 3. Its adjoint, noted

*h**

_{m}is defined

^{1}by and is also compact. Hence the symmetric part , the sum of two compact operators, is also compact. Furthermore, we have 〈

**V**,

*h*(

_{m}**V**)〉 = 〈

**V**,

*h*(

^{s}_{m}**V**)〉, as can be easily verified. It is also self-adjoint since, clearly,

*h*=

^{s}_{m}*h**

^{s}_{m}.

We recall the following property of the spectrum of a compact self-adjoint operator in a Hilbert space (see, e.g., Dieudonné 1960).

*The spectrum of a compact, self-adjoint operator of a Hilbert space is countable and real. Each nonzero spectral value is an eigenvalue, and the dimension of the corresponding eigenspace is finite.*

We have the following:

*A sufficient condition for the stability of a stationary solution to equation 2.7 is that all the eigenvalues of the linear compact, self-adjoint operator*

*h*^{L,s}_{m}be less than 1, where*h*^{L,s}_{m}is defined by*where*:

*h*^{L,s}_{m}is the symmetric part of the linear compact operatorThe proof of this theorem is a generalization to the continuum case of a result obtained by Matsuoka (1992).

*D*

**S**

_{m})

^{−1}

**S**and rewrite equation 2.7 for a homogeneous input

**I**

_{ext}as follows:

**U**be a stationary solution of equation 2.7. Also let

**V**be the unique solution of the same equation with initially some condition (see proposition 2). We introduce the new function

**X**=

**V**−

**U**, which satisfies where the vector

**(**

*Θ***X**) is given by . Consider now the functional We note that This is because (Taylor expansion with integral remainder) and by construction of the vector . This implies that the functional Δ(

**X**) is strictly positive for all and Δ(0) = 0. It also implies—and this is used in the sequel—that

*z*Θ

_{i}(

*z*) ⩾ Θ

_{i}(

*z*)

^{2}.

Note that we have the following:

*If the condition of theorem 2 is satisfied, the homogeneous solution of equation**2.7**is unique.*

Indeed, the result of theorem 2 is independent of the particular stationary solution **U** that is chosen in the proof.

### 4.2. The Activity-Based Model.

**W**

_{mc}=

*D*

**S**

_{m}

**W**and the linear compact operator

*k*from to :

_{m}*A sufficient condition for the stability of a solution to equation*

*2.8*

*is that all the eigenvalues of the linear compact operator*

*k*be of a magnitude less than 1, where^{L}_{m}*k*is defined by^{L}_{m}**U**be a stationary solution of equation 2.8 for a stationary external current

**I**

_{ext}(

**r**). As in the proof of theorem 2, we introduce the new function

**X**=

**V**−

**U**, where

**V**is the unique solution of the same equation with initial conditions , an element of . We have Using a first-order Taylor expansion with an integral remainder, this equation can be rewritten as Consider now the functional . Its time derivative is We replace

**X**

_{t}(

**r**,

*t*)) by its value in this expression to obtain where the nonlinear operator σ

_{m}is defined by a diagonal matrix whose diagonal elements are between 0 and 1. We rewrite in a slightly different manner, introducing the operator

*k*: From the Cauchy-Schwarz inequality and the property of σ

^{L}_{m}_{m}(

**X**), we obtain A sufficient condition for to be negative is therefore that for all

**Y**.

## 5. Numerical Experiments

In this section and the next, we investigate the question of effectively (i.e., numerically) computing stationary solutions of equation 2.7, which is equivalent to computing solutions of equation 3.2. Similar results are obtained for equation 2.8.

In all of our numerical experiments, we assume the sigmoidal functions *S _{i}*,

*i*= 1, …,

*n*introduced in definition 1 to be of the form 2.1.

### 5.1. Algorithm.

Note that this method allows to computing the solution of equation 3.2 only when it admits a unique solution and *f* is contracting. However, it could admit more than one solution (recall it always has a solution; see theorem 1) or *f* could be noncontracting. Another method has to be found in these cases.

In our case, , where Ω is an open-bounded set of and *M* = *f _{v}*. According to lemmas 2 and 1, if

*DS*‖

_{m}**W**‖

_{F}< 1,

*f*is contracting.

_{v}**V**

_{p},

*p*⩾ 0 is approximated by a piecewise constant function

**V**

_{p,h}, where

*h*is the step of the approximation, defined by a finite number of points . In order to avoid difficulties because

**V**

_{p,h}∈

**L**

^{2}

_{n}(Ω), hence defined almost everywhere, we assume that

**W**and

**I**are smooth. It is not a restriction because every function of

**L**

^{2}

_{n}(Ω) can be approximated by a smooth function. As the bump solution is smooth, as soon as

**W**,

**I**are smooth, we can use the multidimensional gaussian quadrature formula (Press, Flannery, Teukolsky, & Vetterling, 1988; Stoer & Bulirsch, 1972) with

*N*points (in the examples below, usually

*N*= 20) on each axis. In order to interpolate the values of the bump from a finite (and small) number of its values

**V**

_{n}(

**r**

_{h,j,Gauss}), we use Nyström's method (Hazewinkel, 2001), as follows: where the

*g*s are the weights of the gaussian quadrature method and the points

_{j}**r**

_{p,j,Gauss}are chosen according to the gauss quadrature formula. It is to be noted that the choice of a particular quadrature formula can make a huge difference in accuracy and computing time (see section A.2).

The following theorem ensures that lim_{h→0}**V**^{f}_{h} = **V**^{f}:

*Assume that lim _{h→0}W_{h} = W in L^{2}_{n×n}(Ω × Ω) then*

*with*

*a*= ‖_{h}**−***W**W*_{h}‖_{F}.The proof is an adaptation of (Krasnosel'skii, Vainikko, Zabreiko, & Stetsenko, 1972, theorem 19.5).

### 5.2. Examples of Bumps.

We show four examples of the application of the previous numerical method to the computation of bumps for various values of *n* and *q*.

*n*populations (

**V**= [

*V*

_{1}, …,

*V*]

_{n}^{T},

**W**∈

**L**

^{2}

_{n×n}(Ω × Ω)), some excitatory and some inhibitory. Ω = [ − 1, 1]

^{q}. We characterize in section 6.2 how the shape of Ω influences that of the bumps. The connectivity matrix is of the form with a

*q*×

*q*symmetric positive definite matrix. The weights α

_{ij},

*i*,

*j*= 1, …,

*n*form an element

**of , and is an element of . The weights**

*α***are chosen so that**

*α**DS*‖

_{m}**W**‖

_{F}< 1. The sign of α

_{ij},

*i*≠

*j*, indicates whether population

*j*excites or inhibits population

*i*. The bumps are computed using the algorithm described in the previous section.

#### 5.2.1. First Example: *n* = 2, *q* = 2, Constant Current.

#### 5.2.2. Second Example: *n* = 2, *q* = 2, Nonconstant Current.

**I**is still equal to 0 for its second coordinate and is not constant but equal to its previous value, −0.3, to which we have added a circularly symmetric 2D gaussian centered at the point of coordinates (0.5, 0, 5) of the square Ω with standard deviation 0.18 and maximum value 0.2. It is interesting to see how the shape of the previous bump is perturbed. The matrix

**is the same as in the first example. The matrix**

*α***T**is equal to corresponding to a spatially broader interaction for the first population and narrower for the second.

#### 5.2.3. Third Example: *n* = 3, *q* = 2, Constant Current.

Figure 5 shows an example of a bump for three neural populations—two excitatory and one inhibitory—in two dimensions. We use the following values of the parameters:

#### 5.2.4. Fourth Example: *n* = 2, *q* = 3, Constant Current.

_{3}is the 3 × 3 identity matrix.

## 6. Sensitivity of the Bump to Variations of the Parameters

In this section, we characterize how the solutions of equation 3.2 vary with the parameters that appear in the equation. These parameters are of two types: first, we have a finite number of real parameters such as the external currents, the weights in the connectivity matrix **W**, or the parameters of the sigmoids, and, second, the shape of the domain Ω, a potentially infinite-dimensional parameter.

We focus on the voltage-based model; the analysis in the activity-based case is very similar. We start with a set of general considerations in the finite-dimensional case, which we then apply to the various cases. We then tackle the more difficult case of the dependency with respect to the shape of Ω.

As *f _{v}* is a smooth function of the parameters (

**I**,

**,**

*α***S**, …), one can show (by extending Banach's theorem) that the fixed-point

**V**

^{f}inherits the smoothness of

*f*.

_{v}### 6.1. The Finite-Dimensional Parameters.

Note that ∂_{λ}**S**(**V**^{f}) is the partial derivative of the vector **S** with respect to the scalar parameter λ evaluated at **V** = **V**^{f}.

#### 6.1.1. Sensitivity of the Bump to the Exterior Current.

*I*

_{1}, we find: This inequality is to be understood component by component. It predicts the influence of

*I*

_{1}on

**V**

^{f}. For example, with the parameters

**and**

*α***T**used in Figure 3 but with an external current equal to 0, we obtain the bump shown in Figure 7 (top) with the derivatives shown at the bottom of the same figure. We also show in Figure 8 of

*V*

_{1}and

*V*

_{2}along the diagonal and the

*x*-axis for different values of

*I*

_{1}close to 0. The reader can verify that the values increase with

*I*

_{1}, as predicted.

#### 6.1.2. Sensitivity of the Bump to the Weights *α*.

*α*

The fixed point is a decreasing function of the inhibitory parameter *b* (see Figure 9).

The other cases are similar.

#### 6.1.3. Sensitivity of the Bump to the Thresholds.

_{i},

*i*= 1, 2 we have from definition 1 of

**S**and with the notations of proposition 6, where

**e**

_{1}= [1, 0]

^{T},

**e**

_{2}= [0, 1]

^{T}. We show in Figure 10 some cross sections of the bump

**V**

^{f}obtained for the same values of the parameters as in Figure 3 and three values of the threshold vector.

#### 6.1.4. Sensitivity of the Bump to the Slope of the Sigmoid.

*s*,

_{i}*i*= 1, 2, we have from definition 1 of

**S**and with the notations of proposition 6, where the matrix

**S**is given by Figure 11 shows the two coordinates ∂

_{s}

*V*

^{f}_{1}and ∂

_{s}

*V*

^{f}_{2}for

*s*

_{1}=

*s*

_{2}=

*s*of the derivative of the bump

**V**

^{f}at

*s*= 1 obtained for the same values of the other parameters as in Figure 3, except the intensity which is equal to 0.

### 6.2. Sensitivity of the Bump to Variations of the Shape of the Domain *Ω*.

*Ω*

We expect the bump to be somewhat dependent on the shape of Ω. It would nonetheless be desirable that this dependency would not be too strong for the modeling described in this letter to have some biological relevance. Indeed, if the bumps are metaphors of persistent states, we expect them to be relatively robust to the actual shape of the cortical part being modeled. For example, if we take Ω to be a representation of the primary visual cortex *V*1 whose shape varies from individual to individual, it would come as a surprise if the shape of a bump induced by the same spatially constant stimulus were drastically different.

Technically, in order to study the dependence of **V**^{f} with respect to Ω, we need to assume that Ω is smooth; its border ∂Ω is a smooth curve (*q* = 2) or surface (*q* = 3) unlike the previous examples, where Ω was the square [ − 1, 1]^{2}. But a difficulty arises from the fact that the set of regular domains is not a vector space—hence, the derivative of a function (the bump) with respect to a domain has to be defined with some care. The necessary background is found in section A.3.

*Let us assume that Ω is a smooth, bounded domain of*

*. If*

**is in***W**W*^{1,2}_{n×n}(Ω × Ω),*I*_{ext}is in*W*^{1,2}_{n}(Ω) (see section A.1 for a definition) the material derivative (see section A.3 for a definition)*V*^{f}_{m}(**, Ω) of the bump***r**V*^{f}satisfies the following equation:*where*

*D*,_{i}*i*= 1, 2 indicates the derivative with respect to the*i*th variable and 〈*J*′(Ω),**〉 is the Gâteaux derivative of***X**J*(Ω) with respect to the vector field**:***X**where*

**is defined in the proof below. We have***X**where*

*d***is the surface element on the smooth boundary ∂Ω of Ω, and***a***is its unit inward normal vector.***n**T*: We define In the first integral, we make the change of variable

**r**′ →

**r**′ + τ

**X**and obtain We have Hence, for τ sufficiently small det

*J*

_{τ}>0. Moreover: and where

*D*,

_{i}*i*= 1, 2 indicates the derivative with respect to the

*i*th argument. Thus, we have Because and ∫

_{Ω}

**W**(

**r**,

**r**′, Ω)

**S**(

**V**

^{f}(

**r**′, Ω))

*d*

**r**′ =

**V**

^{f}(

**r**, Ω) −

**I**

_{ext}(

**r**), we obtain equation 6.2. The value of 〈

*J*′(Ω),

**X**〉 is obtained from corollary 2.

#### 6.2.1. Numerical Application for the Domain Derivatives.

**V**

^{f}for Ω equal to the unit disc

*D*(0, 1) and in Figure 13 the one for Ω equal to the ellipse

^{3}of equation . The values of the weight parameters

**are the same as in Figure 3 and**

*α***I**= [0, 0]

^{T}. The matrix

**T**is equal to Note that because the diagonal elements are not equal for

*T*

_{11},

*T*

_{12}, and

*T*

_{13},

**W**is not circularly symmetric, and so is the bump in Figure 12 despite the fact that Ω is circularly symmetric.

**X**corresponding to the transformation

*T*(1) transforms the disc

*D*(0, 1) into the ellipse ,

**X**(

**r**) = [(

*a*− 1)

*r*

_{1}, 0]

^{T}.

**T**

_{ij}are symmetric): Thus, we can use a simple modification of the algorithm that computes

**W**·

**S**(

**V**) to obtain the previous expression.

*J*(Ω) and 〈*J*′(Ω), **X**〉 are computed with a gauss quadrature formula. For a circle in polar coordinates, **N**(**r**′) = **r**′.

Let us be a bit more precise. In the case shown in Figure 12, we choose **I** = 0. Using Banach's theorem, we compute **V**^{f}_{Gauss} for *N* = 30 and use Nyström's interpolation to compute **V**^{f}_{Nys} for *n* = 100 (for example) points on each axis.

**V**

^{f}

_{Gauss}, we compute

**V**

^{f}

_{m,Gauss}for

*N*points. But the equation for

**V**

^{f}

_{m}reads Having computed a Nyström interpolation of

*n*points for

**V**

^{f}

_{m}=

**W**.

*D*

**S**(

**V**

^{f}).

**V**

^{f}

_{m}+ 〈

*J*′(Ω),

**X**〉, we again use a Nyström interpolation with the last equation to compute

**V**

^{f}

_{m,Nystrom}for

*n*points on each axis.

We used this numerical method in every previous example related to the computation of derivatives.

## 7. Conclusion and Perspectives

We have studied two classes (voltage and activity based) of neural continuum networks in the context of modeling macroscopic parts of the cortex. In both cases, we have assumed an arbitrary number of interacting neuron populations, either excitatory or inhibitory. These populations are spatially related by nonsymmetric connectivity functions representing cortico-cortical, local connections. External inputs are also present in our models to represent nonlocal connections, for example, with other cortical areas. The relationship between (average) membrane potential and activity is described by nondegenerate sigmoidal nonlinearities, that is, not by Heaviside functions, which have often been considered instead in the literature because of their (apparent) simplicity.

The resulting nonlinear integro-differential equations are of the Hammerstein type (Hammerstein, 1930) and generalize those proposed by Wilson and Cowan (1973).

Departing from most of the previous work in this area, we relax the usual assumption that the domain of definition where we study these networks is infinite, that is, equal to or , and we explicitly consider the biologically much more relevant case of a bounded subset Ω of , obviously a better model of a piece of cortex. The importance of this has been emphasized by Nunez (2005) for a number of years, and we fully agree with him on this.

Using methods of functional analysis, we have studied the existence and uniqueness of a stationary (i.e., time-independent) solution of these equations in the case of a stationary input. These solutions are often referred to as persistent states, or bumps, in the literature.

We have proved that under very mild assumptions on the connectivity functions, such solutions always exist (this is due in part to the fact that we do not use Heaviside functions and mostly to the fact that we consider the cortex as a bounded set).

We have provided sufficient conditions on the connectivity functions for the solution to be absolutely stable, that is, independent of the initial state of the network. These conditions can be expressed in terms of the spectra of some functional operators, which we prove to be compact, that arise naturally from the equations describing the network activity.

We have also studied the sensitivity of the solutions to variations of such parameters as the connectivity functions, the sigmoids, the external inputs, and the shape of the domain of definition of the neural continuum networks. This last analysis is more involved than the others because of the infinite-dimensional nature of the shape parameter. An analysis of the bifurcations of the solutions when the parameters vary over large ranges requires the use of techniques of bifurcation analysis for infinite-dimensional systems and is out of the scope of this letter.

We believe, and we hope by now to have convinced the reader, that the functional analysis framework that we have used in this letter is the right one to try to answer some of the mathematical questions raised by models of connected networks of nonlinear neurons. We also believe that some of these also begin to answer biological questions since these networks models, despite the admitedly immense simplifications they are built from, are nonetheless metaphors of real neural assemblies.

## Appendix A: Notations and Background Material

### A.1. Matrix Norms and Spaces of Functions.

*n*×

*n*real matrixes. We consider the Frobenius norm on , and consider the space

**L**

^{2}

_{n×n}(Ω × Ω) of the functions from Ω × Ω to whose Frobenius norm is in

*L*

^{2}(Ω × Ω). If

**W**∈

**L**

^{2}

_{n×n}(Ω × Ω), we note ‖

**W**‖

^{2}

_{F}= ∫

_{Ω×Ω}‖

**W**(

**r**,

**r**′)‖

^{2}

_{F}

*d*

**r**

*d*

**r**′. Note that this implies that each element

*W*,

_{ij}*i*,

*j*= 1, …,

*n*in in

*L*

^{2}(Ω × Ω). We note the set

**L**

^{2}

_{n}(Ω) of square-integrable mappings from Ω to and the corresponding norm. We have the following:

*W*is in

_{ij}*L*

^{2}(Ω × Ω),

*W*(

_{ij}**r**, .) is in

*L*

^{2}(Ω) for almost all

**r**, thanks to Fubini's theorem. So

*W*(

_{ij}**r**, .)

*x*(.) is integrable for almost all

_{j}**r**from what we deduce that

**Y**is well defined for almost all

**r**. Next, we have and (Cauchy-Schwarz) from which it follows that (Cauchy-Schwarz again, discrete version): from which it follows that

**Y**is in

**L**

^{2}

_{n}(Ω) (thanks again to Fubini's theorem) and

We also use the following:

We have , where |Ω| is the Lebesgue measure of Ω (its area). Similarly, .

In theorem 6, we use the Sobolev spaces **W**^{1,2}_{n}(Ω) and **W**^{1,2}_{n×n}(Ω × Ω). **W**^{1,2}_{n}(Ω) is the set of functions such that each component *X _{i}*,

*i*= 1, …,

*n*is in

*W*

^{1,2}(Ω), the set of functions of

*L*

^{2}(Ω) whose first-order derivatives exist in the weak sense and are also in

*L*

^{2}(Ω) (see Evans, 1998). Similarly

**W**

^{1,2}

_{n×n}(Ω × Ω) is the set of functions such that each component

*X*,

_{ij}*i*,

*j*= 1, …,

*n*is in

*W*

^{1,2}(Ω).

### A.2. Choice of the Quadrature Method.

We emphasize the importance of the choice of a specific quadrature formula using the following example: where we compare a 0th-order finite elements methods with Gauss's method (the parameters of the Gauss quadrature formula are computed with a precision of 10^{−16} using Newton's method):

Method . | Value . |
---|---|

Exact | 2.350 402 387 287 603 … |

0th order (N = 1000) | 2.351 945 … |

Finite element | |

Gauss (N = 5) | 2.350 402 386 46 … |

Method . | Value . |
---|---|

Exact | 2.350 402 387 287 603 … |

0th order (N = 1000) | 2.351 945 … |

Finite element | |

Gauss (N = 5) | 2.350 402 386 46 … |

The Gauss method is far more powerful and allows us to compute bumps in 3D for an arbitrary number of populations.

### A.3. Shape Derivatives.

*T*and

*T*

^{−1}one to one. In detail, where the functional space is the set of mappings such that they and their first-order derivatives are in . In detail,

Given a shape function , for , let us define . The key point is that since is a Banach space, we can define the notion of a derivative with respect to the domain Ω as:

*F**is Gâteaux differentiable with respect to Ω if and only if**is Gâteaux differentiable with respect to T.*

In order to compute Gâteaux derivatives with respect to *T*, we introduce a family of deformations (*T*(τ))_{τ⩾0} such that for τ ⩾ 0, *T*(0) = *Id*, and . From a practical point of view, there are many ways to construct such a family, the most famous one being the Hadamard deformation (Hadamard, 1968), which goes as follows.

We also introduce:

and

The following theorem whose proof can be found, in Delfour and Zolésio (2001) and Sokolowski and Zolésio (1992) relates the Gâteaux derivative and the shape derivative:

The following corollary is used in the proof of theorem 6:

## Notes

^{1}

By definition, 〈**V**_{1}, *h _{m}*(

**V**

_{2})〉 = 〈

*h**

_{m}(

**V**

_{1}),

**V**

_{2}〉, for all elements

**V**

_{1},

**V**

_{2}of .

^{2}

**W** · *D***S**(**V**^{f}) is the Frechet derivative of the operator *f _{v}* at the point

**V**

^{f}of .

^{3}

represents the ellipse lying along the first axis of coordinates with semimajor axis *a* and semiminor axis *b*.