Analysis of a hyperbolic geometric model for visual texture perception

We study the neural field equations introduced by Chossat and Faugeras to model the representation and the processing of image edges and textures in the hypercolumns of the cortical area V1. The key entity, the structure tensor, intrinsically lives in a non-Euclidean, in effect hyperbolic, space. Its spatio-temporal behaviour is governed by nonlinear integro-differential equations defined on the Poincaré disc model of the two-dimensional hyperbolic space. Using methods from the theory of functional analysis we show the existence and uniqueness of a solution of these equations. In the case of stationary, i.e. time independent, solutions we perform a stability analysis which yields important results on their behavior. We also present an original study, based on non-Euclidean, hyperbolic, analysis, of a spatially localised bump solution in a limiting case. We illustrate our theoretical results with numerical simulations. AMS subject classifications: 30F45, 33C05, 34A12, 34D20, 34D23, 34G20, 37M05, 43A85, 44A35, 45G10, 51M10, 92B20, 92C20.


Introduction
The selectivity of the responses of individual neurons to external features is often the basis of neuronal representations of the external world. For example, neurons in the primary visual cortex (V1) respond preferentially to visual stimuli that have a specific orientation [1][2][3], spatial frequency [4], velocity and direction of motion [5], color [6]. A local network in the primary visual cortex, roughly 1 mm 2 of cortical surface, is assumed to consist of subgroups of inhibitory and excitatory neurons each of which is tuned to a particular feature of an external stimulus. These subgroups are the so-called Hubel and Wiesel hypercolumns of V1. We have introduced in [7] a new approach to model the processing of image edges and textures in the hypercolumns of area V1 that is based on a nonlinear representation of the image first order derivatives called the structure tensor [8,9]. We suggested that this structure tensor was represented by neuronal populations in the hypercolumns of V1. We also suggested that the time evolution of this representation was governed by equations similar to those proposed by Wilson and Cowan [10]. The question of whether some populations of neurons in V1 can represent the structure tensor is discussed in [7] but cannot be answered in a definite manner. Nevertheless, we hope that the predictions of the theory we are developing will help deciding on this issue.
Our present investigations were motivated by the work of Bressloff, Cowan, Golubitsky, Thomas and Wiener [11,12] on the spontaneous occurence of hallucinatory patterns under the influence of psychotropic drugs, and its extension to the structure tensor model. A further motivation was the following studies of Bressloff and Cowan [4,13,14] where they study a spatial extension of the ring model of orientation of Ben-Yishai [1] and Hansel, Sompolinsky [2]. To achieve this goal, we first have to better understand the local model, that is the model of a 'texture' hypercolumn isolated from its neighbours.
The aim of this paper is to present a rigorous mathematical framework for the modeling of the representation of the structure tensor by neuronal populations in V1. We would also like to point out that the mathematical analysis we are developing here, is general and could be applied to other integro-differential equations defined on the set of structure tensors, so that even if the structure tensor were found to be not represented in a hypercolumn of V1, our framework would still be relevant. We then concentrate on the occurence of localized states, also called bumps. This is in contrast to the work of [7] and [15] where 'spatially' periodic solutions were considered. The structure of this paper is as follows. In Section 2 we introduce the structure tensor model and the corresponding equations. We also link our model to the ring model of orientations. In Section 3 we use classical tools of evolution equations in functional spaces to analyse the problem of the existence and uniqueness of the solutions of our equations. In Section 4 we study stationary solutions which are very important for the dynamics of the equation by analysing a nonlinear convolution operator and making use of the Haar measure of our feature space. In Section 5, we push further the study of stationary solutions in a special case and we present a technical analysis involving hypergeometric functions of what we call a hyperbolic radially symmetric stationary-pulse in the high gain limit. Finally, in Section 6, we present some numerical simulations of the solutions to verify the findings of the theoretical results.

The model
By definition, the structure tensor is based on the spatial derivatives of an image in a small area that can be thought of as part of a receptive field. These spatial derivatives are then summed nonlinearly over the receptive field. Let I (x, y) denote the original image intensity function, where x and y are two spatial coordinates. Let I σ 1 denote the scale-space representation of I obtained by convolution with the Gaussian kernel g σ (x, y) = 1 2πσ 2 e −(x 2 +y 2 )/(2σ 2 ) : The gradient ∇I σ 1 is a two-dimensional vector of coordinates I σ 1 x , I σ 1 y which emphasizes image edges. One then forms the 2 × 2 symmetric matrix of rank one T 0 = ∇I σ 1 (∇I σ 1 ) T , where T indicates the transpose of a vector. The set of 2 × 2 symmetric positive semidefinite matrices of rank one will be noted S + (1, 2) throughout the paper (see [16] for a complete study of the set S + (p, n) of n × n symmetric positive semidefinite matrices of fixed-rank p < n). By convolving T 0 componentwise with a Gaussian g σ 2 we finally form the tensor structure as the symmetric matrix: where we have set for example: Since the computation of derivatives usually involves a stage of scale-space smoothing, the definition of the structure tensor requires two scale parameters. The first one, defined by σ 1 , is a local scale for smoothing prior to the computation of image derivatives. The structure tensor is insensitive to noise and details at scales smaller than σ 1 . The second one, defined by σ 2 , is an integration scale for accumulating the nonlinear operations on the derivatives into an integrated image descriptor. It is related to the characteristic size of the texture to be represented, and to the size of the receptive fields of the neurons that may represent the structure tensor.
By construction, T is symmetric and non negative as det(T ) ≥ 0 by the inequality of Cauchy-Schwarz, then it has two orthonormal eigenvectors e 1 , e 2 and two non negative corresponding eigenvalues λ 1 and λ 2 which we can always assume to be such that λ 1 ≥ λ 2 ≥ 0. Furthermore the spatial averaging distributes the information of the image over a neighborhood, and therefore the two eigenvalues are always positive. Thus, the set of the structure tensors lives in the set of 2 × 2 symmetric positive definite matrices, noted SPD(2, R) throughout the paper. The distribution of these eigenvalues in the (λ 1 , λ 2 ) plane reflects the local organization of the image intensity variations. Indeed, each structure tensor can be written as the linear combination: T = λ 1 e 1 e T 1 + λ 2 e 2 e T 2 = (λ 1 − λ 2 )e 1 e T 1 + λ 2 e 1 e T 1 + e 2 e T 2 = (λ 1 − λ 2 )e 1 e T 1 + λ 2 I 2 , where I 2 is the identity matrix and e 1 e T 1 ∈ S + (1, 2). Some easy interpretations can be made for simple examples: constant areas are characterized by λ 1 = λ 2 ≈ 0, straight edges are such that λ 1 λ 2 ≈ 0, their orientation being that of e 2 , corners yield λ 1 ≥ λ 2 0. The coherency c of the local image is measured by the ratio c = λ 1 −λ 2 λ 1 +λ 2 , large coherency reveals anisotropy in the texture.
We assume that a hypercolumn of V1 can represent the structure tensor in the receptive field of its neurons as the average membrane potential values of some of its membrane populations. Let T be a structure tensor. The time evolution of the average potential V (T , t) for a given column is governed by the following neural mass equation adapted from [7] where we allow the connectivity function W to depend upon the time variable t and we integrate over the set of 2 × 2 symmetric definite-positive matrices: The nonlinearity S is a sigmoidal function which may be expressed as: where μ describes the stiffness of the sigmoid. I ext is an external input.
The set SPD(2) can be seen as a foliated manifold by way of the set of special symmetric positive definite matrices SSPD(2) = SPD(2) ∩ SL(2, R). Indeed, where D is the Poincaré Disk, see, for example, [7]. As a consequence we use the following foliation of SPD(2) : SPD(2) hom = D × R + * , which allows us to write for all T ∈ SPD(2), T = (z, ) with (z, ) ∈ D × R + * . T , z and are related by the relation det(T ) = 2 and the fact that z is the representation in D ofT ∈ SSPD(2) with T = T .
It is well-known [17] that D (and hence SSPD (2)) is a two-dimensional Riemannian space of constant sectional curvature equal to −1 for the distance noted d 2 defined by The isometries of D, that are the transformations that preserve the distance d 2 are the elements of unitary group U(1, 1). In Appendix A we describe the basic structure of this group. It follows, for example, [7,18], that SDP(2) is a three-dimensional Riemannian space of constant sectional curvature equal to −1 for the distance noted d 0 defined by As shown in Proposition B.0.1 of Appendix B it is possible to express the volume element dT in (z 1 , z 2 , ) coordinates with z = z 1 + iz 2 : We note dm(z) = dz 1 dz 2 (1−|z| 2 ) 2 and equation (2) can be written in (z, ) coordinates: We get rid of the constant 8 √ 2 by redefining W as 8 In [7], we have assumed that the representation of the local image orientations and textures is richer than, and contains, the local image orientations model which is conceptually equivalent to the direction of the local image intensity gradient. The richness of the structure tensor model has been expounded in [7]. The embedding of the ring model of orientation in the structure tensor model can be explained by the intrinsic relation that exists between the two sets of matrices SPD(2, R) and S + (1, 2). First of all, when σ 2 goes to zero, that is when the characteristic size of the structure becomes very small, we have T 0 g σ 2 → T 0 , which means that the tensor T ∈ SPD(2, R) degenerates to a tensor T 0 ∈ S + (1, 2), which can be interpreted as the loss of one dimension. We can write each T 0 ∈ S + (1, 2) as T 0 = xx T = r 2 uu T , where u = (cos θ, sin θ) T and (r, θ ) is the polar representation of x. Since, x and −x correspond to the same T 0 , θ is equated to θ + kπ , k ∈ Z. Thus S + (1, 2) = R + * × P 1 , where P 1 is the real projective space of dimension 1 (lines of R 2 ). Then the integration scale σ 2 , at which the averages of the estimates of the image derivatives are computed, is the link between the classical representation of the local image orientations by the gradient and the representation of the local image textures by the structure tensor. It is also possible to highlight this explanation by coming back to the interpretation of straight edges of the previous paragraph. When λ 1 λ 2 ≈ 0 then T ≈ (λ 1 − λ 2 )e 1 e T 1 ∈ S + (1, 2) and the orientation is that of e 2 . We denote by P the projection of a 2 × 2 symmetric definite positive matrix on the set S + (1, 2) defined by: where T is as in equation (1). We can introduce a metric on the set S + (1, 2) which is derived from a well-chosen Riemannian quotient geometry (see [16]). The resulting Riemannian space has strong geometrical properties: it is geodesically complete and the metric is invariant with respect to all transformations that preserve angles (orthogonal transformations, scalings and pseudoinversions). Related to the decomposition S + (1, 2) = R + * × P 1 , a metric on the space S + (1, 2) is given by: The space S + (1, 2) endowed with this metric is a Riemannian manifold (see [16]). Finally, the distance associated to this metric is given by: where we normalize to 1 the volume element for the θ coordinate. Let now τ = P(T ) be a symmetric positive semidefinite matrix. The average potential V (τ, t) of the column has its time evolution that is governed by the following neural mass equation which is just a projection of equation (2) on the subspace S + (1, 2): In (r, θ ) coordinates, (4) is rewritten as: This equation is richer than the ring model of orientation as it contains an additional information on the contrast of the image in the orthogonal direction of the prefered orientation. If one wants to recover the ring model of orientation tuning in the visual cortex as it has been presented and studied by [1,2,19], it is sufficient i) to assume that the connectivity function is time-independent and has a convolutional form: and ii) to look at semi-homogeneous solutions of equation (4), that is, solutions which do not depend upon the variable r. We finally obtain: where: It follows from the above discussion that the structure tensor contains, at a given scale, more information than the local image intensity gradient at the same scale and that it is possible to recover the ring model of orientations from the structure tensor model.
The aim of the following sections is to establish that (3) is well-defined and to give necessary and sufficient conditions on the different parameters in order to prove some results on the existence and uniqueness of a solution of (3).

The existence and uniqueness of a solution
In this section we provide theoretical and general results of existence and uniqueness of a solution of (2). In the first subsection (Section 3.1) we study the simpler case of the homogeneous solutions of (2), that is, of the solutions that are independent of the tensor variable T . This simplified model allows us to introduce some notations for the general case and to establish the useful Lemma 3.1.1. We then prove in Section 3.2 the main result of this section, that is the existence and uniqueness of a solution of (2). Finally we develop the useful case of the semi-homogeneous solutions of (2), that is, of solutions that depend on the tensor variable but only through its z coordinate in D.

Homogeneous solutions
A homogeneous solution to (2) is a solution V that does not depend upon the tensor variable T for a given homogenous input I (t) and a constant initial condition V 0 . In (z, ) coordinates, a homogeneous solution of (3) is defined by: where: Hence necessary conditions for the existence of a homogeneous solution are that: • the double integral (6) is convergent, • W (z, , t) does not depend upon the variable (z, ). In that case, we write W (t) instead of W (z, , t).
In the special case where W (z, , z , , t) is a function of only the distance d 0 between (z, ) and (z , ): the second condition is automatically satisfied. The proof of this fact is given in Lemma D.0.2 of Appendix D. To summarize, the homogeneous solutions satisfy the differential equation:

A first existence and uniqueness result
Equation (3) defines a Cauchy's problem and we have the following theorem.
Proof It is a direct application of Cauchy's theorem on differential equations. We consider the mapping f : J × R → R defined by: It is clear that f is continuous from J × R to R. We have for all x, y ∈ R and t ∈ J : Since, W is continuous on the compact interval J , it is bounded there by C > 0 and: We can extend this result to the whole time real line if I and W are continuous on R. Proof We have already shown the following inequality: Then f is locally Lipschitz with respect to its second argument. Let V be a maximal solution on J 0 and we denote by β the upper bound of J 0 . We suppose that β < +∞.
Then we have for all t ≥ 0: This implies that the maximal solution V is bounded for all t ∈ [0, β], but Theorem C.0.2 of Appendix C ensures that it is impossible. Then, it follows that necessarily β = +∞.

Simplification of (6) in a special case
Invariance In the previous section, we have stated that in the special case where W was a function of the distance between two points in D × R + * , then W (z, , t) did not depend upon the variables (z, ). As already said in the previous section, the following result holds (see proof of Lemma D.0.2 of Appendix D). Mexican hat connectivity In this paragraph, we push further the computation of W in the special case where W does not depend upon the time variable t and takes the special form suggested by Amari in [20], commonly referred to as the 'Mexican hat' connectivity. It features center excitation and surround inhibition which is an effective model for a mixed population of interacting inhibitory and excitatory neurons with typical cortical connections. It is also only a function of d 0 (T , T ).
In detail, we have: where: In this case we can obtain a very simple closed-form formula for W as shown in the following lemma.
where erf is the error function defined as: Proof The proof is given in Lemma E.0.3 of Appendix E.

General solution
We now present the main result of this section about the existence and uniqueness of solutions of equation (2). We first introduce some hypotheses on the connectivity function W . We present them in two ways: first on the set of structure tensors considered as the set SPD(2), then on the set of tensors seen as D × R + * . Let J be a subinterval of R. We assume that: Equivalently, we can express these hypotheses in (z, ) coordinates:

Functional space setting
We introduce the following mapping f g : (t, φ) → f g (t, φ) such that: Our aim is to find a functional space F where (3) is well-defined and the function f g maps F to F for all ts. A natural choice would be to choose φ as a L p (D × R + * )-integrable function of the space variable with 1 ≤ p < +∞. Unfortunately, the homogeneous solutions (constant with respect to (z, )) do not belong to that space. Moreover, a valid model of neural networks should only produce bounded membrane potentials. That is why we focus our choice on the functional space F = L ∞ (D × R + * ). As D × R + * is an open set of R 3 , F is a Banach space for the norm: φ F = sup z∈D sup ∈R + * |φ(z, )|.

The existence and uniqueness of a solution of (3)
We rewrite (3) as a Cauchy problem: Proof We prove that f g is continuous on J × F . We have and therefore Because of condition (H2) we can choose |t − s| small enough so that W(t) − W(s) L 1 is arbitrarily small. This proves the continuity of f g . Moreover it follows from the previous inequality that: This ensures the Lipschitz continuity of f g with respect to its second argument, uniformly with respect to the first. The Cauchy-Lipschitz theorem on a Banach space yields the conclusion. [21]. The main differences are that, first, we allow the connectivity function to depend upon the time variable t and, second, that our space features is no longer a R n but a Riemanian manifold. In their article, Potthast and Graben also work with a different functional space by assuming more regularity for the connectivity function W and then obtain more regularity for their solutions.

Remark 3.2.1 Our result is quite similar to those obtained by Potthast and Graben in
Proof We have just seen in the previous proof that f g is globally Lipschitz with respect to its second argument: then Theorem C.0.3 of Appendix C gives the conclusion.

The intrinsic boundedness of a solution of (3)
In the same way as in the homogeneous case, we show a result on the boundedness of a solution of (3).

Proposition 3.2.3
If the external current I ext belongs to C(R + , F ) and is bounded in time sup t∈R + I ext (t) F < +∞ and W satisfies hypotheses (H1bis)-(H3bis) with J = R + , then the solution of (10) is bounded for each initial condition V 0 ∈ F . Let us set: Proof Let V be a solution defined on R + . Then we have for all t ∈ R + * : The following upperbound holds We can rewrite (11) as: The inequality (12) shows that for t large enough this yields a contradiction. Therefore there exists t 0 > 0 such that and hence The following corollary is a consequence of the previous proposition.

Semi-homogeneous solutions
A semi-homogeneous solution of (3) is defined as a solution which does not depend upon the variable . In other words, the populations of neurons is not sensitive to the determinant of the structure tensor, that is to the contrast of the image intensity. The neural mass equation is then equivalent to the neural mass equation for tensors of unit determinant. We point out that semi-homogeneous solutions were previously introduced in [7] where a bifurcation analysis of what they called H-planforms was performed. In this section, we define the framework in which their equations make sense without giving any proofs of our results as it is a direct consequence of those proven in the general case. We rewrite equation (3) in the case of semi-homogeneous solutions: where We have implicitly made the assumption, that W sh does not depend on the coordinate . Some conditions under which this assumption is satisfied are described below and are the direct transductions of those of the general case in the context of semihomogeneous solutions.
Let J be an open interval of R. We assume that: Note that conditions (C1)-(C2) and Lemma 3.1.1 imply that for all z ∈ D, From now on, F = L ∞ (D) and the Fischer-Riesz's theorem ensures that L ∞ (D) is a Banach space for the norm: ψ ∞ = inf{C ≥ 0, |ψ(z)| ≤ C for almost every z ∈ D}. This solution, defined on the subinterval J of R can in fact be extended to the whole real line, and we have the following proposition. We can also state a result on the boundedness of a solution of (13): The open ball B ρ of F of center 0 and radius ρ is stable under the dynamics of equation (13). Moreover it is an attracting set for this dynamics

Stationary solutions
We look at the equilibrium states, noted V 0 μ of (3), when the external input I and the connectivity W do not depend upon the time. We assume that W satisfies hypotheses (H1bis)-(H2bis). We redefine for convenience the sigmoidal function to be: so that a stationary solution (independent of time) satisfies: We define the nonlinear operator from F to F , noted G μ , by: Finally, (14) is equivalent to:

Study of the nonlinear operator G μ
We recall that we have set for the Banach space F = L ∞ (D × R + * ) and Proposition 3.2.1 shows that G μ : F → F . We have the further properties: Proposition 4.1.1 G μ satisfies the following properties: The first property was shown to be true in the proof of Theorem 3.3.1. The second property follows from the following inequality: We denote by G l and G ∞ the two operators from F to F defined as follows for all V ∈ F and all (z, ) ∈ D × R + * : and where H is the Heaviside function.
It is straightforward to show that both operators are well-defined on F and map F to F . Moreover the following proposition holds.
Proof It is a direct application of the dominated convergence theorem using the fact that:

The convolution form of the operator G μ in the semi-homogeneous case
It is convenient to consider the functional space F sh = L ∞ (D) to discuss semihomogeneous solutions. A semi-homogeneous persistent state of (3) is deduced from (14) and satisfies: where the nonlinear operator G sh μ from F sh to F sh is defined for all V ∈ F sh and z ∈ D by: We define the associated operators, G sh l , G sh ∞ : We rewrite the operator G sh μ in a convenient form by using the convolution in the hyperbolic disk. First, we define the convolution in a such space. Let O denote the center of the Poincaré disk that is the point represented by z = 0 and dg denote the Haar measure on the group G = SU(1, 1) (see [22] and Appendix A), normalized by: for all functions of L 1 (D). Given two functions f 1 , f 2 in L 1 (D) we define the convolution * by: Proof We only prove the result for G μ . Let z ∈ D, then: and for all g ∈ SU(1, 1), d 2 (z, z ) = d 2 (g · z, g · z ) so that: Let b be a point on the circle ∂D. For z ∈ D, we define the 'inner product' z, b to be the algebraic distance to the origin of the (unique) horocycle based at b through z (see [7]). Note that z, b does not depend on the position of z on the horocycle. The Fourier transform in D is defined as (see [22]): Proof For all λ ∈ R and b = e iθ ∈ ∂D, We recall that for all φ ∈ R r φ is the rotation of angle φ and we have W sh (r φ · z) = W sh (z), dm(z) = dm(r φ · z) and z, b = r φ · z, r φ · b , then: We now introduce two functions that enjoy some nice properties with respect to the Hyperbolic Fourier transform and are eigenfunctions of the linear operator G sh l .
Proof We begin with b = 1 ∈ ∂D and use the horocyclic coordinates. We use the same changes of variables as in Lemma 3.1.1: By rotation, we obtain the property for all b ∈ ∂D.
For the second property [22,Lemma 4.7] shows that: A consequence of this proposition is the following lemma.

Lemma 4.2.2
The linear operator G sh l is not compact and for all μ ≥ 0, the nonlinear operator G sh μ is not compact.
Proof The previous Proposition 4.2.2 shows that G sh l has a continuous spectrum which iimplies that is not a compact operator.
Let U be in F sh , for all V ∈ F sh we differentiate G sh μ and compute its Frechet derivative: If we assume further that U does not depend upon the space variable z, U(z) = U 0 we obtain: If G sh μ was a compact operator then its Frechet derivative D(G sh μ ) U 0 would also be a compact operator, but it is impossible. As a consequence, G sh μ is not a compact operator.

The convolution form of the operator G μ in the general case
We adapt the ideas presented in the previous section in order to deal with the general case. We recall that if H is the group of positive real numbers with multiplication as operation, then the Haar measure dh is given by dx x . For two functions f 1 , f 2 in L 1 (D × R + * ) we define the convolution by: We recall that we have set by definition: W(z, ) = W (d 2 (z, 0), | log( )|).
Proof Let (z, ) be in D × R + * . We follow the same ideas as in Proposition 4.2.1 and prove only the first result. We have We next assume further that the function W is separable in z, and more precisely that W(z, ) = W 1 (z)W 2 (log( )) where W 1 (z) = W 1 (d 2 (z, 0)) and W 2 (log( )) = W 2 (| log( )|) for all (z, ) ∈ D × R + * . The following proposition is an echo to Proposition 4.2.2.
where W 2 is the usual Fourier transform of W 2 .
Proof The proof of this proposition is exactly the same as for Proposition 4.2.2. Indeed: A straightforward consequence of this proposition is an extension of Lemma 4.2.2 to the general case:

The set of the solutions of (14)
Let B μ be the set of the solutions of (14) for a given slope parameter μ: We have the following proposition. Proof Due to the properties of the sigmoid function, there always exists a constant solution in the case where I ext is constant. In the general case where I ext ∈ F , the statement is a direct application of the Banach fixed point theorem, as in [23].  U(1, 1), the group of the isometries of D (see Appendix A), and the condition μS m W g 0 < α is satisfied, then the unique stationary solution will also be invariant under the action of the same subgroup. We refer the interested reader to our work [15] on equivariant bifurcation of hyperbolic planforms on the subject.
When the condition μS m W g 0 < α is satisfied we call primary stationary solution the unique solution in B μ .

Stability of the primary stationary solution
In this subsection we show that the condition μS m W g 0 < α guarantees the stability of the primary stationary solution to (3). Theorem 4.5. 1 We suppose that I ∈ F and that the condition μS m W g 0 < α is satisfied, then the associated primary stationary solution of (3) is asymtotically stable.
∂ t X(z, , t) = −αX(z, , t) ⇒ ∂ t e αt X(z, , t) If we set: G(t) = e αt X(t) ∞ , then we have: and G is continuous for all t ≥ 0. The Gronwall inequality implies that: and the conclusion follows.

Spatially localised bumps in the high gain limit
In many models of working memory, transient stimuli are encoded by featureselective persistent neural activity. Such stimuli are imagined to induce the formation of a spatially localised bump of persistent activity which coexists with a stable uniform state. As an example, Camperi and Wang [24] have proposed and studied a network model of visuo-spatial working memory in prefontal cortex adapted from the ring model of orientation of Ben-Yishai and colleagues [1]. Many studies have emerged in the past decades to analyse these localised bumps of activity [25][26][27][28][29], see the paper by Coombes for a review of the domain [30]. In [25,26,28], the authors have examined the existence and stability of bumps and multi-bumps solutions to an integro-differential equation describing neuronal activity along a single spatial domain. In [27,29] the study is focused on the two-dimensional model and a method is developed to approximate the integro-differential equation by a partial differential equation which makes possible the determination of the stability of circularly symmetric solutions. It is therefore natural to study the emergence of spatially localized bumps for the structure tensor model in a hypercolumn of V1. We only deal with the reduced case of equation (13) which means that the membrane activity does not depend upon the contrast of the image intensity, keeping the general case for future work.
In order to construct exact bump solutions and to compare our results to previous studies [25][26][27][28][29], we consider the high gain limit μ → ∞ of the sigmoid function. As above we denote by H the Heaviside function defined by H (x) = 1 for x ≥ 0 and H (x) = 0 otherwise. Equation (13) is rewritten as: We have introduced a threshold κ to shift the zero of the Heaviside function. We make the assumption that the system is spatially homogeneous that is, the external input I does not depend upon the variables t and the connectivity function depends only on the hyperbolic distance between two points of D : W (z, z ) = W (d 2 (z, z )). For illustrative purposes, we will use the exponential weight distribution as a specific example throughout this section: The theoretical study of equation (20) has been done in [21] where the authors have imposed strong regularity assumptions on the kernel function W, such as Hölder continuity, and used compactness arguments and integral equation techniques to obtain a global existence result of solutions to (20). Our approach is very different, we follow that of [25][26][27][28][29]31] by proceeding in a constructive fashion. In a first part, we define what we call a hyperbolic radially symmetric bump and present some preliminary results for the linear stability analysis of the last part. The second part is devoted to the proof of a technical Theorem 5.1.1 which is stated in the first part. The proof uses results on the Fourier transform introduced in Section 4, hyperbolic geometry and hypergeometric functions. Our results will be illustrated in the following Section 6. For convenience, we note M(z, K) the integral K W (z, z ) dm(z ) with K = {z ∈ D|V (z) ≥ κ}. The relation V (z) = κ holds for all z ∈ ∂K.

Definition 5.1.1 V is called a hyperbolic radially symmetric stationary-pulse solution of (20) if V depends only upon the variable r and is such that:
and is a fixed point of equation (20): where I ext (r) = Ie From symmetry arguments there exists a hyperbolic radially symmetric stationarypulse solution V (r) of (20), furthermore the threshold κ and width ω are related according to the self-consistency condition where The existence of such a bump can then be established by finding solutions to (23) The function N(ω) is plotted in Figure 1  Anticipating the stability results of Section 5.3, we obtain that when N (ω) < 0 then the corresponding solution is stable.
We end this subsection with the usefull and technical following formula.

where W(λ) is the Fourier Helgason transform of W(z)
with α + β + 1 = ρ and F is the hypergeometric function of first kind.

Remark 5.1.3 Let us point out that this result can be linked to the work of Folias and
Bressloff in [31] and then used in [29]. They constructed a two-dimensional pulse for a general, radially symmetric synaptic weight function.

where J ν (x) is the Bessel function of the first kind andw is the real Fourier transform of w. In our case, instead of the Bessel function, we find (ν,ν) λ (r) which is linked to the hypergeometric function of the first kind.
We now show that for a general monotonically decreasing weight function W , the function M(r, ω) is necessarily a monotonically decreasing function of r. This will ensure that the hyperbolic radially symmetric stationary-pulse solution (22) is also a monotonically decreasing function of r in the case of a Gaussian input. The demonstration of this result will directly use Theorem 5.1.1. It is result of elementary hyperbolic trigonometry that d 2 tanh(r), tanh(r )e iθ = tanh −1 tanh(r) 2 + tanh(r ) 2 − 2 tanh(r) tanh(r ) cos(θ ) 1 + tanh(r) 2 tanh(r ) 2 − 2 tanh(r) tanh(r ) cos(θ ) (25) we let ρ = tanh(r), ρ = tanh(r ) and define It follows that , We conclude that if ρ > tanh(ω) then for all 0 ≤ ρ ≤ tanh(ω) and 0 ≤ θ ≤ 2π 2 ρ − ρ cos(θ ) + 2ρρ ρ − ρ cos(θ ) > 0, which implies M(r, ω) < 0 for r > ω, since W < 0. To see that it is also negative for r < ω, we differentiate equation (24) with respect to r: The following formula holds for the hypergeometric function (see Erdelyi in [32]): Substituting in the previous equation giving ∂M ∂r we find: implying that: Consequently, ∂M ∂r (r, ω) < 0 for r < ω. Hence V is monotonically decreasing in r for any monotonically decreasing synaptic weight function W .
As a consequence, for our particular choice of exponential weight function (21), the radially symmetric bump is monotonically decreasing in r, as it will be recover in our numerical experiments in Section 6.

Proof of Theorem 5.1.1
The proof of Theorem 5.1.1 goes in four steps. First we introduce some notations and recall some basic properties of the Fourier transform in the Poincaré disk. Second we prove two propositions. Third we state a technical lemma on hypergeometric functions, the proof being given in Lemma F.0.4 of Appendix F. The last step is devoted to the conclusion of the proof.

First step
In order to calculate M(r, ω), we use the Fourier transform in D which has already been introduced in Section 4. First we rewrite M(r, ω) as a convolution product: Proof We start with the definition of M(r, ω) and use the convolutional form of the integral: In [22], Helgason proves an inversion formula for the hyperbolic Fourier transform and we apply this result to W: the last equality is a direct application of Lemma 4.2.1 and we can deduce that Finally we have: which is the desired formula.
It appears that the study of M(r, ω) consists in calculating the convolution product

Proposition 5.2.2 For all
Proof Let z = k · O for k ∈ G we have:

Second step
In this part, we prove two results: is a radial function, that is, it depends only upon the variable r. • the following equality holds for z = tanh(r)e iθ :

Proposition 5.2.3 If z = k · O and z is written tanh(r)e iθ with r = d 2 (z, O) in hyperbolic polar coordinates the function λ * 1 B h (0,ω) (z) depends only upon the variable r.
Proof If z = tanh(r)e iθ , then z = rot θ a r · O and k −1 = a −r rot −θ . Similarly z = rot θ a r · O. We can write thanks to the previous Proposition 5.2.2: which, as announced, is only a function of r.
We now give an explicit formula for the integral B h (0,ω) λ (a −r · z ) dm(z ).

Proposition 5.2.4
For all z = tanh(r)e iθ we have: Proof We first recall a formula from [22].

Lemma 5.2.1 For all g ∈ G the following equation holds:
Proof See [22].
It follows immediately that for all z ∈ D and r ∈ R we have: We integrate this formula over the hyperbolic ball B h (0, ω) which gives: and we exchange the order of integration: We note that the integral B h (0,ω) e (iλ+1) z ,b dm(z ) does not depend upon the variable b = e iφ . Indeed: and indeed the integral does not depend upon the variable b: Finally, we can write: because λ = −λ (as solutions of the same equation). This completes the proof that:

Third step
We state a useful formula.

Lemma 5.2.2
For all ω > 0 the following formula holds: Proof See Lemma F.0.4 of Appendix F.

The main result
At this point we have proved the following proposition thanks to Propositions 5.2.1 and 5.2.4.

Proposition 5.2.5
If z = tanh(r)e iθ ∈ B h (0, ω), M(r, ω) is given by the following formula: We are now in a position to obtain the analytic form for M(r, ω) of Theorem 5.1.1. We prove that Indeed, in hyperbolic polar coordinates, we have: On the other hand: This yields and we use Lemma 5.2.2 to establish (24).

Linear stability analysis
We now analyse the evolution of small time-dependent perturbations of the hyperbolic stationary-pulse solution through linear stability analysis. We use classical tools already developped in [29,31].

Spectral analysis of the linearized operator
Equation (20) is linearized about the stationary solution V (r) by introducing the timedependent perturbation: This leads to the linear equation: We separate variables by setting φ(z, t) = φ(z)e βt to obtain the equation: Introducing the hyperbolic polar coordinates z = tanh(r)e iθ and using the result: we obtain: Note that we have formally differentiated the Heaviside function, which is permissible since it arises inside a convolution. One could also develop the linear stability analysis by considering perturbations of the threshold crossing points along the lines of Amari [20]. Since we are linearizing about a stationary rather than a traveling pulse, we can analyze the spectrum of the linear operator without the recourse to Evans functions.
Discrete spectrum If we are not in the previous case we have to study the solutions of the integral equation (28). This equation shows that φ(r, θ) is completely determined by its values φ(ω, θ) on the circle of equation r = ω. Hence, we need only to consider r = ω, yielding the integral equation: The solutions of this equation are exponential functions e γ θ , where γ satisfies: By the requirement that φ is 2π -periodic in θ , it follows that γ = in, where n ∈ Z.
Thus the integral operator with kernel W has a discrete spectrum given by: β n is real since: Hence, We can state the folliwing proposition: Provided that for all n ≥ 0, β n < 0 then the hyperbolic stationary pulse is stable.
We now derive a reduced condition linking the parameters for the stability of hyperbolic stationary pulse.
Reduced condition Since W • tanh −1 (r) is a positive function of r, it follows that: Stability of the hyperbolic stationary pulse requires that for all n ≥ 0, β n < 0. This can be rewritten as: Using the fact that β n ≤ β 0 for all n ≥ 1, we obtain the reduced stability condition: From (22) we have: We have previously established that M r (ω) > 0 and I (ω) is negative by definition. Hence, letting D(ω) = |I (ω)|, we have By substitution we obtain another form of the reduced stability condition: We also have: and showing that the stability condition (29) is satisfied when N (ω) < 0 and is not satisfied when N (ω) > 0.

Numerical results
The aim of this section is to numerically solve (13) for different values of the parameters. This implies developing a numerical scheme that approaches the solution of our equation, and proving that this scheme effectively converges to the solution.
Since equation (13) is defined on D, computing the solutions on the whole hyperbolic disk has the same complexity as computing the solutions of usual Euclidean neural field equations defined on R 2 . As most authors in the Euclidean case [26,27,29,31], we reduce the domain of integration to a compact region of the hyperbolic disk. Practically, we work in the Euclidean ball of radius a = 0.5 and center 0. Note that a Euclidean ball centered at the origin is also a centered hyperbolic ball, their radii being different.
We have divided this section into four parts. The first part is dedicated to the study of the discretization scheme of equation (13). In the following two parts, we study the solutions for different connectivity functions: an exponential function, Section 6.2, and a difference of Gaussians, Section 6.3.

Discretization scheme
We discretize R in order to turn (30) into a finite number of equations. For this and obtain the (N + 1)(M + 1) equations: which define the discretization of (30): whereṼ (t) ∈ M N +1,M+1 (R),Ṽ (t) i,j = V (r i , θ j , t). Similar definitions apply toĨ andṼ 0 . Moreover: M n,p (R) is the space of the matrices of size n × p with real coefficients. It remains to discretize the integral term. For this as in [33], we use the rectangular rule for the quadrature so that for all (r, θ ) ∈ R we have: We end up with the following numerical scheme, where V i,j (t) (resp. I i,j (t)) is an approximation ofṼ i,j (t) (resp.Ĩ i,j ),

Discussion
We discuss the error induced by the rectangular rule for the quadrature. Let f be a function which is f C 2 where m and n are the number of subintervals used and f C 2 = |α|≤2 sup [a,b]×[c,d] |∂ α f | where, as usual, α is a multi-index. As a consequence, if we want to control the error, we have to impose that the solution is, at least, C 2 in space.
Four our numerical experiments we use the specific function ode45 of Matlab which is based on an explicit Runge-Kutta formula (see [34] for more details on Runge-Kutta methods).
We can also establish a proof of the convergence of the numerical scheme which is exactly the same as in [33] excepted that we use the theorem of continuous dependence of the solution for ordinary differential equations.

Purely excitatory exponential connectivity function
In this subsection, we give some numerical solutions of (13) in the case where the connectivity function is an exponential function, w(x) = e − |x| b , with b a positive parameter. Only excitation is present in this case. In all the experiments we set α = 0.1 and S(x) = 1 1+e −μx with μ = 10.

Constant input
We fix the external input I (z) to be of the form: In all experiments we set I = 0.1 and σ = 0.05, this means that the input has a sharp profile centered at 0. We show in Figure 2 plots of the solution at time T = 2,500 for three different values of the width b of the exponential function. When b = 1, the whole network is highly excited, whereas as b changes from 1 to 0.1 the amplitude of the solution decreases, and the area of high excitation becomes concentrated around the external input.

Variable input
In this paragraph, we allow the external current to depend upon the time variable. We have: where z 0 (t) = r 0 e i 0 t . This is a bump rotating with angular velocity 0 around the circle of radius r 0 centered at the origin. In our numerical experiments we set r 0 = 0.4, 0 = 0.01, I = 0.1 and σ = 0.05. We plot in Figure 3 the solution at different times T = 100, 150, 200, 250.

High gain limit
We consider the high gain limit μ → ∞ of the sigmoid function and we propose to illustrate Section 5 with a numerical simulation. We set α = 1, κ = 0.04, ω = 0.18. We fix the input to be of the form: with I = 0.04 and σ = 0.05. Then the condition of existence of a stationary pulse (23) is satisfied, see Figure 1. We plot a bump solution according to (23) in Figure 4.  (13) in the case of an exponential connectivity function with b = 0.1 at different times with a time-dependent input, see text.

Excitatory and inhibitory connectivity function
We give some numerical solutions of (13) in the case where the connectivity function is a difference of Gaussians, which features an excitatory center and an inhibitory surround: We illustrate the behaviour of the solutions when increasing the slope μ of the sigmoid. We set the sigmoid S(x) = 1 1+e −μx − 1 2 so that it is equal to 0 at the origin and we choose the external input equal to zero, I (z, t) = 0. In this case the constant function equal to 0 is a solution of (13). For small values of the slope μ, the dynamics of the solution is trivial: every solution asymptotically converges to the null solution, as shown in top left hand corner of Figure 5 with μ = 1. When increasing μ, the stability bound, found in Section 4.5 is no longer satisfied and the null solution may no longer be stable. In effect this solution may bifurcate to other, more interesting solutions. We plot in Figure 5, some solutions at T = 2,500 for different values of μ (μ = 3, 5, 10, 20 and 30). We can see exotic patterns which feature some interesting symmetries. The formal study of these bifurcated solutions is left for future work.

Conclusion
In this paper, we have studied the existence and uniqueness of a solution of the evolution equation for a smooth neural mass model called the structure tensor model. This model is an approach to the representation and processing of textures and edges in the visual area V1 which contains as a special case the well-known ring model of orientations (see [1,2,19]). We have also given a rigorous functional framework for the study and computation of the stationary solutions to this nonlinear integrodifferential equation. This work sets the basis for further studies beyond the spatially periodic case studied in [15], where the hypothesis of spatial periodicity allows one to replace the unbounded (hyperbolic) domain by a compact one, hence making the functional analysis much simpler.
We have completed our study by constructing and analyzing spatially localised bumps in the high-gain limit of the sigmoid function. It is true that networks with Heaviside nonlinearities are not very realistic from the neurobiological perspective and lead to difficult mathematical considerations. However, taking the high-gain limit is instructive since it allows the explicit construction of stationary solutions which is impossible with sigmoidal nonlinearities. We have constructed what we called a hyperbolic radially symmetric stationary-pulse and presented a linear stability analysis adapted from [31]. The study of stationary solutions is very important as it conveys information for models of V1 that is likely to be biologically relevant. Moreover our study has to be thought of as the analog in the case of the structure tensor model to the analysis of tuning curves of the ring model of orientations (see [1,2,19,35]). However, these solutions may be destabilized by adding lateral spatial connections in a spatially organized network of structure tensor models; this remains an area of future investigation. As far as we know, only Bressloff and coworkers looked at this problem (see [3,4,[11][12][13][14]).
Finally, we illustrated our theoretical results with numerical simulations based on rigorously defined numerical schemes. We hope that our numerical experiments will lead to new and exciting investigations such as a thorough study of the bifurcations of the solutions of our equations with respect to such parameters as the slope of the sigmoid and the width of the connectivity function. for example, [17]. The direct isometries (preserving the orientation) in D are the elements of the special unitary group, noted SU(1, 1), of 2 × 2 Hermitian matrices with determinant equal to 1. Given: an element of SU(1, 1), the corresponding isometry γ in D is defined by: Orientation reversing isometries of D are obtained by composing any transformation (32) with the reflection κ : z →z. The full symmetry group of the Poincaré disc is therefore: Let us now describe the different kinds of direct isometries acting in D. We first define the following one parameter subgroups of SU(1, 1): Note that rot φ · z = e iφ z and also a r · O = tanh r, with O being the center of the Poincaré disk that is the point represented by z = 0. The group K is the orthogonal group O(2). Its orbits are concentric circles. It is possible to express each point z ∈ D in hyperbolic polar coordinates: z = rot φ a r · O = tanh re iφ and r = d 2 (z, 0).
The orbits of A converge to the same limit points of the unit circle ∂D, b ±1 = ±1 when r → ±∞. They are circular arcs in D going through the points b 1 and b −1 .
The orbits of N are the circles inside D and tangent to the unit circle at b 1 . These circles are called horocycles with base point b 1 . N is called the horocyclic group. It is also possible to express each point z ∈ D in horocyclic coordinates: z = n s a r · O, where n s are the transformations associated with the group N (s ∈ R) and a r the transformations associated with the subroup A (r ∈ R).

A.1 Iwasawa decomposition
The following decomposition holds, see [36]: This theorem allows us to decompose any isometry of D as the product of at most three elements in the groups, K, A and N . wherex i , i = 1, 2, 3, is given by: The determinant of the Jacobian of the transformation (x 1 , x 2 , x 3 ) → ( , z 1 , z 2 ) is found to be equal to: Hence, the volume element in ( , z 1 , z 2 ) coordinates is We suppose that f ∈ C(J × O, F ) and is locally Lipschitz with respect to its second argument. Then for all (t 0 , V 0 ) ∈ J ×O, there exists τ > 0 and V ∈ C 1 (]t τ , t 0 +τ [, O) unique solution of (34).
Lemma C.0.1 Under hypotheses of Theorem C.0.1, if V 1 ∈ C 1 (J 1 , O) and V 2 ∈ C 1 (J 2 , O) are two solutions and if there exists t 0 ∈ J 1 ∩ J 2 such that V 1 (t 0 ) = V 2 (t 0 ) then: This lemma shows the existence of a larger interval J 0 on which the initial value problem (34) has a unique solution. This solution is called the maximal solution. Theorem C.0. 3 We suppose f ∈ C(J × F , F ) and is globally Lipschitz with respect to its second argument. Then for all (t 0 , V 0 ) ∈ J × F , there exists a unique V ∈ C 1 (J, F ) solution of (34). Proof We work in (z, ) coordinates and we begin by rewriting the double integral (6) for all (z, ) ∈ R + * × D: The change of variable → yields: And it establishes that W does not depend upon the variable . To finish the proof, we show that the following integral does not depend upon the variable z ∈ D: where f is a real-valued function such that (z) is well defined. We express z in horocyclic coordinates: z = n s a r · O (see Appendix A) and (35) becomes: The relation n −xe 2r a r · O = a r n −x · O (proved, for example, in [22]) yields: (1−|z | 2 ) 2 , which shows that (z) does not depend upon the variable z, as announced.

Appendix E: Proof of Lemma 3.1.2
In this section we prove the following lemma. where erf is the error function defined as: Proof We consider the following double integrals: (36) so that: Since the variables are separable, we have: One can easily see that: We now give a simplified expression for i . We set f i (x) = e −x 2 /(2σ 2 i ) and then we have, because of Lemma 3.1.1: The change of variable x = arctanh(r) implies dx = dr 1−r 2 and yields: Because of the above definition of λ , this reduces to π ω 0 λ tanh(r) sinh(2r) dr.