A From a double-ring to a single layer velocity integration network

Velocity integration in a multilayer neural field model of spatial working memory1

Abstract

We analyze a multilayer neural field model of spatial working memory, focusing on the impact of interlaminar connectivity, spatial heterogeneity, and velocity inputs. Models of spatial working memory typically employ networks that generate persistent activity via a combination of local excitation and lateral inhibition. Our model is comprised of a multilayer set of equations that describes connectivity between neurons in the same and different layers using an integral term. The kernel of this integral term then captures the impact of different interlaminar connection strengths, spatial heterogeneity, and velocity input. We begin our analysis by focusing on how interlaminar connectivity shapes the form and stability of (persistent) bump attractor solutions to the model. Subsequently, we derive a low-dimensional approximation that describes how spatial heterogeneity, velocity input, and noise combine to determine the position of bump solutions. The main impact of spatial heterogeneity is to break the translation symmetry of the network, so bumps prefer to reside at one of a finite number of local attractors in the domain. With the reduced model in hand, we can then approximate the dynamics of the bump position using a continuous time Markov chain model that describes bump motion between local attractors. While heterogeneity reduces the effective diffusion of the bumps, it also disrupts the processing of velocity inputs by slowing the velocity-induced propagation of bumps. However, we demonstrate that noise can play a constructive role by promoting bump motion transitions, restoring a mean bump velocity that is close to the input velocity.

m

ultilayer networks, neural fields, stochastic differential equations, bump attractors

{AMS}

68Q25, 68R10, 68U05

1 Introduction

Spatial working memory tasks test the brain’s ability to encode information for short periods of time [30, 54]. A subject’s performance during such tasks can be paired with brain recordings to help determine how neural activity patterns represent memory during a trial [32]. In general, working memory involves the retention of information for time periods lasting a few seconds [6]. More specifically, spatial working memory involves the short term storage of a spatial variable, such as idiothetic location [16] or a location on a visual display [24]. A well tested theory of spatial information storage on short timescales involves the generation of persistent activity that encodes input during the retention interval [74]. Network models of this activity typically involve local excitation and broader inhibition, producing localized activity packets referred to as bump attractors [19, 45]. These models have recently been validated using recordings from oculomotor delayed-response tasks in monkeys [77] and from grid cell networks of freely moving rats [79]. This suggests that studying network mechanisms for generating reliable neural activity dynamics can provide insight into how the brain robustly performs spatial working memory tasks.

In addition to the short term storage of location, several networks of the brain can integrate velocity signals to update a remembered position [52]. Angular velocity of the head is used by the vestibular system to update memory of heading direction [72]. Furthermore, intracellular recordings from goldfish demonstrate that eye position can be tracked by neural circuits that integrate saccade velocity to update memory of eye orientation [1]. Velocity integration has also been identified in place cell and grid cell networks, which track an animal’s idiothetic location [76, 34, 31]. While these networks each possess distinct circuit mechanisms for integrating and storing information, the general dynamics of their stored position variables tends to be similar [51]. Neuronal networks that support a continuous (line) or approximately continuous (chain) attractor of solutions constitute a unifying framework for modeling these different systems [43]. One can then consider the effect of noise, network architecture, or erroneous inputs on the accuracy of position memory [81, 62, 15].

One important feature of spatial working memory, often overlooked in models, is its distributed nature [37]. Most models focus on the dynamics of persistent activity representing position memory in a single-layer network [81, 19, 47]. However, extensive evidence demonstrates working memory for visuo-spatial and idiothetic position is represented in several distinct modules in the brain that communicate via long-range connectivity [23, 67]. There are many possible advantages conferred by such a modular organization of networks underlying spatial memory. One well tested theory notes different network layers can represent position memory on different spatial scales, leading to higher accuracy within small-scale layers and wider range in large-scale layers [14]. Furthermore, the information contained in spatial working memory is often needed to plan motor commands, so it is helpful to distribute this signal across sensory, memory, and motor-control systems [64]. Another advantage of generating multiple representations of position memory is that it can stabilize the memory through redundancy [69]. For instance, coupling between multiple layers of a working memory network can reduce the effects of noise perturbations, as we have shown in previous work [39].

In addition to being distributed, the networks that generate persistent activity underlying spatial working memory also appear to be heterogeneous. For instance, prefrontal cortical networks possess a high degree of variation in their synaptic plasticity properties as well as their cortical wiring [60, 75]. Furthermore, there is heterogeneity in the way place cells from different hippocampal regions respond to changes in environmental cues [3, 48]. Along with such between-region variability, there is local variability in the sequenced reactivations of place cells that mimic the activity patterns that typically occur during active exploration [55]. In particular, these reactivations are saltatory, rather than smoothly continuous, so activity focuses at a discrete location in the network before rapidly transitioning to a discontiguous location. Such activity suggests that the underlying network supports a series of discrete attractors, rather than a continuous attractor [12].

Given the spatially distributed and heterogeneous nature of neural circuits encoding spatial working memory, we will analyze tractable models that incorporate these features. We are particularly interested in how the architecture of multilayer networks impacts the quality of the encoded spatial memory. In previous work, we examined networks whose interlaminar connectivity was weak and/or symmetric, ignoring the effects of spatial heterogeneity in constituent layers [39, 40]. In this work, we will depart from the limit of weak coupling, and derive effective equations for the dynamics of bumps whose positions encode a remembered location. Through the use of linearization and perturbation theory, we can thus determine how both the spatial heterogeneity of individual layers and the coupling between layers impact spatial memory storage. In previous work, we found that spatial heterogeneity can help to stabilize memories of a stationary position [42], but such heterogeneities also disrupt the integration of velocity inputs [58]. Thus, it is important to understand the advantages and drawbacks of heterogeneities, and quantify how they trade off with one another.

We focus on a multilayer neural field model of spatial working memory, with arbitrary coupling between layers and spatial heterogeneity within layers. Furthermore, as we are interested in both the retention of memory and the integration of input, we incorporate a velocity-based modulation to the recurrent connectivity which is non-zero when the network receives a velocity signal [81]. The stationary bump solutions of this network are analyzed in Section 3. Since the effects of velocity input and heterogeneity are presumed to be weak, the stationary bump solutions only depend upon the connectivity between layers. Analyzing the stability of bumps, we can determine the marginally stable modes of these bump solutions which will be susceptible to noise perturbations. Subsequently, we derive a one-dimensional stochastic equation that describes the response of the bump solutions to noise, velocity input, and spatial heterogeneity. With this approximation in hand, we can determine the effective diffusion and velocity of bumps using asymptotic methods, which compare well with numerical simulations of the full model (Section 4). Lastly, we analyze more nuanced architectures in Section 5, whose bump solutions possess multiple marginally stable modes. As a result, we find we must derive multi-dimensional stochastic equations to describe their dynamics in response to noise. Our work examines in detail the effects of modular network architecture on the coding of spatial working memory.

2 Multilayer neural field with spatial heterogeneity

Figure 1: Schematic of multilayer network features. (A) Purely excitatory long-range interlaminar connectivity is activated by regions of high activity such as the bump attractor solution in layer 1 (), projecting to similarly tuned locations in layer 2, reinforcing the position of the activity bump there (). (B) Different network topologies as specified by the weight functions (, ) are explored in two layer networks (feedforward, symmetric, and asymmetric) as well as three layers (directed ring, inward star). (C) Local heterogeneities within each layer introduced into the recurrent weight functions , Eq. (1), generate preferred locations for the bump attractor solutions to the model Eq. (2). We consider a variety of networks, which possess different attractor structures in each of their constituent layers. Continuous attractors possess marginally stable bump solutions at each location around the ring, while chains of discrete attractors possess stable nodes (blue dots) where bumps prefer to reside separated by saddles (red circles). (D) Velocity integration via the asymmetric integral term involving in Eq. (1) causes bump attractor solutions to move about the domain, incrementing position in proportion to the velocity amplitude.

Neural field models of persistent activity have been used extensively to understand the relationships between network properties and spatiotemporal activity dynamics [21, 10]. Stable bump attractors arise as solutions to these models when network connectivity is locally excitatory and broadly inhibitory, and these solutions are translationally invariant when connectivity is also strictly distance-dependent [2, 25]. However, the incorporation of multiple neural field layers and spatial heterogeneity can break the translation invariance of single network layers, so that bumps have preferred positions within their respective layer [27, 29, 41, 42]. Our analysis focuses on a multilayer neural field model with general connectivity between layers. Spatial heterogeneity within layers, velocity input, and noise are all assumed to be weak ():

(1)

where denotes the average neural synaptic input at location at time in network layer , and is a convolution. Note that we have restricted the spatial domain to be one-dimensional and periodic. There are several experimental examples of spatial working memory which operate on such a domain including oculomotor delayed-response tasks for visual memory [30, 77] as well as spatial navigation along linear tracks [7, 80]. While we suspect that several of our findings extend to two-dimensional spatial domains [57], we reserve such analysis for future work. Recurrent synaptic connectivity within layers is given by the collection of kernels , and we allow these functions to be spatially heterogeneous, rather than simply distance-dependent. We thus define them as

(2)

where the impact of the heterogeneity is weak (), and is only dependent on the distance . As opposed to recurrent connectivity, we assume the interlaminar connectivity (, ) is homogeneous, so we can always write . The homogeneous portion of the recurrent connectivity in each layer is locally excitatory and laterally inhibitory: e.g., the unimodal cosine function

(3)

which we use in some of our computations. Similarly, we will often consider a cosine shaped excitatory weight function for interlaminar connectivity:

(4)

We introduce homogeneous, distance-dependent kernels for the connectivity between layers. This is motivated by recent experimental work demonstrating that several brain areas involved in spatial working memory are reciprocally coupled to one another [20, 23], and these areas all tend to have similar topographically organized delay period activity [68, 38, 30]. Thus, we expect that topologically organized connectivity would be re-enforced via Hebbian plasticity rules [44, 59]. Such connectivity functions tend to generate stationary bump solutions within each layer [36, 46, 41, 39], and we will analyze these solutions in some detail in Section 3. Other lateral inhibitory functions, such as sums of multiple cosine modes, will also generate stationary bump solutions but they do not qualitatively alter the dynamics of the system.

Note, the general form of the weight functions allows us to explore a variety of network topologies, and their impact on the dynamics of bump attractors. For example, it is clear that a feedforward network (Fig. 1A) will primarily be governed by the dynamics of the upstream layer. However, the dynamics of bumps in more intricate networks (Fig. 1B) are more nuanced. Applying both linear stability analysis and perturbation theory to bumps in Section 3, we can explore the specific impacts of different conformations of . Furthermore, we expect the heterogeneities arising in local connectivity Eq. (2) will interact with interlaminar connectivity to shape the overall dynamics of bumps (Fig. 1C).

The impact of neural activity via synaptic connectivity is thus given via the integral terms, where a nonlinearity is applied to the synaptic input variables:

and such sigmoids are analogous to the types of saturating nonlinearities that arise from mean field analyses of spiking population models [13, 61]. For analytical tractability, we often consider the high gain limit () in our examples, resulting in the Heaviside nonlinearity [2, 21]

(5)

The effects of velocity inputs are accounted for by the second integral term in Eq. (1), based on a well tested model of the head direction system [72] as well as spatial navigation models that implement path integration [65, 51]. While some of these models use multiple layers to account for different velocity directions [78, 15], the essential dynamics are captured by a single-layer with recurrent connections modulated by velocity input [81, 58]. Since we are studying motion along a one-dimensional space, the weak () velocity input to each neural activity layer is given by a scalar variable which can be positive (for rightward motion) or negative (for leftward motion) as shown in Fig. 1D. We derive a reduction of the double ring model (one ring for each velocity direction) of velocity integration presented in [78] to a single layer for velocity (positive or negative) in the Appendices. The connectivity functions targeting each layer should be interpreted as interactions that are shaped by an incoming velocity signal to that layer. Essentially, this connectivity introduces asymmetry into the weight functions, which will cause shifts in the position of spatiotemporal solutions. Typically, this weight function is chosen to be of the form , in single layers [81]. In the absence of any heterogeneity, such a layer will have bumps that propagate at velocity precisely equal to  [58]. As shown in the Appendices A and B, we can extend this previous assumption to incorporate velocity-related connectivity that respects the interlaminar structure of the network, so that

(6)

As we demonstrate in Section 3.3, this results in bump solutions that propagate with velocity .

Dynamic fluctuations are a central feature of neural activity, and they can often serve to corrupt task pertinent signals, creating error in cognitive tasks [26]. The error in spatial working memory tasks tends to build steadily in time, in ways that suggest the process underlying the memory may evolve according to a continuous time random walk [77, 8]. As there is no evidence of long timescale correlations in the underlying noise process, we are satisfied to model fluctuations in our model using a spatially correlated white noise process:

where is the spatial filter of the noise in layer and is a spatially and temporally white noise increment. We define the mean and covariance of the vector :

(7)

where is the even symmetric spatial correlation term, and is the Dirac delta function.

Subsequently, we will analyze the existence and stability of stationary bump solutions to Eq. (1) in Section 3.1. Since we will perform this analysis under the assumption of spatially homogeneous synaptic weight functions ( in Eq. (2)), these solutions will be marginally stable to perturbations that shift their position. However, once we incorporate noise, heterogeneity, and velocity inputs in Section 3.3, we can perturbatively analyze their effects by linearizing about the stationary bump solutions. The low-dimensional stochastic system we derive will allow us to study the impact of multilayer architecture on the processing of velocity inputs in Section 4.

3 Bump attractors in a multilayer neural field

Our analysis begins by constructing stationary bump solutions to Eq. (1) for an arbitrary number of layers and even, translationally-symmetric synaptic weight functions . Note, there are a few recent studies that have examined the existence and stability of stationary bump solutions to multilayer neural fields [27, 29, 39]. In particular, Folias and Ermentrout studied bifurcations of stationary bumps in a pair of lateral inhibitory neural field equations [29]. They identified solutions in which bumps occupied the same location in each layer (syntopic) as well as different locations (allotopic), and they also demonstrated traveling bumps and oscillatory bumps that emerged from these solutions. However, they did not study the general problem of an arbitrary number of layers, and their analysis of networks with asymmetric coupling was relatively limited. Since the solutions will form the basis of our subsequent perturbation analysis of heterogeneity and noise, we will outline the existence and stability analysis of bumps first, for an arbitrary number of layers . The reader is advised to consult the works of Folias and Ermentrout for a more detailed characterization of the possible bifurcations of stationary patterns in a pair of neural fields [27, 29]. We also note that, while we are restricting our analysis to the case of one-dimensional domains, we expect our results to extend to two or more dimensions as demonstrated in [57]. Furthermore, previous experiments in rats have probed the behavior and neurophysiological underpinnings of spatial navigation along linear tracks [7, 80]. Thus, we believe the model we analyze here would be pertinent to these cases in which the environment is nearly one-dimensional. After we characterize the stability of stationary bump solutions, we will consider the effects of weak perturbations to these solutions, which will help reveal how noise, heterogeneity, and interlaminar coupling shape the network’s processing of velocity inputs.

3.1 Existence of bump solutions

In the absence of a velocity signal and heterogeneity (), we can characterize stationary solutions to Eq. (1), given by . Conditions for the existence of stable stationary bumps in single layer neural fields have been well-characterized [2, 47, 33, 10], but much remains in terms of understanding how the form of would impact the existence and stability of bumps in a multilayer network. Furthermore, the stationary equations for bump solutions are a form of the well-studied Hammerstein equation [35, 4], and bump stability is characterized by Fredholm integral equations of the second kind [5]. For our purposes, we will construct bumps under the assumption that they exist. Then, we will employ self-consistency, to determine solution validity. This is straightforward in the case of a Heaviside nonlinearity , Eq. (5), but we can derive some results for general nonlinearities . First, note that, in the case of translationally symmetric kernels , we obtain the following convolution relating stationary solutions in each layer to one another:

(8)

In later analysis, we will also find the formula for the spatial derivative useful:

(9)

Next, since each must be periodic in , we can expand it in a Fourier series

(10)

where is the maximal integer index of a mode for bumps in all layers. Indeed, there will be a finite number of terms in the Fourier series, Eq. (10), under the assumption that the weight functions all have a finite Fourier expansion. Since most typical smooth weight functions are well approximated by a few terms in a Fourier series [73], we take this assumption to be reasonable. Once we do so, we can construct solvable systems for the coefficients of the bumps, Eq. (10), and their stability as in [17]. For even symmetric weight kernels, we can write

so that Eq. (8) implies that

(11a)
(11b)

Since the noise-free, heterogeneity-free system is translationally invariant, there is a family of solutions with center of mass at any location on . Furthermore, the evenness of the weight functions we have chosen implies the resulting system is reflection symmetric, so we can restrict our examination to even solutions, so for all , so Eq. (10) becomes

(12)

Plugging the formula Eq. (12) into Eq. (11), we find

(13)

The coefficients can be found using numerical root finders [73]. However, for particular functions and , we can project the system Eq. (13) to a much lower-dimensional set of equations, which can sometimes be solved analytically.

For instance, consider the Heaviside nonlinearity , Eq. (5). In this case, stationary bump solutions centered at are assumed to have superthreshold activity on the interval in each layer ; i.e. for . Applying this assumption to the stationary Eq. (8) yields

Self-consistency then requires that , as originally pointed out by Amari [2], which allows us to write

(14)

Again, Eq. (14) is a system of nonlinear equations, which can be solved numerically via root-finding algorithms. However, as opposed to the integral terms in Eq. (13), the integrals in Eq. (14) are tractable, which makes for a more straightforward implementation of a root-finder. If we utilize the canonical cosine weight functions, Eq. (3) and (4), we find we can carry out the integrals in Eq. (14) to yield:

(15)

Henceforth, we mostly deal with the specific case of cosine weight connectivity, although we suspect our results extend to the case of more general weight functions. This allows us to define connectivity simply using the scalar strength values of the interlaminar coupling, which comprise the off-diagonal entries of the following matrix: for . As discussed in Section 2, and specifically Fig. 1B, we categorize the network graphs of primary interest to our work here into the main cases of a two-layer network and some specific cases of a network with more layers. We now briefly demonstrate how such graph structures can impact the stationary solutions, as it foreshadows the impact on the non-equilibrium dynamics of the network.

Figure 2: Bump half-width plots for two-layer () networks with Heaviside nonlinearity , Eq. (5), and cosine coupling functions Eq. (3) and (4), as given by Eq. (15). (A,B,C) Bifurcation diagrams for half-width of bumps in the red layer shown in network diagram above. (A) Half-width of the bump in layer in a symmetric network, plotted as a function of threshold , as given by Eq. (17). Stable (solid) and unstable (dashed) branches of double bumps annihilate in a saddle-node (SN) bifurcation at low threshold , and at high threshold . Stable and unstable branches of single bumps also emerge from a SN for sufficiently high , while the stable branch annihilates with a branch of double bumps for low enough . (B) Half-width as a function of in an asymmetric network, as given by Eq. (17). (C) Bump half-width in a feedforward network, given by the single Eq. (18), shows both single bumps and double bump branches annihilate at the same upper threshold . (D) Critical coupling and input strength for needed to instantiate a single bump in layer 1 or a double bump solution, in a feedforward network, where . Shaded regions are generated by numerically simulating Eq. (1), and thick blue lines are calculated theoretically (see ‘Critical input needed for activation of bumps’ in main text). (E) Half-width of the layer 1 bump of a double bump solution for a recurrent network with over a range of coupling strength and threshold . Partitions demonstrate that a stable 1-bump solution also coexists in a subregion of the domain. No 2-bumps exist in the white region. (F) Half-width of the layer 2 bump of a double bump solution for a feedforward network over a range of coupling strength and threshold . Stable 1-bumps exist below the magenta line. For sufficiently large coupling and low threshold, only the ‘all-on’ solution exists in layer 2.

Two-layer symmetric network (). In this case, we can derive a few analytical results concerning the bifurcation structure of stationary bump solutions. However, to identify the half-widths and , it is typically necessary to solve Eq. (15) numerically to produce the plots shown in Fig. 2A. First of all, for double bump solutions, in which both layers possess stationary superthreshold activity, if we assume symmetric solutions, so that , then we can write Eq. (15) as

(16)

We cannot solve the transcendental Eq. (16) explicitly for the bump half-width . In order to gain some insight, we can identify the range over which solutions to the equations exist. This can be determined explicitly by finding the turning points of the right hand side of Eq. (16) (See blue dots in Fig. 2A,B,C), corresponding to the extrema of the function between which solutions exist. Thus, we can determine the location of these turning points, which are saddle-node (SN) bifurcations, by differentiating the right hand side :

so by requiring , we have

matching the locations of the double bump SN bifurcations (blue dots) shown in Fig. 2A.

Furthermore, SN bifurcations associated with the coalescing of stable single bump branches with unstable double bump branches (purple dots in Fig. 2A,B,C) can be determined using a threshold condition. For instance, given a layer 1 bump with half-width , we require the stationary solution in layer 2 () remains subthreshold (, ). Given in Eq. (15), single bump solutions in layer 1 satisfy , so are solutions with corresponding to the stable bump [41]. Thus, we require , so selecting for the maximal value of and plugging in , we have an explicit equation for the critical interlaminar strength above which there are no stable single bump solutions: , providing an implicit equation for the SN locations in Fig. 2A,B,C, and corresponding to the magenta curves in Fig. 2E,F.

Two-layer asymmetric network (). Double bump solution half-widths tend to differ in this case , obeying the pair of implicit equations

(17a)
(17b)

which we solve numerically to generate the branches plotted in Fig. 2B, as well as the surface plot in Fig. 2E. Note, however, it is still possible to determine the range of values in which stable single bump solutions exist in layer using the requirement , as derived in the symmetric network case.

Two-layer feedforward network (). This is a special case of the asymmetric network, where the nonlinear system, Eq. (15), defining the bump half-widths reduces to:

which can further be reduced to a single implicit equation for the half-width in the target layer 2 (see schematic in Fig. 2C):

(18)

which can be solved using numerical root finding to yield the curves in Fig. 2C,F.

‘All-on’ solutions in the two-layer network. Given excitatory interlaminar connections, it is possible to generate ‘all-on’ solutions in one and sometimes two layers of the network. An ‘all-on’ solution is one in which a layer has a stationary solution that is entirely superthreshold, for all . In the case of a feedforward network (Fig. 2C,F), the target layer 2 will have an ‘all-on’ solution when the minimal value of given a stable bump solution in layer 1. As a result, an ‘all-on’ solution in layer 2 would have the form

so requiring yields

obtaining equality along the grey line plotted in Fig. 2F. For recurrent networks, we can easily identify the threshold curves above which double ‘all-on’ solutions exist. These have the simpler forms:

so we need to require that and .

Critical input needed for activation of bumps. We are studying multilayer networks wherein we assume bump solutions can be instantiated by an external input. However, it is important to identify the critical input needed to nucleate and maintain such bumps in the two layers of the network. As demonstrated in Fig. 2A,B,C, there are multiple stable stationary solutions across a range of threshold and coupling values .

We wish to demonstrate that it is possible to instantiate a two bump solution given only an input, , to layer 1, and we focus exclusively on the feedforward network. This single layer will only have subthreshold activity if everywhere. If input is superthreshold (), stationary bump solutions, driven by an input in layer 1, are then given [36, 28]: . Thus, bumps driven by inputs just beyond the critical level , will have half-widths approximately satisfying . These bumps will have half-widths then given by the implicit equation . For values of with , we can show that this function is monotone decreasing, since when . In this case, . Therefore, as will tend to increase as is decreased from , so for we expect a single stable stationary bump solution in layer 1 (See also [28, 41]). This suggests either single or double bump solutions will emerge as long as , as shown in Fig. 2D. Increasing the strength of input will only serve to further stabilize this stationary bump. To determine the critical strength needed to propagate this bump forward to layer 2, we must solve for the half-width in , and require that layer 2 is driven superthreshold, so that . This admits an explicit inequality , so we need only solve for numerically to obtain the vertical boundary between single and double bump solutions in Fig. 2D.

3.2 Linear stability of bumps

Linear stability of the bump solutions , Eq. (8), can be determined by analyzing the evolution of small, smooth, separable perturbations such that . We expect () to be neutrally stable to translating perturbations , arising from the translation symmetry of Eq. (1) given . On the other hand, the bump may be linearly stable or unstable to perturbations of its half-width  [2]. The results we derive here for such perturbations are what determine the stability of branches plotted in Fig. 2.

To begin, consider , where thus describes perturbations to the shape of the bump that may evolve temporally. Plugging this into the full neural field Eq. (1) with , , and , we can apply the stationary Eq. (8), and subsequently write the equation as

(19)

Due to the linearity of the equation, we may apply separation of variables to each , such that  [66, 73]. Substituting into Eq. (19), we have for each :

(20)

Thus, each side of Eq. (20) depends exclusively on a different variable, or , so both must equal a constant . Therefore, for all , suggesting perturbations will grow indefinitely as for , indicating instability. While oscillatory instabilities are plausible ( with ), given specific forms of interlaminar coupling (e.g., combinations of interlaminar excitation and inhibition [29]), we did not identify such instabilities in the mutual excitatory layered networks we studied (Fig. 2). Thus, we expect instabilities emerging where will typically be of saddle-node type (). Furthermore, the equation for is now given for all , as

(21)

Eigenvalues are thus determined by consistent solutions for , to Eq. (21). One such solution is for , as can be shown by applying Eq. (9). As mentioned above, this demonstrates the neutral stability of bump solutions to translating perturbations, due to the translational invariance of Eq. (1).

Further analysis in the case of a general firing rate function can be difficult. However, if we consider the Heaviside nonlinearity given by Eq. (5), we obtain a specific case of Eq. (21), which is easier to analyze [2, 22]:

(22)

where we have made use of the fact that

(23)

since , and we have assigned

(24)

under the assumption that is monotone decreasing in , which is the case for cosine weight functions, Eq. (3) and (4). Note, all eigenfunctions of Eq. (22) that satisfy the conditions , for all , have associated eigenvalue given by so , which does not contribute to any instabilities. To specify other eigensolutions, we examine cases where for at least one . In such cases, we can obtain expressions for the eigenvalues by examining Eq. (22) at the points . In this case, the eigenfunctions are defined by their values at the threshold crossing points: for each . Thus, defining these unknown values , we simplify Eq. (22) to a linear system of equations of the form

(25)

where the elements of the blocks of the matrix are . We make use of the result  [71], which implies the set of eigenvalues of is the union of the set of the eigenvalues of and . Subsequently, the eigenvalues of will be . We now outline a few examples in which we can compute these eigenvalues analytically.

Figure 3: Linear stability of bumps in a feedforward two-layer network, demonstrated by simulations of the model Eq. (1) with . (A) When the bump in layer 1 is shifted, the bump in layer 2 (dashed line) relaxes to the new position of bump 1 (solid line). (B) When the bump in layer 2 is shifted, it relaxes back to the fixed position of bump 1. (C) When both bumps are shifted, both retain their new position, respecting the translation symmetry of the underlying Eq. (1).

Two-layer feedforward network. Assuming , layer 1 sends input to layer 2, but receives no feedback from layer two. Linear stability associated with the stationary bump solutions to this model is then determined in part by computing the eigenvalues of:

(26)

where . Since the matrices in Eq. (26) are triangular, their eigenvalues are given by their diagonal entries. Applying Eq. (24), , we can express eigenvalues of as:

Neutral stability with respect to the eigenfunction ensures the existence of . Concerning the other three eigenvalue formulae, the terms will be positive by our assumptions on the weights made after Eq. (24), so , corresponding to the fact that the bump in layer 2 is linearly stable to translating perturbations, when the position of the bump in layer 1 is held fixed. In Fig. 3, we show that the upstream layer (1) governs the long term location of both bumps. The layer 2 bump always relaxes to the layer 1 bump’s location. In a related way, the eigenvalues and correspond to expansions/contractions of the bump widths in layers 1 and 2, respectively. Typically, there are two bump solutions in a single-layer network, whose width perturbations correspond with the eigenvalue : one that is narrow and unstable to such perturbations (), and another that is wide and stable to such perturbations ([2, 22, 41]. Lastly, the bump in layer 2, driven by activity in layer 1 is influenced by features of layers 1 and 2, as shown in the formula for . When , we expect , and the bump will be stable to width perturbations.

Exploding star network. Another example architecture involves a single layer with feedforward projections to multiple () layers. In this case, for and . Only perturbations that shift the bump in layer 1 have a long term impact on the position of bumps in the network. The translation modes of bumps in layers have associated negative eigenvalues, as we will show, which is a generalization of the two-layer feedforward case. Linear stability is computed by first determining the eigenvalues of:

(27)

Subtracting one from the eigenvalues of the matrices in Eq. (27) and applying the formula for , Eq. (24), we find eigenvalues, given by for , where

As in the two-layer network, bumps are neutrally stable to perturbations of the form , corresponding to . In addition, we expect for since . As mentioned above, we would expect wide bumps to be stable to expansion/contraction perturbations, whose stability is described by the eigenvalues for .

Figure 4: Linear stability of bumps in a two-layer symmetric recurrent network with , demonstrated by simulations of the model Eq. (1) with . Bumps initially at the same position are perturbed to examine the resulting evolution of their positions. (A) When the bump in layer 1 (solid line) is shifted, both bumps relax to the average of their initially perturbed position. (B) When the bump in layer 2 (dashed line) is shifted, again, both bumps relax to an intermediate position. (C) When both bumps are shifted to the same location, both retain their new position.

Two-layer recurrent network. In the case of a fully recurrent network, where for all , all matrix entries are nonzero: . First, note that is an eigenvalue of , since , so we denote as the eigenvalue describing the translation symmetry of bumps. To gain further insight, we can also compute the other three eigenvalues:

For a symmetric recurrent network: , , and (, ), these formulas reduce to and . Bumps are linearly stable to perturbations that move them apart (Fig. 4A,B), and neutrally stable to translations that move them to the same location (Fig. 4C).

Figure 5: Linear stability of bumps in a layer imploding star network with for , . (A) When the bump in layer 1 (solid line) is shifted, the bump in layer 3 (dashed line) relaxes to a position between the layer 1 and layer 2 (dotted line) bump. (B) When both bumps in layers 1 and 2 are perturbed to a new location, the bump in layer 3 relaxes to that new location. (C) When only the bump in layer 3 is perturbed, it relaxes back to the locations of the bumps in layer 1 and 2.

Imploding star graphs. Additional dimensions of neutral stability can arise in the case of more than two layers, depending on the graph defining interlaminar connectivity. For instance, if there are multiple layers that receive no feedback from other layers, then for , so the first rows of the matrix are the canonical unit vectors . Thus, there are at least unity eigenvalues of , implying has multiplicity at least in Eq. (25), corresponding to the neutral stability of bumps in the layers that receive no feedback. We consider such an example when :

Eigenvalues of Eq. (25) are then and for , and

Bumps in both layers 1 and 2 are neutrally stable to translations (Fig. 5A,B), whereas the bump in layer 3 is linearly stable to translation, relaxing to a weighted average of the positions of the layers 1 and 2 bumps (Fig. 5C). Adding dimensions to the space of translation symmetric perturbations will change the low-dimensional approximation that captures the dynamics of multilayer bump solutions in response to noise perturbations (Compare Sections 3.3 and 5).

Directed loop of layers. As a last example, we consider an -layer directed loop, wherein each layer provides feedforward synaptic input to a subsequent layer. As a result, there is a band of nonzero interlaminar coupling along for (replace with ). Again, there is a zero eigenvalue in Eq. (25), since

Our desired result can be demonstrated by computing the determinant of the bidiagonal matrix: