Attractor and saddle node dynamics in heterogeneous neural fields

Background: We present analytical and numerical studies on the linear stability of spatially non-constant stationary states in heterogeneous neural fields for specific synaptic interaction kernels. Methods: The work shows the linear stabiliy analysis of stationary states and the implementation of a nonlinear heteroclinic orbit. Results: We find that the stationary state obeys the Hammerstein equation and that the neural field dynamics may obey a saddle-node bifurcation. Moreover our work takes up this finding and shows how to construct heteroclinic orbits built on a sequence of saddle nodes on multiple hierarchical levels on the basis of a Lotka-Volterra population dynamics. Conclusions: The work represents the basis for future implementation of meta-stable attractor dynamics observed experimentally in neural population activity, such as Local Field Potentials and EEG.


Background
Neural field models, such as [1,2] are continuum limits of large-scale neural networks.Typically, their dynamic variables describe either mean voltage [2] or mean firing rate [1,3] of a population element of neural tissue (see [4,5] for recent reviews).
The present article considers the paradigmatic Amari equation [2] describing the spatiotemporal dynamics of mean potential V (x, t) over a cortical d-dimensional manifold ⊂ R d : where K(x, y) is the spatial synaptic connectivity between site y ∈ and site x ∈ , and S is a nonlinear, typically sigmoidal, transfer function.This model neglects external inputs for simplicity but without constraining the generality of the subsequent results.Possible synaptic time scales are supposed to be included in the kernel function K and can be introduced by a simple scaling of time.
In general, the connectivity kernel K(x, y) fully depends on both sites x and y, which case is referred to as spatial heterogeneity.If the connectivity solely depends on the difference between x and y, i.e.K(x, y) = K(x − y), the kernel is called spatially homogeneous http://www.epjnonlinearbiomedphys.com/content/2/1/4 [2].Furthermore, if the connectivity depends on the distance between x and y only, i.e.K(x, y) = K(||x − y||), with ||x|| as some norm in , the kernel is spatially homogeneous and isotropic [6].
Spatially homogeneous (respectively isotropic) kernels have been intensively studied in the literature due to their nice analytical properties.In this case, the evolution equations have exact solutions such as bumps [2,7], breathers [8][9][10] or traveling waves [7,11].Moreover, such kernels allow the application of the technique of Green's functions for deriving partial neural wave equations [7,12,13].
The present work focuses on spatially heterogeneous neural fields which have been discussed to a much lesser extent in previous studies than homogeneous neural fields [14][15][16][17][18][19][20][21][22].This study resumes these attempts by investigating stationary states of the Amari equation (1) with heterogeneous kernels and their stability.Such a theory would be mandatory for modeling transient neurodynamics as is characteristic, e.g., for human cognitive phenomena [23], human early evoked potentials [24] or, among many other phenomena, for bird songs [25].
The article is structured in the following way.In the "Results" section we present new analytical results on stationary solutions of the Amari equation (1) and their stability in the presence of heterogeneous connectivity kernels.Moreover, we present numerical simulation results for the kernel construction and its stability analysis.The "Methods" section is devoted to construct such kernels through dyadic products of desired stationary states (cf. the previous work of Veltz and Faugeras [26]).A subsequent linear stability analysis reveals that these stationary solutions could be either attractors or saddles, depending on the chosen parametrization.Finally, we present a way to connect such saddle state solutions via heteroclinic sequences [27,28] in order to construct transient processes.

Results
In this section we present the main results of our study on heterogeneous neural fields.

Analytical study
The Amari equation (1) has a trivial solution V 0 (x) = 0 and non-trivial solutions V 0 (x) = 0 that obey the Hammerstein integral equation [29] V Inspired by Hebbian learning rules for the synaptic connectivity kernel K(x, y) which found successful applications, e.g., in bi-directional associative memory [30], we consider symmetric spatially heterogeneous kernels K(x, y) = K(y, x) that can be constructed from dyadic products of the system's non-trivial stationary states Together with Eq. ( 2), this choice yields the additional condition for non-trivial stationary states which is a nonlinear integral equation of Fredholm type.Since 0 < S(x) < 1 for a logistic transfer function, a necessary condition for non-trivial stationary states is which indicates immediately a method to find a non-trivial solution numerically as shown below.
Small deviations u(x, t) = V (x, t) − V 0 (x) from a non-trivial stationary state V 0 (x) obey the linear integro-differential equation where L(x, y) = K(x, y)S (y) and S (y) = dS[V 0 (y)] /dV 0 (y).A linear stability analysis of ( 6) carried out in the next section shows that a non-trivial stationary state V 0 (x) is either a fixed point attractor, or (neglecting a singular case) a saddle with one-dimensional unstable manifold.Such saddles can be connected to form stable heteroclinic sequences.

Numerical study
To gain deeper insight into possible stationary solutions of the Amari equation ( 1) and their stability, the subsequent section presents the numerical integration of equation ( 1) in one spatial dimension for a specific spatial synaptic connectivity kernel.Since previous experimental studies [31] have revealed Gaussian distributed probability densities of neuron interactions in the visual cortex of rats, it is reasonable to look for spatially discretized stationary states in the family of Gaussian functions parameterized by the amplitude W 0 , the variance σ 2 , the noise level κ and the spatial discretization interval x.By virtue of this parametrization of the discrete Hammerstein equation, it is sufficient to fit the model parameters optimally in such a way that the Hammerstein equation holds.Figure 1(a) illustrates a noisy kernel K and (b) shows the corresponding stationary state V 0 (x) for certain parameters.
For each stationary state, one obtains a kernel L(x, y) of the linear stability analysis whose spectrum characterizes the stability of the system in the vicinity of the stationary state.If the eigenvalue with the maximum real part 1 > 1, then the stationary state V 0 (x) is exponentially unstable whereas 1 < 1 guarantees exponential stability.Figure 1(b) shows the eigenmode e 1 (x) corresponding to the eigenvalue with maximum real part which has a similar shape as the stationary state.
Moreover Figure 2 presents parameters for which V 0 (x) fulfills the Hammerstein equation (2), i.e. for which V 0 (x) is the stationary solution.We observe that some parameter sets exhibit a change of stability, i.e. the eigenvalue with maximum real part may be 1 > 1 or 1 < 1 for certain parameter subsets.Taking a closer look at the stability of V 0 (x), the computation of the eigenvalues k , k = 1, . . ., n reveals a dramatic gap in the spectrum: the eigenvalue with maximum real part 1 is well isolated from the rest of the spectrum { k>1 } with | k>1 | < 10 −14 .This is in accordance to the discussion of Eq. ( 17) on the linear spectrum.In contrast, for the system starting close to an unstable stationary state, cf.point (2) in Figure 2, the field activity moves away from V 0 (x) and approaches a new stationary  Here the stationary state V 0 (x) is exponentially stable with the parameters of point (1) in Figure 2. The simulation starts close to the stationary state V(x, 0) = V 0 (x) + ξ(x) with the random numbers ξ(x) taken from a zero-mean Gaussian distribution with variance 0.3.The plot on the left hand side shows the deviation from the stationary state V 0 (x).The plot on the right hand side shows the stationary states V 0 (x) (solid line) and the final field activity at large times V(x, t = 125) (dashed line) which is almost identical to V 0 (x).The spatial domain has length L = 1 with n = 300.
state close to but different from V 0 (x), cf. Figure 4.This new stationary state obeys the Hammerstein equation (2).
Recalling the presence of the trivial stable solution V = 0, the activity shown in Figure 4 indicates the bistability of the system for the given parameter set.
Figure 5 supports this bistability for the same parameter set but different initial conditions, which presents the jump from the unstable stationary state V 0 (x) to the trivial stable stationary state V = 0.The choice whether the system approaches the upper or lower stable stationary state depends on the initial condition of the simulation and is random for random initial conditions as implemented in Figures 3, 4 and 5. Hence, this example reveals the existence of a saddle-node bifurcation in heterogeneous neural fields.
Finally, we would like to stress that the analysis presented above does not depend on the smoothness of the kernel and stationary state.For a strong noise level in the synaptic interaction kernel K, the analytical discussion above still describes the stationary state and the linear stability quite well as shown in Figure 6 for a stable stationary state V 0 (x) close to the stability threshold.

Heteroclinic orbits
The previous section has shown that heterogeneous neural fields may exhibit various stationary states with different stability properties.In particular we found that stationary states could be saddles with one-dimensional unstable manifolds that could be connected to stable heteroclinic sequences (SHS: [27,28]) which is supported by experimental Here the stationary state V 0 (x) is exponentially unstable with the parameters of point (2) in Figure 2. The plot on the left hand side shows the deviation from the stationary state V 0 (x), the right hand side panel presents the stationary states V 0 (x) (solid line) and the final field activity at large times V(x, t = 125) (dashed line) which is close to V 0 (x).Other parameters are taken from Figure 3. evidence [24,32,33].We present in the following paragraphs our main findings on heterogeneous neural fields exhibiting heteroclinic sequences and also hierarchies of such sequences.

One level heteroclinic sequence
It is possible to expand the integral in the Amari equation (1) into a power series yielding with kernels The solution of Eq. ( 8) represents a heteroclinic sequence that connects saddle points {V k (x)} along their respective stable and unstable manifolds.Its transient evolution is described as winnerless competition in a Lotka-Volterra population dynamics governed by interaction weights ρ ik between neural populations k and i and their respective growth rates σ i .
In Eq. ( 9) the {V + k (x)} comprise a bi-orthogonal system of the saddles {V k (x)}.Therefore, the kernel K 1 (x, y) describes a Hebbian synapse between sites y and x that has been trained with pattern sequence V k .This finding confirms the previous result of [18].Moreover the three-point kernel K 2 (x, y, z) further generalizes Hebbian learning to interactions between three sites x, y, z ∈ .Note that the kernels K i are linear combinations of dyadic product kernels, similar to those introduced in (3).Thus, our construction of heteroclinic orbits straightforwardly results in Pincherle-Goursat kernels used in [26].

Multi-level hierarchy of heteroclinic sequences
Now we assume that the neural field supports a hierarchy of stable heteroclinic sequences in the sense of [25,34].For the general case one has to construct integral kernels for a much wider class of neural field equations which can be written as where the new temporal kernel G describing synaptic-dendritic filtering is usually the Green's function of a linear differential operator Q, such that G(t, τ ) = G(t − τ ) and the temporal integration in Eq. ( 10) is temporal convolution.Equation ( 10) can be simplified by condensing space x and time t into spacetime s = (x, t).Then, (10) becomes with a tensor product kernel H(s, s For only two levels in such a hierarchy, we obtain with kernels σ (1) σ (2) k 1 (x) .http://www.epjnonlinearbiomedphys.com/content/2/1/4 Here, the V (ν) k ν (x) denote the k ν -th saddle in the ν-th level of the hierarchy (containing n ν stationary states).Saddles are chosen again in such a way that they form a system of bi-orthogonal modes V (ν) + k ν (x) whose behavior is determined by Lotka-Volterra dynamics with growth rates σ (ν) k ν > 0 and time-dependent interaction weights ρ (ν) k ν k ν (t) = 1 that are given by linear superpositions of templates r (ν) k ν j ν l ν+1 .Additionally, τ (ν) are the characteristic time scales for level ν.Levels are temporally well separated through τ (ν)  τ (ν+1) [35].Interestingly these kernels are time-independent.Since neural field equations can be written in the same form as Eq. ( 12), this result shows that hierarchies of Lotka-Volterra systems are included in the neural field description.Again we point out that the resulting kernels are linear combinations of dyadic products as introduced in Eq. ( 3), see also the work of Veltz and Faugeras [26].

Discussion
This study considers spatially heterogeneous neural fields describing the mean potential in neural populations according to the Amari equation [2].To our best knowledge this work is one of the first deriving the implicit conditions for stationary solutions and identifying the corresponding stability constraint as the Hammerstein integral equation.Moreover, as one of the first studies our work derives conditions for the linear stability of such stationary solutions and derives an analytical expression for stability subjected to the properties of heterogeneous synaptic interaction.The analytical results are complemented by numerical simulations illustrating the stability and instability of heterogeneous stationary states.We point out that the results obtained extend previous studies both on homogeneous neural fields and heterogeneous neural fields as other studies in this context have done before [21].
By virtue of the heterogeneity of the model, it is possible to consider more complex spatiotemporal dynamical behavior than the dynamics close to a single stationary state.We show how to construct hierarchical heteroclinic sequences built of multiple stationary states each exhibiting saddle node dynamics.The work demonstrates in detail how to construct such sequences given the stationary states and their saddle node dynamics involved.Motivated by a previous study on hierarchical heteroclinic sequences [25,34], we constructed such sequences in heterogeneous neural fields of the Amari type.Our results indicate that such a hierarchy may be present in a single heterogeneous neural field whereas previous studies [25,34] consider the presence of several neural populations to describe heteroclinic sequences.The kernels obtained from heteroclinic saddle node dynamics are linear superpositions of tensor product kernels, known as Pincherle-Goursat kernels in the literature [26].

Conclusion
The present work is strongly related to the literature hypothesizing the presence of chaotic itinerant neural activity, cf.previous work by [36,37].This concept of chaotic itinerancy is attractive but yet lacking a realistic neural model.We admit that the present work represents just a first initial starting point for further model analysis of sequential neural activity in heterogeneous neural systems.It opens up an avenue of future research and promises to close the gap between the rather abstract concept of sequential, i.e. temporally transient, neural activity and corresponding mathematical neural models.http://www.epjnonlinearbiomedphys.com/content/2/1/4

Stationary states and their stability
In order to learn more about the spatiotemporal dynamics of a neural field, in general it is reasonable to determine stationary states and to study their linear stability.This already gives some insight into the nonlinear dynamics of the system.Since the dynamics depends strongly on the interaction kernel and the corresponding stationary states, it is necessary to work out conditions for the existence and uniqueness of stationary states and their stability.

Analytical study
For stationary solutions, the left hand side of Eq. ( 1) vanishes and we obtain the Hammerstein equation ( 2) [29].It has a trivial solution V 0 (x) = 0 and further non-trivial solutions V 0 (x) = 0 under certain conditions.The existence and number of solutions of Eq. ( 2) depends mainly on the operator A = I + KF [38] and its monotonicity [39], with the unity operator I, the linear integral operator Ku = K(x, y)u(y)dy and the Nemytskij operator Fu(x) = S(u(x)).For instance, a simple criterion for the existence of at least a single non-trivial solution is the symmetry and positive definiteness of the kernel K(x, y) and the condition S(u) ≤ C 1 u + C 2 [29] .Moreover previous studies have proposed analytical [40] and numerical [41] methods to find solutions of the Hammerstein equation ( 2).
To illustrate the non-uniqueness of solutions of the Hammerstein equation, let us expand the stationary solution into a set of bi-orthogonal spatial modes {φ n (x)}, {ψ n (x)} with constant coefficients V n , U n and with the Kronecker symbol δ n,m .Then the Hammerstein equation recasts to with Since f m is a nonlinear function of V n , Eq. ( 13) has multiple solutions {V n } for a given bi-orthogonal basis.Hence, V 0 (x) may not be unique.
Considering small deviations u(x, t) = V (x, t) − V 0 (x) from a non-trivial stationary state V 0 (x), these deviations obey the linear equation ( 6) with L(x, y) = K(x, y)S (y) and S (y) = dS[V 0 (y)] /dV 0 (y).http://www.epjnonlinearbiomedphys.com/content/2/1/4 A solution of ( 6) is then u(x, t) = exp(λt)e(x), λ ∈ C with mode e(x).Inserting this solution into Eq.( 6) yields the continuous spectrum of the corresponding linear operator L(x, y) determined implicitly by the eigenvalue equation Under the above assumption of a dyadic product kernel (3), the eigenfunctions {e k (x)} of the kernel L(x, y), defined through with eigenvalues k = λ k + 1 ∈ C form an orthonormal system with respect to the scalar product with weight S (y).This follows from the dyadic product kernel (3) through From (17) we deduce two important results: • a certain eigenmode e k 0 (x) is proportional to the stationary state V 0 (x) with scaling factor V 0 , e k 0 S / k 0 , and • other eigenmodes in the eigenbasis e k =k 0 (x) are orthogonal to V 0 (x) yielding Hence for dyadic product kernels (3) the spectrum of L includes one eigenmode with eigenvalue 1 = 0, i.e. λ 1 = −1, while all other eigenmodes are stable with k =1 = 0 and thus λ k =1 = −1.Therefore, a non-trivial stationary state could become either an asymptotically stable fixed point, i.e. an attractor, for n < 1, for all n, or a saddle with a one-dimensional unstable manifold, for 1 > 1 and n =1 < 1 (neglecting the singular case 1 = 1).Finally we have to justify the necessary condition (5) which follows from 0 < S(x) < 1 and inserted into the Hammerstein equation ( 2) for a dyadic product kernel (3).

Numerical study
In order to investigate stationary states of the Amari equation numerically, we choose a finite one-dimensional spatial domain of length L and discretize it into a regular grid of n intervals with grid interval length x = L/n.Then the kernel function is K(x = n x, y = m x) = K nm / x and Eq. ( 1) reads where V n (t) = V (n x, t).The corresponding Hammerstein equation ( 2) is given by Taking into account the insight from the discussion of Eq. ( 3) and its consequences for the eigenmodes, we have chosen the spatial kernel to We employ an Euler-forward integration scheme for the temporal evolution with discrete time step t = 0.05.
We render the stationary state random by adding the noise term η i which are random numbers taken from a normal distribution with zero mean and unit variance.We point out that we choose η i such like | n i=1 η i | < 0.05 in order to not permit an amplitude increase in the dynamics by noise, but just as a modulation.The sigmoid function is chosen as S(V ) = 1/(1 + exp(−α(V − θ))) parameterized by the slope parameter α and the mean threshold θ.

Heteroclinic orbits
The general neural field equation (10) supplies the Amari equation (1) as a special case for G(t) = e −t (t) ( (t) as Heaviside's step function).Then Q = ∂ t + 1 is the Amari operator for the Green's function G(t).For second-order synaptic dynamics [42] and for the filter properties of complex dendritc trees [43], more complicated kernels or differential operators, respectively, have to be taken into account.
As a first step for constructing heteroclinic sequences for (1) or (10), we expand the nonlinear transfer function S in Eq. ( 11) into a power series about a certain state V (s) [15], Inserting V (s) from Eq. ( 11) yields By eliminating the nonlinearities S[V (s j )] using the power series above again, we get Inserting this expression into Eq.( 11) leads to a generalized Volterra series with a sequence of integral kernels H m that can be read off after some further tedious rearrangements [44].

One-level hierarchy
In a previous work, we have derived a one-level hierarchy of stable heterogenic sequences [44], which we briefly recapitulate here.We assume a family of n stationary states V k (x), 1 ≤ k ≤ n that are to be connected along a heteroclinic sequence.Each state is assumed to be realized by a population in the neural field governed by (10), that is characterized by a population activity α k (t) ∈[ 0, 1].Then the overall field quantity is obtained through an order parameter expansion [15,45] V These population amplitudes result from a winnerless competition in a generalized Lotka-Volterra system approximating a Wilson-Cowan model [46,47] with growth rates σ k > 0, and interaction weights ρ kj > 0, ρ kk = 1, that are trained by the algorithm of [27] and [28] for the desired sequence of transitions.
For the following construction, we also assume that the system of stationary states {V k (x)} is linearly independent such that there is a bi-orthogonal system {V + k (x)}, obeying Then, we obtain from ( 22) Next, we take the derivation of ( 22) with respect to time t, by exploiting ( 23) Adding (22), we obtain from which we eliminate all occurrences of ξ by means of ( 25).This yields Moreover, the nonlinear spatial integral transformation in the Amari equation ( 1) may be written as a generalized Volterra series ( 21) Comparison of the first three terms with (28) yields the result (9).

Figure 3
Figure3presents the spatiotemporal evolution of the heterogeneous neural field starting close to the stable stationary state V 0 (x), see point (1) in Figure2.As expected, the field activity remains in the vicinity of the stable state.In contrast, for the system starting close to an unstable stationary state, cf.point (2) in Figure2, the field activity moves away from V 0 (x) and approaches a new stationary

Figure 2
Figure 2 Parameters which guarantee the solution of the Hammerstein equation (2) which are the stationary solutions of Eq. (1) and the line styles encode their exponential stability (stable: solid line, unstable: dashed line).Further parameters are κ = 0.0, θ = 3.0.

Figure 3
Figure 3 Numerically simulated spatio-temporal dynamics of the heterogeneous neural field (left) and stationary states (right).Here the stationary state V 0 (x) is exponentially stable with the parameters of point (1) in Figure2.The simulation starts close to the stationary state V(x, 0) = V 0 (x) + ξ(x) with the random numbers ξ(x) taken from a zero-mean Gaussian distribution with variance 0.3.The plot on the left hand side shows the deviation from the stationary state V 0 (x).The plot on the right hand side shows the stationary states V 0 (x) (solid line) and the final field activity at large times V(x, t = 125) (dashed line) which is almost identical to V 0 (x).The spatial domain has length L = 1 with n = 300.

Figure 4
Figure 4 Numerically simulated spatio-temporal dynamics of the heterogeneous neural field (left) and stationary states (right).Here the stationary state V 0 (x) is exponentially unstable with the parameters of point (2) in Figure2.The plot on the left hand side shows the deviation from the stationary state V 0 (x), the right hand side panel presents the stationary states V 0 (x) (solid line) and the final field activity at large times V(x, t = 125) (dashed line) which is close to V 0 (x).Other parameters are taken from Figure3.

Figure 5 4 Figure 6
Figure 5 Numerically simulated spatio-temporal dynamics of the heterogeneous neural field (left) and stationary states (right).Here the stationary state V 0 (x) is exponentially unstable with the parameters of point (2) in Figure2.All parameters are identical to Figure3, however the final field activity at large times V(x, t = 125) (dashed line) is close to the trivial solution V = 0.This jump of activity recalls the presence of the trivial stationary solution V = 0 and hence reflects the multi-stability of the heterogeneous neural field.