subsecref subsecname = thmref thmname = theorem lemref lemname = lemma

# Efficient low-dimensional approximation of continuous attractor networks

Alexander Seeholzer${}^{\text{1}}$, Moritz Deger${}^{\text{1,2}}$, Wulfram
Gerstner${}^{\text{1}}$

###### ?abstractname?

Continuous “bump” attractors are an established model of cortical working
memory for continuous variables and can be implemented using various neuron
and network models. Here, we develop a generalizable approach for the
approximation of bump states of continuous attractor networks implemented in
networks of both rate-based and spiking neurons. The method relies on a low-
dimensional parametrization of the spatial shape of firing rates, allowing to
apply efficient numerical optimization methods. Using our theory, we can
establish a mapping between network structure and attractor properties that
allows the prediction of the effects of network parameters on the steady state
firing rate profile and the existence of bumps, and vice-versa, to fine-tune a
network to produce bumps of a given shape.

1: School of Computer and Communication Sciences and School of Life Sciences,
Brain Mind Institute,
École Polytechnique Fédérale de Lausanne, 1015 Lausanne EPFL, Switzerland
2: Institute for Zoology, Faculty of Mathematics and Natural Sciences,
University of Cologne, 50674 Cologne, Germany

_A previous version of this article was published as Chapter 3 of the first
author’s Ph.D. thesis[1]_

###  1 Introduction

Behaving animals commonly need to transiently memorize information about the
environment. For example, as an animal overlooks the visual scenery, locations
of certain salient stimuli need to be recorded and stored. Such information
does not need to be stored in long-term memories. Rather, working memory must
provide a quickly accessible computational substrate for storing information
over short durations. While long-term memory is thought to be stored in the
efficacy of synaptic connections in the brain [2, 3, 4], a possible substrate
for working memory may be transiently stable states of neuronal activity
across cortical networks [5, 6, 7].

As model implementations of this concept, localized spatial profiles of neural
activity have been proposed for the internal representation of sensory stimuli
[8, 9, 10]. First, neurons are associated to the presence of physical
quantities through elevated responses during and after the presentation of
stimuli, akin to receptive fields. For example, the presentation of stimuli at
varying angular positions in the visual field evokes persistent and elevated
firing rates in selective groups of neurons of the prefrontal cortex during
delay periods [11]. Choosing recurrent connection weights (or connection
probabilities) which are stronger between neurons that are responsive to
similar stimuli, together with feedback inhibition limiting the total firing
rates in a network, allows this class of models to display _bumps_ of self-
sustained activity: neuronal activity that is localized in the space of
possible stimuli. Since these states are stable attractive states, and all
possible such states form a continuum, these models are often referred to as
_continuous attractors_. The elevated firing of neurons responsive to similar
stimuli is then seen as the working memory representation of physical quanta
stored in the network, e.g. spatial orientations [12], or angular positions in
the visual field [13]. Similar computational circuits might also serve as the
basis of persistent internal representations in hippocampal areas [14, 15].

Continuous attractor models with simplified shapes of connectivity or neuronal
input-output relations can be analyzed and often exactly solved [8, 9, 10, 16,
17, 18, 19], or may generally be approximated in the linear input-output
regime of balanced networks [20]. However, the inclusion of biologically
plausible nonlinearities, like nonlinear neuronal input-output relations [13,
21], neuronal adaptation [22, 23], or synaptic nonlinearities like short-term
plasticity [18, 24] and saturating NMDA kinetics [13, 25], complicate the
mathematical solution of these systems considerably and make a derivation of
the stable firing rate profile unfeasible. Therefore, such systems are usually
studied by explicit simulations of the underlying dynamics or by numerical
optimization of approximated equations for all neurons [20, 26]. While these
procedures in principle allow the prediction of the network activity as a
function of the parameters, they involve computationally demanding numerical
optimization of high-dimensional systems of equations, possibly as costly as
simulating the full neuronal dynamics. Thus, currently, relating the
microscopic network parameters to the resulting emergent bump states involves
repeated and possibly time consuming simulations. For example, this makes the
matching of the network steady states to physiologically constrained features
tedious.

Here, we present a generalizable approach for the approximation of the
network-wide steady states of continuous attractors. Our approach allows the
prediction of the shape of steady-state firing rate profiles, under nonlinear
neuronal dynamics and varying configurations of the underlying microscopic
system, without having to solve the dynamics of the full, high-dimensional
system. Our novel method relies on a low-dimensional parametrization of the
network’s firing rate profile, which allows us to derive computationally
tractable systems of self-consistent equations describing the attractor
steady-state, akin to mean-field approaches for networks with discrete
attractors [27]. These equations can be used to efficiently predict the
dependence of the firing rate profile on microscopic network parameters.
Importantly, because the dimension of the parameterization of the spatial
activity profile is low, our approach makes optimization of the microscopic
network parameters for the appearance of desired bump profiles feasible. We
apply our method to both networks of simplified rate neurons, and networks of
complex, conductance-based integrate-and-fire neurons with saturating and
voltage-dependent nonlinear NMDA transmission.

###  2 Results

Mean-field approaches (see e.g. [27, 28]) that predict the steady states of
recurrently connected neuronal networks usually rely on dimensionality
reduction. The number of equations describing the dynamics is reduced by
partitioning neurons into groups of “similar” neurons, and deriving
expressions which describe the average statistics for these coupled groups in
the steady states. For example, the simplest such partition consists in
considering the mean firing rates of excitatory and inhibitory neurons
separately, e.g. all excitatory neurons fire with similar mean rates given the
same input:

$\nu_{E}=F_{E}\left(\text{input}_{E}\right).$

If the groups of neurons are now homogeneously coupled, i.e. the connections
between neurons depend only on the groups of the neurons involved, one can
derive the input to neurons of each group in dependence of the firing rates of
the groups only. This leads to a closed system of self-consistency equations
describing the coupled steady-state firing rates:

$\displaystyle\nu_{E}$
$\displaystyle=F_{E}\left(\text{input}_{E}\left[\nu_{E},\nu_{I}\right]\right),$
$\displaystyle\nu_{I}$
$\displaystyle=F_{I}\left(\text{input}_{I}\left[\nu_{E},\nu_{I}\right]\right).$

In the steady states of continuous attractor models (see 1), neurons fire at
different rates, making a clear partition into discrete groups of similarly
firing neurons difficult. Therefore, the solution of such systems usually
relies on the explicit simulation of the neuronal dynamics of all neurons
along the spatial dimension, or a numerical solution of the coupled self-
consistency equations for all neurons.

Here, by using the continuity of the shape of the attractor states, we
demonstrate that continuous attractors are amenable to dimensional reduction,
by parametrizing the attractor state by a low-dimensional family of functions.
In 2.2, we check our method on networks of simple rate neurons, for which the
method might not yield much improvement over simulations or numerical
solutions of the steady states. For the spiking networks considered in 2.3, we
show that our approach speeds up predictions of the steady states considerably
and further makes these networks amenable to the optimization of network
parameters.

####  2.1 General equations for the approximation of stable states in ring-
attractors

?figurename? 1: Stable firing rate profiles in continuous attractors. A:
Neurons (circles) are assigned a position on a ring and connected with
distance-dependent connection weights to all other neurons: the firing of a
neuron (red circle) inhibits all other neurons (global inhibition, blue lines)
and strongly excites neurons close it (local excitation, red lines). B:
Example simulations of a continuous attractor. At times $t=0,1,2s$ the system
is reset to localized activity centered at different positions. With time the
activity bump broadens towards a stable state. Colors indicate firing rates as
in panel C. C: Firing rates of the network shown in panel B, measured close to
the stable states at $t=0.9,1.9,2.9s$. Plots are generated using the rate
model introduced in 4.1.

As a concrete class of continuous attractors, we consider the ring-attractor
model, in which stable bumps of neuronal activity are freely translatable
along all positions on a circle. Ring-attractor models can be constructed by
placing $N$ neurons (rate-based or spiking) at equally spaced angular
positions $\theta$ along the ring (1A) [8, 9, 10]. We choose the angular space
to consist of positions $\theta\in[-\pi,\pi)$, where we identify the ends of
the interval: a neuron at position $\theta=\pi-\epsilon$ is the neighbor of a
neuron at position $\theta=-\pi.$ At short angular distances, recurrent
connections are chosen to be strong and excitatory, while neurons further
apart in angular space effectively inhibit each other’s firing (1A). Due to
the symmetry in connectivity with respect to distance, these networks can form
a continuous manifold of stable states for sufficiently strong connections:
the network activity in response to external inputs converges to firing rate
profiles centered around some angular position. The position can be, for
example, controlled by providing an external input to the network centered
around any desired position (1B). The stereotypical shape of the resulting
firing rate profiles (1C) is invariant with respect to translations in angular
space.

The continuity of the firing rate profile allows us to parametrize the firing
rates in the population by a small number of parameters, and to derive
equations from the underlying model that constrain these parameters. Here,
inspired by shapes observed in simulations, we choose to parameterize the
firing rate profile by a generalized Gaussian function, where we assume
without loss of generality that the distribution is centered at $\theta=0$
(cf. 2A):

$g(\theta)=g_{0}+g_{1}\exp\left(-\left[\frac{\left|\theta\right|}{g_{\sigma}}\right]^{g_{r}}\right).$
(1)

?figurename? 2: Approximation of bump shapes in attractor networks. A:
Parametrization of firing rate profiles by 4 shape parameters: $g_{0}$ and
$g_{1}$ control the baseline and maximal firing rates, $\sigma_{r}$ controls
the width, and $r$ controls the steepness. B: Assuming a spatial profile given
by the parametrization $g$, the input to any neuron can be calculated by
summing all synaptic inputs: in the example, the neuron at position $0$
receives some strong local excitatory input (red thick arrows) and weaker
inhibitory input from more distant neurons (blue thin arrows). The firing rate
prediction (dashed line) for any neuron can then be calculated as a function
of its input. This is illustrated for a neuron at $0$ (circle). C: Self-
consistency errors between the current parametrization (gray line) and the
resulting firing rate prediction (dashed line) are calculated at a small
number of points along the firing rate profile. Intermediate points
(triangles) are positioned dynamically during optimization. D: Optimized bump
parametrization (dashed lines) for systems with three different connectivities
(solutions of the full system are plotted in color). E: Dependence of solution
on intermediate point placement. Left points are given by downward triangles,
rightward points by upward triangles (compare to panel C). Single optimization
runs (light gray lines) together with median parameters (dashed lines, same as
in D) and full solutions (colored lines).

Here, $g_{0}$ controls the baseline firing rate and $g_{0}+g_{1}$ will be the
maximum firing rate of the profile. The parameters $g_{\sigma}$ and $g_{r}$
control the width and steepness of the profile, respectively. If we know the
distance-dependent connectivity $w$ between neurons and their input-output
relation $F$, we can predict the expected neuronal firing at any position
$\theta$ in the population (cf. 2B). Crucially, we want the firing rate
profile $g(\theta)$ to be generated by the neuronal dynamics – we thus
identify $g(\theta)$ at the point $\theta$ with the firing rate $\nu(\theta)$
of a neuron at position $\theta$. Finally, we replace the synaptic input to
the neuron at position $\theta$ with the contributions from all neurons firing
at rates $g$ along the ring. For any given position $\theta_{i}$ along the
ring, this yields a self-consistent equation in the function $g$:

$\displaystyle g(\theta_{i})$
$\displaystyle=F(\text{input}_{\theta_{i}}\left[g\right])$
$\displaystyle=F\left(\int_{-\pi}^{\pi}d\varphi
w(\varphi-\theta_{i})g(\varphi)\right),$ (2)

with the corresponding _self-consistency error_

$\displaystyle\text{Err}_{i}$ $\displaystyle\equiv
g(\theta_{i})-F(\text{input}_{\theta_{i}}\left[g\right]).$ (3)

In principle, this procedure can yield up to $N$ coupled error functions, one
for each of the $N$ neurons. One could then minimize the quadratic error
$\sum_{i}\text{Err}_{i}^{2}$ with respect to the parameters
$\\{g_{0},g_{1},g_{\sigma},g_{r}\\}$ to find an approximate solution of the
system. However, since the evaluation of each error function can be costly
(e.g. in spiking networks, see 2.3.2), we propose a low-dimensional
approximation to constrain the set of $4$ free parameters: we pick only $4$
points $\theta_{1},\dots,\theta_{4}$, at which we evaluate the errors. This
assumes that the shapes of firing rate profiles maintained by the network are
well approximated by the function $g(\theta)$, which we found to be the case
for all networks considered here. This leaves the choice of points
$\theta_{i}$ to evaluate. To ensure that errors are evaluated across different
firing rates, we set the position of these points to cover a range of function
values $h_{i}=g(\theta_{i})$: we choose the top of the distribution
$\theta_{1}=0$ with $h_{1}=g_{1}+g_{0}$, as well as the lowest point
$\theta_{4}=\pi$ with $h_{4}=g_{0}$ (circles in 2C). The remaining
intermediate points (triangles in 2C) are dynamically positioned (see 4.3.1 in
Methods): their position depends on the function $g$ such that they always
sample intermediate function values $h_{i}$.

####  2.2 Approximation of ring-attractor profiles in rate models

The proposed method is, in principle, applicable to any neuron model with a
defined input-output relation $F$. The shape of the stable attractor profiles
will, however, depend on the concrete choice of neuron model and the
microscopic parameters, in particular the parameters governing the
connectivity between neurons. To test the ability of the low-dimensional
approximation proposed in the last section to correctly predict the shapes of
firing rate profiles, we implemented the ring-attractor model introduced above
(see 1) in a network of rate-based neurons with $\tanh$ input-output function
and a generalized Gaussian recurrent connectivity (see 4.1 in Methods). For
the rate neuron models chosen here, the self-consistency errors (LABEL:error-
func) are given by (see 4.3.2 in Methods):

$\displaystyle\text{Err}_{i}$
$\displaystyle=g(\theta_{i})-\frac{\nu_{\text{max}}}{2}\left[1+\tanh\left(\frac{\tau_{s}}{s_{0}}\frac{N}{2\pi}\text{input}_{\theta_{i}}\left[g\right]\right)\right],$
$\displaystyle=g(\theta_{i})-\frac{\nu_{\text{max}}}{2}\left[1+\tanh\left(\frac{\tau_{s}}{s_{0}}\frac{N}{2\pi}\int_{-\pi}^{\pi}d\varphi
w(\varphi-\theta_{i})g(\varphi)\right)\right],$ (4)

where $\tau_{s}$ is the time-constant of synaptic inputs, $\nu_{\text{max}}$
is the maximal firing rate and $s_{0}$ is an input scale.

#####  2.2.1 Prediction of stable firing rate profiles

To approximate the firing rate profiles that these networks admit as self-
consistent solutions, we minimize the error functions (LABEL:rate-error) with
respect to the parameters $\\{g_{0},g_{1},g_{\sigma},g_{r}\\}$ as free
variables. We find that, for a range of connectivity parameters (see 2), the
predicted shapes converge to unique solutions. This solution matches the
steady state of the microscopic network simulations accurately (2D). This is
the case both for attractor states that lie in the linear regime of the
neuronal input-output relations (2D, red line) as well as for highly nonlinear
attractor dynamics in which neuronal firing reaches saturation values, leading
to plateau-shaped firing rate profiles (2D, blue and green lines).

As discussed above, the placement of intermediate sampling points (2C,
triangles) is not constrained by theory, but remains a free parameter of our
approach. We chose these points at positions $\theta_{i}$ such that they
sample given function values $h_{i}=g(\theta_{i})$. To investigate the
dependence of the prediction on the placement of intermediate sampling points,
we calculated several predictions while randomly varying the choices of
$h_{i}$ (2E, triangles). We find that this hardly affects the converged
solutions.

#####  2.2.2 Optimization of network parameters

?figurename? 3: Application: Optimization of network connectivity. A: Stable
network firing rate profiles resulting from 10 optimization runs of the
network local connectivity profiles. The optimization target profiles are
plotted in colors, together with single optimization runs (light gray lines,
almost overlapping with colored lines) and one example highlighted
optimization result (dashed lines) each. Left optimization points are given by
downward triangles, rightward points by upward triangles (compare to 2C,E). B:
Results of the 10 optimized connectivity profiles (light gray lines) together
with one example profile (colored lines). Solutions for System $2$ are
degenerate, while others are fairly unique.

While our low-dimensional system of self-consistent equations can be used for
the prediction of the steady-state firing rates, they can also be used in the
inverse way, to optimize any of the network parameters. We demonstrate this
here, using the shape of the recurrent connectivity as an example. However,
such optimizations can include further parameters of the network model (see
2.3.3).

To optimize the network connectivity parameters, we keep the parametrization
parameters $\left\\{g_{0},g_{1},g_{\sigma},g_{r}\right\\}$ fixed to the
desired values of the shape of the firing rate profile. We then optimize the
self-consistent equations (LABEL:mfeq-selfconsist) for values of the recurrent
connectivity parameters $\left\\{w_{0},w_{1},w_{\sigma},w_{r}\right\\}$ which
lead to solutions of the equations and produce the desired bump profile. In 3
we show the results of this procedure for the three systems also investigated
in 2. The procedure yields network connectivities that fulfill the desired
properties (3A), largely independently of the points $\theta_{i}$ chosen for
the evaluation of the errors. Importantly, for some shapes the solutions show
a degeneracy (3B, gray lines), in the sense that several connectivity
parameter sets are found that produce the same stable firing rate profile.

####  2.3 Approximation of ring-attractor profiles in spiking networks

?figurename? 4: Application: Prediction of firing rate profiles in spiking
continuous attractor networks. A: Recurrently connected spiking network of a
population of excitatory (E, red triangles) and inhibitory (I, blue circles)
neurons. Networks are fully connected with uniform weights, except for E-E
connections (dashed red line), which are distance-dependent. All neurons
receive excitatory input with spikes generated by homogeneous Poisson
processes. B: Example simulation: E neurons fire asynchronously and
irregularly until an external cue is given (centered at $0$ for $0.2s$
starting at $t=1.8s$). After this stimulus, a bump of elevated activity
sustains itself around the point of stimulation ($w_{+}=2.0$). C: Distance-
dependent E-E connectivity as a function of the parameter $w_{+}$ (maximal
strength or recurrent connections). Values range between $w_{+}=1$ (red) and
$w_{+}=3$ (blue). D: Appearance of the bump firing rate profile as a function
of the connectivity parameter $w_{+}$. Theoretical predictions (lines) and
simulation results (maximum of fit of $g$ to firing rates recorded over $1.5$s
of delay activity, mean over $5$ repetitions, errors show $95\%$ CI) for
maximal firing rates of E (colored) and I (gray) neurons. Colors similar to
those in panel C. E: Shape of the optimized firing rates profiles (solid
colored lines) compared to mean firing rates measured from spiking simulations
(thin colored lines) (mean $\pm 0.5\text{STD}$ of firing rates recorded from
$1.5s$ of delay activity in $5$ repetitions) for all values of $w_{+}$ in
panel C and D (similar colors).

In complex spiking neuron models, the steady-state input-output relations
often involve integral functions [27, 29] that are not amenable to further
theoretical analysis. Further, the introduction of nonlinear NMDA transmission
at excitatory synapses [13, 30] complicates the analysis of such models
considerably: voltage-dependent gating of the maximal NMDA conductance (by the
voltage-dependence of the $\text{Mg}^{2+}$ block) and saturation of the NMDA
at high conductances necessitates the numerical solution of a $4$-dimensional
system of self-consistent equations for a relatively simple 2 population mean-
field model [27] (see below).

Here, we demonstrate that firing rate profiles in a continuous attractor
network of spiking neurons with such nonlinear NMDA transmission are amenable
to the same approach as described above, which still involves evaluating
comparatively few equations. The spiking network we implement (similar to [13]
and used with variations in e.g. [21, 31, 32, 33, 34], see 4.2 in Methods)
consists of two fully connected populations of conductance-based integrate-
and-fire neurons: a population of inhibitory neurons with unstructured all-to-
all connectivity, and a population of excitatory neurons, with distance-
dependent recurrent excitatory connections (4A). In addition, all neurons
receive excitatory background input mediated by spikes generated by Poisson
processes. These network can be tuned such that they possess a bi-stability
(4B): a uniform state with spontaneous spiking activity in the excitatory
population (the inhibitory population is always uniformly spiking) coexists
with an “evoked” spatially inhomogeneous bump-state that appears after an
external cue input is given to a subgroup of excitatory neurons (stimulus is
present at $t=1.8-2.0s$ in 4B).

#####  2.3.1 Self-consistent equations for networks of spiking neurons

For the excitatory population, we again parametrize the spatial profile of
firing rates by (LABEL:g-gen-gaussian), which allows us to derive self-
consistent equations for any neuron in the excitatory population. We construct
self-consistent equations for the excitatory firing rates at positions
$\theta_{i}$ as in (LABEL:mfeq-selfconsist). However, these now will depend
additionally on the inhibitory firing rate $\nu_{I}$. Also, the voltage-
dependence of the differential equations leads to an additional self-
consistent equation for the mean-voltage $\bar{V}$. For any position $\theta$,
the excitatory self-consistent equations are of the form (see 4.3.3 in Methods
for detailed expressions):

$\displaystyle g(\theta)$
$\displaystyle=F(\text{input}_{\theta}\left[g\right],\nu_{I},\bar{V}(\theta))$
$\displaystyle\equiv F\left(\int_{-\pi}^{\pi}d\varphi
w(\varphi-\theta)\psi\left(g(\varphi)\right),\nu_{I},\bar{V}(\theta)\right),$
(5) $\displaystyle\bar{V}(\theta)$
$\displaystyle=G\left(\text{input}_{\theta}\left[g\right],\nu_{I},\bar{V}(\theta)\right).$
(6)

The function $\psi(g)$ expresses the mean synaptic activation under
presynaptic Poisson spiking at rate $g$. For accuracy, we chose to measure
$\psi$ numerically for the model of nonlinear NMDA conductance of the
recurrent excitatory synapses given in the network (see 4.3.3 in Methods).

To constrain the free parameters of $g(\theta)$, we again pick $4$ points
$\theta_{i}\in\left\\{\theta_{1},\dots,\theta_{4}\right\\}$, each now yielding
2-dimensional error functions

$\displaystyle Err_{i}$
$\displaystyle=\left(\begin{array}[]{c}g(\theta_{i})-F\left(\text{input}_{\theta_{i}}\left[g\right],\nu_{I},\bar{V}(\theta_{i})\right)\\\
\bar{V}(\theta_{i})-G\left(\text{input}_{\theta_{i}}\left[g\right],\nu_{I},\bar{V}(\theta_{i})\right)\end{array}\right).$
(9)

The resulting 8 equations are optimized for the $4$ parameters of the
parametrization $g$, as well as the additional $4$ variables
$\bar{V}(\theta_{i})$. The inhibitory population, on the other hand, is
assumed to be homogeneous. Its activity can be described by a single mean
firing rate $\nu_{I}$ and the average voltage in the inhibitory population
$\bar{V}_{I}$, resulting in a pair of additional self-consistency errors that
constrain these two variables:

$\displaystyle Err_{I}=$
$\displaystyle\left(\begin{array}[]{c}\nu_{I}-F(\text{input}_{I},\nu_{I},\bar{V}_{I})\\\
\bar{V}_{I}-G(\text{input}_{I},\nu_{I},\bar{V}_{I})\end{array}\right),$ (12)

where
$\text{input}_{I}=\frac{1}{2\pi}\int_{-\pi}^{\pi}d\varphi\psi\left(g(\varphi)\right)$
is the mean recurrent excitatory input to inhibitory neurons.

As mentioned above, in a certain range of parameters the spiking system
possesses two dynamically stable states (4B): the uniform state and the
“evoked” bump state. In this bistable regime, the associated self-consistency
Equations (5)-(6) must have an an additional unstable solution [35]. Even for
parameters in which the bump-state is the only stable state of the system, the
uniform state will still be a (unstable) solution of the self-consistent
equations. Accordingly, numerical solutions of the errors Eqs. (9)-(12)
sometimes converge to the uniform state or an unstable intermediate solution,
even if a stable bump state at higher firing rates exists. In the following we
consider only the solutions with the highest spatial modulation found under
repeated solutions (see Discussion).

#####  2.3.2 Prediction of firing rate profiles from network properties

?figurename? 5: Efficient optimization of self-consistency errors. A:
Evaluations of self-consistency errors until convergence for three different
optimization methods ($w_{+}=2.5$) ($10$ repetitions, errors show $95\%$ CI).
_L-BFGS-B_ and _SLSQP_ both minimize the sum of squared errors, _Root_
minimizes the error vector directly. Right-hand axis shows wall-clock times of
optimization procedures. B: Averaged root mean square error (RMSE) between
optimized firing rate profiles and mean firing rates measured from simulations
($w_{+}=2.5$, data from 4C,D,E). RMSE was calculated for all ($800$) neurons
($n=10$, as in panel A). C: Optimized firing rate profiles, together with
placement of sampling points for two combinations of optimizers and points
(Root on $4$ points, SLSQP on $12$ points). Gray line shows mean firing rate
measured from simulations ($w_{+}=2.5$, as in 4C,D,E).

Above, we derived error functions constraining the parametrization of firing
rate profiles for spiking networks (LABEL:spike-selfcons-err)-(12). Here, we
use these to predict the dependence of the spatial shape of the firing rate
profile on of a bifurcation parameter $w_{+}$, which is the maximal strength
of recurrent excitatory connections. At $w_{+}=1$ the connection profile is
homogeneous, while at larger values local connections are stronger (4C). The
strength of long range connections is calculated by a normalization condition
(see 4.2.3 in Methods).

As $w_{+}$ is increased (4C) above a critical value, a spatially inhomogeneous
bump state appears in simulations of the spiking network (4E). Our theory
predicts this dependence of the network state on the bifurcation parameter,
while approximating to a large degree the changing shape of the rate profile
as the parameter is increased (4D,E). The firing rates of the inhibitory
population and their increase with the parameter $w_{+}$ are also well
described (4D, black dots and lines).

As mentioned above, the error functions (LABEL:spike-selfcons-err) could, in
principle, be evaluated at an arbitrary number of points. To constrain the $4$
parameters of the parametrization $g(\theta)$, we chose only the minimal
number of points. This reduces the necessary number of evaluations of the
errors $\text{Err}_{i}$ (5A). Further, since in this case the dimensions of
the optimization variables $\left\\{g_{0},g_{1},g_{\sigma},g_{r}\right\\}$ and
the error vector coincide, application of a more efficient numerical
optimization method (_Root_ , see 4.3.4) allows for faster optimization (5A,
green bar), which reduces the needed time from $\sim 30s$ to close to $2s$
(5A, right hand axis). We observe, however, that adding additional points does
slightly influence the resulting prediction (5B), where more sampling points
placed in the flanks of the bump tend to reduce slightly the predicted maximal
firing rates (5C, orange points).

In a second experiment, we show that the theory can be used to efficiently
predict the effect that changing network parameters have on the shape of the
resulting firing rate profile. Similar to the simulated experiment in [31,
Fig. 3], we systematically reduced the strength of recurrent excitatory-to-
excitatory ($g_{\text{EE}}$) and inhibitory-to-excitatory ($g_{\text{EI}}$)
conductances of the network from the baseline of the network presented in 4.
Such changes of the ratio of excitation to inhibition have been hypothesized
to occur under cortical disinhibition observed in schizophrenia [31, 36].
Recovering the result presented in the study, we see that the width of the
bump profile111Note that the network presented here generally has a wider
profile than the one investigated in [31]. depends mostly on the ratio of
recurrent conductances, and thus undergoes significant widening under
disinhibition. As we have shown in 5A, the optimization procedure for each
datapoint is comparatively fast and thus enables these type of parameter scans
for wide ranges of values under many parameters.

These results show that our approach can be used to accurately describe the
firing rate profiles of bump-attractor networks of recurrently connected
spiking excitatory and inhibitory neurons, across a range of parameters. While
evaluating the error function at more points can lead to slightly increased
accuracy of the prediction, the impact on optimization performance is
significant, increasing the needed time by an order of magnitude.

#####  2.3.3 Optimization of network parameters

?figurename? 6: Applications: Shape prediction over parameter ranges /
Optimization of network parameters. All panels use the spiking network model
of 4. A Theoretically predicted width of the firing rate profiles (full width,
as predicted by the parameter $2\cdot g_{\sigma}$) under varying strengths of
recurrent excitatory-to-excitatory ($g_{\text{EE}}$) and inhibitory-to-
excitatory ($g_{\text{EI}}$) synaptic connections. Compare to the similar
experiment performed in spiking network simulations in Murray2012 [31, Fig.
3]. The bifurcation parameter was kept fixed at $w_{+}=2.5$. B Network
parameters optimized for 3 different shapes of the stable firing rate profile
(see 3 for details). Solid dark lines are theoretical predictions using the
same method as in 2.3.2. Lighter lines and shaded areas are mean $\pm
0.5\text{STD}$ of firing rates recorded over $3.5s$ of delay activity in 5
repetitions.

As demonstrated above for a rate-based network, our low-dimensional
approximation of continuous attractors allows the optimization of network
parameters. Here, we demonstrate that this approach extends well to continuous
attractors implemented in recurrently connected spiking neural networks with
nonlinear synaptic transmission. As in the case of the rate network, this is
achieved by fixing some desired properties of the firing rate profile, while
minimizing the error functions Eqs. (9)-(12) with respect to several network
parameters. Here, we included the shape parameters $w_{\sigma}$ and $w_{+}$ of
the distance-dependent connectivity, as well as the strengths of all recurrent
synaptic connections:
$g_{\text{EE}},\,g_{\text{IE}},\,g_{\text{EI}},\,g_{\text{II}}$ (cf. 3 for
details).

To find networks that admit a given shape of the firing rate profile, we first
fixed the firing rate modulation $g_{1}$ and width $g_{\sigma}$, while
optimizing the connectivity parameters $w_{+},w_{\sigma}$ as well as the
strength of all recurrent excitatory and inhibitory transmission. In total,
there were $17$ variables (see 3 for listings and optimization results) which
were optimized (see 4.3.4 for details). Varying these free parameters allows
us to optimize the remaining parameters of the spiking networks through a
range of shapes of the stable firing rate profile, from rather thin bumps of
high activity (6, red) to wide bumps of low activity (6, blue). To check
whether we could optimize the spiking network to show saturated flat-top
shapes at low firing rates, we fixed $g_{1}=25Hz$ and the sharpness parameter
to a high value $g_{r}=8$, while optimizing the parameter $g_{\sigma}$ (6,
green).

Similar optimization results for all three bump shapes were achieved by
imposing additionally a low basal excitatory firing rate
$\nu_{E}^{\text{basal}}$, which constrains firing rates in the uniform
(spatially homogeneous) state for $w_{+}=1$. For the optimizations above, the
basal rates were unconstrained and varied between $1Hz$ and $5Hz$. Thereby,
the number of optimization variables is increased to $20$, since additionally
a basal inhibitory firing rate and basal inhibitory and excitatory mean
voltages need to be introduced to calculate the self-consistency error of the
basal firing rate $\nu_{E}^{\text{basal}}$. It should be noted, that the value
of this constraint affects the possible bump shapes: for example, setting
$\nu_{E}^{\text{basal}}=1Hz$ did not yield a converging optimization for the
blue and green curves of 6 – this could be alleviated by relaxing to the
higher value $\nu_{E}^{\text{basal}}=3Hz$.

###  3 Discussion

We have presented a framework for the approximation of self-sustained firing-
rate profiles in continuous attractor neural networks. Analytical computation
of the steady states of continuous attractors is often not possible, since it
involves the solution of high-dimensional systems of nonlinear equations,
making numerical solutions necessary. Moreover, the spatially inhomogeneous
firing rate profiles of these networks prohibit dimensional reduction of these
equations by separation of neurons into homogeneous populations, as is usually
done in mean-field approaches [27, 28]. Here, we propose a simple approach,
consisting in approximating the continuous firing rate profiles by a family of
functions $g$ with only $4$ parameters. These parameters are constrained by
equations expressing the microscopic dynamics of the neurons and synapses
involved in the model, and can be optimized to find admissible solutions. As
we have shown, this can be used for the efficient mapping of the effects that
different network parameters have on the bump shape. Next to predicting the
emergent steady states of attractors, the utility of the low-dimensional
approximation is that the derived self-consistent equations are efficiently
optimizable: we were able to use standard numerical methods to constrain the
parameters of spiking networks to show desired firing-rate profiles.

In the main text we have formulated our approach as generally as possible, to
emphasize that the method does not rely on a specific neuronal (rate or
spiking) model, as long as a prediction of firing rates given the synaptic
input can be derived. Therefore, the theory could be extended easily to other
neuron models, or connectivities where connection probabilities are distance-
dependent (e.g. [37]). The approach presented here could also be used to
predict the activity of two-dimensional attractor models implemented in
“sheets” of neurons, that are often used in the context of hippocampal
networks [14, 38, 39]. Assuming isotropy of the connectivity (if connection
strengths depend only on the Euclidean distance between neurons), a two-
dimensional generalized Gaussian function with $g_{0},g_{1},g_{r}$ and a
single parameter $g_{\sigma}$ as before could be used to approximate activity
states. For non-isotropic networks, further width parameters could be
introduced and constrained by sampling at additional points.

We have shown that our approach is amenable to the inclusion of synaptic
nonlinearities like saturating NMDA transmission (which is captured by the
synaptic activation function $\psi$, cf. (LABEL:spike-selfcons-rate)). Other
sources of nonlinear synaptic transmission, for example the activity-dependent
short-term plasticity of synapses [24] that is often investigated in the
context of working memory models [18, 37, 40], can be similarly incorporated:
calculation (or numerical estimation) of a compound function $\psi$ that
describes the steady-state values of synaptic activation under all nonlinear
processes affecting synaptic transmission would suffice to adapt the
theory222See Chapter 4 of [1], where we apply the same method to networks with
short-term synaptic plasticity.. Similar to the estimation of the mean-voltage
in spiking networks demonstrated here, the theory could also be extended to
incorporate adaptation effects on the steady-state firing rates of neurons
[22].

For the prediction of firing rates in spiking networks we have adapted a
theory for the description of mean-field firing rates of conductance based
integrate-and-fire neurons [27] to predict spatial firing rate profiles in
recurrently connected networks [13]. Mean firing rates and mean voltages in
this theory are generally expectation values over ensembles of neurons that
can be assumed to have homogeneous activity. Strictly speaking, here we
violate this assumption by taking the rate prediction of theoretical ensembles
as the prediction of the firing rates (and mean voltages) of single neurons at
given positions along the ring attractor. However, since we are investigating
the stationary state of networks, the mean firing rates calculated from this
theory can also be interpreted as the time averaged firing rate of single
neurons [41]. Further, the approximation of the recurrent synaptic inputs in
the steady state usually relies on the averaging over presynaptic ensembles of
homogeneously firing neurons. Nonetheless, the theory still works quite well,
which might be due to the fact that we are calculating the synaptic drive as
an integral over a continuum of presynaptic neurons, thereby effectively
averaging out deviations of the synaptic drive that are to be expected in
single neuron samples from such ensembles.

Although bump attractor states are interesting from a functional point of
view, they are not the only solutions to ring-like networks. As we have seen,
spiking networks can show bi-stability, in which both a stable uniform state
and a stable evoked bump-state co-exist. Here, we have neglected solutions of
our theory that converged to the uniform state (or intermediate unstable
solutions) if a bump-state at larger firing rates was also found as a
solution. While it goes beyond the scope of the current work, our theory could
be extended to also predict the dynamical stability of the states that are
found: similar to approaches in networks that admit a splitting into discrete
populations [28] one could measure the magnitude of perturbations to the
firing rates at the points $\theta_{i}$ resulting from perturbations to the
parametrization: the general perturbation
$g(g_{0},g_{1},g_{\sigma},g_{r})\rightarrow g(g_{0}+\delta g_{0},g_{1}+\delta
g_{1},g_{\sigma}+\delta g_{\sigma},g_{r}+\delta g_{r})$ would translate into
both perturbations of the rates
$g(\theta_{i})=\nu_{0}(\theta_{i})+\delta\nu(\theta_{i})$ and the prediction
$F(\text{input}_{\theta_{i}}\left[g_{0}\right],\bar{V}(\theta_{i}))+\delta
F(\theta_{i})$. Stability can then be determined by comparing the scale of the
input perturbations $\delta\nu(\theta_{i})$ to the predicted output
perturbations $\delta F(\theta_{i})$ [28]. It is worth noting, that
calculating such linear perturbations will involve the derivative of the
synaptic activations $\psi$ (cf. (LABEL:spike-selfcons-rate)), which will also
have to be estimated numerically in cases where no analytical formula for
$\psi$ can be found.

Our choice of parametrization of the firing rate profile was heuristic, guided
by the shapes observed in numerical simulations. The framework presented here,
though, could be used with any other family of functions parametrized by a
small number of parameters, which should be adapted to the shapes to be
approximated. For example, multi-modal ring-attractor profiles resulting from
narrower connectivities (see e.g. [42], or [32] for a spiking network similar
to the one investigated here) can not be approximated by the unimodal family
chosen here. Since the topology of ring-attractors is periodic, a natural
candidate for such a generalization would be the family of finite Fourier
series. However, the nonzero frequency components necessary to faithfully
approximate shapes that deviate far from simple (e.g. cosine-shaped) unimodal
distributions might require a large number of Fourier coefficients for
parametrization.

In this report we mostly chose as many positions $\theta_{i}$ along the
attractor manifold as there are free parameters of the profile $g$ to be
constrained. In principle, the number of free parameters can be chosen
independently of the number of positions, by performing numerical optimization
on a dimension-agnostic sum of squared errors. We have shown that matching the
error dimension to the number of parameters permits using efficient
optimization methods that significantly speed up the optimization. However, we
have also investigated under-determined (see 6B) and over-determined (see
5A-C) systems, which also converged to similar solutions. Finally, when
optimizing the network parameters for desired spatial profiles (3 and 6),
choosing optimization goals outside the space of possible solutions of the
network dynamics did not allow the procedure to converge. Thus, our approach
could be used to estimate the boundaries of the solution space for a given
neural network, by starting the optimization at a known solution and varying
the shape parameters until convergence fails.

###  4 Methods

####  4.1 Rate network model

We study a network of $N=100$ recurrently connected rate neurons indexed by
$i\in\\{0,...,N-1\\}$, where each neuron $i$ is described by a single variable
$\nu_{i}(t)$ which denotes the firing rate of the neuron [43].

The neuronal firing rate is given by a nonlinear input-output function $F$ of
a synaptic variable $s_{i}$ (assuming that the membrane time constant is
considerably faster than the synaptic variable):

$\displaystyle\nu_{i}(t)$
$\displaystyle=F(s_{i}(t))=\frac{\nu_{\text{max}}}{2}\left[1+\tanh\left(\frac{s_{i}(t)}{s_{0}}\right)\right].$
(13)

Here, $\nu_{\text{max}}=50\text{Hz}$ sets the maximal firing frequency, and
$s_{0}=1$ is the (dimensionless) scale of the synaptic input.

The input to neuron $i$ is mediated through the synaptic variable $s_{i}$:

$\displaystyle\dot{s}_{i}(t)$
$\displaystyle=-\frac{s_{i}(t)}{\tau_{s}}+\sum_{j=0}^{N-1}w_{ij}\nu_{j}(t),$
(14)

where $\dot{s}(t)=\frac{d}{dt}s(t)$ denotes the temporal derivative,
$\tau_{s}=100ms$ is the synaptic time constant, and $w_{ij}$ are the recurrent
connection weights (see below), and $0\leq i\leq N-1$.

Neurons are organized at circular positions
$\theta_{i}=i\cdot\frac{2\pi}{N}-\pi\in[-\pi,\pi)$ with identified boundaries,
such that neuron $0$ is the direct neighbor of neuron $N-1$. The recurrent
connections depend only on the distance between neurons in the resulting
angular space: the connection $w_{ij}$ from neuron $j$ to neuron $i$ is given
by a generalized Gaussian function, with 4 free parameters controlling its
shape:

$\displaystyle w_{ij}$
$\displaystyle=w(\theta_{i}-\theta_{j})=\frac{1}{N}\left(w_{0}+w_{1}\exp\left[-\left(\frac{\left|\min\left(\left|\theta_{i}-\theta_{j}\right|,2\pi-\left|\theta_{i}-\theta_{j}\right|\right)\right|}{w_{\sigma}}\right)^{w_{r}}\right]\right).$
(15)

The parameters $\\{w_{0},w_{1},w_{\sigma},w_{r}\\}$ used for the networks of 2
and 3 are given in 2.

####  4.2 Spiking network model

Spiking simulations are based on a reimplementation of a popular ring-
attractor model of visuospatial working memory [13] in the NEST simulator
[44]. Parameters were modified from the original publication to produce the
results shown in 4, 5, and 6 (see 1 for parameter values). For completeness we
restate the definition of the model here.

#####  4.2.1 Neuron model

Neurons are modeled by leaky integrate-and-fire dynamics with conductance
based synaptic transmission [13, 27]. The network consists of recurrently
connected populations of $N_{E}$ excitatory and $N_{I}$ inhibitory neurons,
both additionally receiving external spiking input with spike times generated
by $N_{\text{ext}}$ independent, homogeneous Poisson processes, with mean
rates $\nu_{ext}$. Following [13], we assume that external excitatory inputs
are mediated by fast AMPA receptors, while recurrent excitatory currents are
mediated only by slower NMDA channels.

The neuronal dynamics for neurons in both excitatory and inhibitory
populations are governed by the following system of differential equations
indexed by $i\in\\{0,...,N_{E/I}-1\\}$ (with different sets of parameters for
each population):

$\displaystyle C_{m}\dot{V}_{i}(t)$ $\displaystyle=$ $\displaystyle-
I_{i}^{L}(t)-I_{i}^{\text{Ext}}(t)-I_{i}^{\text{I}}(t)-I_{i}^{\text{E}}(t),$
(16) $\displaystyle I_{i}^{P}$ $\displaystyle=$ $\displaystyle
g_{P}\,s_{i}^{P}(V_{i}(t),t)\,\left(V_{i}(t)-V_{P}\right),$

where $P\in\\{\text{\text{L,Ext,I,E}}\\}$. Here, $C_{\text{m}}$ is the
membrane capacitance and $V_{\text{L}},V_{\text{E}},V_{\text{I}}$ are the
reversal potentials for leak, excitatory currents, and inhibitory currents,
respectively. The parameters $g_{P}$ for $P\in\\{\text{\text{L,Ext,I,E}}\\}$
are fixed scales for leak (L), external input (Ext) and recurrent excitatory
(E) and inhibitory (I) synaptic conductances, which are dynamically gated by
the (possibly voltage dependent) gating variables $s_{i}^{P}(V,t)$. In the
main text we refer to the conductance scales of excitatory neurons by the
“strength of synaptic connections” $g_{\text{EE}}=g_{\text{E}}$ and
$g_{\text{EI}}=g_{\text{I}}$. Similarly, for inhibitory neurons we refer to
the conductance scales by the “strengths” $g_{\text{IE}}=g_{\text{E}}$ and
$g_{\text{II}}=g_{\text{I}}$. The gating variables $s_{i}^{P}$ are described
in detail below, however we set the leak conductance gating variable to
$s_{i}^{L}=1$.

The model neuron dynamics ((LABEL:base_deqs)) are integrated until their
voltage reaches a threshold $V_{\text{thr}}$. At any such time, the respective
neuron emits a spike and its membrane potential is reset to the value
$V_{\text{res}}$. After each spike, voltages are clamped to $V_{\text{res}}$
for a refractory period of $\tau_{\text{ref}}$ (see 1 for parameter values).

#####  4.2.2 Synaptic gating variables

The synaptic gating variables $s_{i}^{P}(t)$ for
$P\in\\{\text{\text{Ext,I}}\\}$ for external and inhibitory currents are
exponential traces of the firing times $t_{j}$ of all presynaptic neurons $j$:

$\dot{s}_{i}^{P}(t)=-\frac{s_{i}^{P}(t)}{\tau_{P}}+\sum_{j\in\text{pre}(P)}\sum_{t_{j}}\delta\left(t-t_{j}\right),$
(17)

where the sum runs over all neurons presynaptic to the neuron $i$ regarding
the connection $P$.

For the recurrent excitatory gating variables $s_{i}^{\text{E}}$ a nonlinear
NMDA model is used [45]. This model has second order kinetics for NMDA channel
activation [25], which result in a saturation of channels. Together with a
voltage dependence $\text{Mg}(V_{i})$ of the conductance (due to the release
of the $\text{Mg}^{2+}$ block, see [46]) this yields the following dynamics:

$\displaystyle s_{i}^{\text{E}}(V,t)$ $\displaystyle=$
$\displaystyle\text{Mg}(V_{i})\sum_{j=1}^{N_{E}}w_{ij}^{\text{E}}y_{j}(t),$
(18) $\displaystyle\dot{y}_{j}$ $\displaystyle=$
$\displaystyle-\frac{y_{j}}{\tau_{\text{E}}}+\alpha
x_{j}(t)\left(1-y_{j}\right),$ (19) $\displaystyle\dot{x}_{j}$
$\displaystyle=$
$\displaystyle-\frac{x_{j}}{\tau_{\text{E,rise}}}+\sum_{t_{j}}\delta\left(t-t_{j}\right),$
(20) $\displaystyle\text{Mg}\left(V\right)$ $\displaystyle=$
$\displaystyle\frac{1}{1+\gamma\exp\left(-\beta V\right)}.$ (21)

See 1 for parameter values used in simulations.

#####  4.2.3 Network connectivity

All connections except for the recurrent excitatory connections are all-to-all
and uniform. The recurrent excitatory connections are chosen to be distance-
dependent. As in the rate model, each neuron of the excitatory population with
index $i\in\\{0,...,N_{E}-1\\}$ is assigned an angular position
$\theta_{i}=i\cdot\frac{2\pi}{N_{E}}-\pi\in[-\pi,\pi)$. Recurrent excitatory
NMDA connections $w_{ij}^{E}$ from neuron $j$ to neuron $i$ are then given by
the Gaussian function $w^{E}(\theta)$:

$\displaystyle w_{ij}^{E}$
$\displaystyle=w^{E}(\theta_{i}-\theta_{j})=w_{0}+\left(w_{+}-w_{0}\right)\exp\left(-\left[\min\left(\left|\theta_{i}-\theta_{j}\right|,2\pi-\left|\theta_{i}-\theta_{j}\right|\right)\right]^{2}\frac{1}{2\sigma_{w}^{2}}\right).$

Additionally, for each neuron we keep the integral over all recurrent
connection weights normalized, resulting in the normalization condition
$\frac{1}{2\pi}\int_{-\pi}^{\pi}d\varphi w^{E}(\varphi)=1.$ This normalization
ensures that varying the maximum weight $w_{+}$ will not change the total
recurrent excitatory input if all excitatory neurons fire at the same rate.
Here, we choose $w_{+}$ as a free parameter and constrain the baseline
connection weight to

$\displaystyle w_{0}$
$\displaystyle=\frac{w_{+}\sigma_{w}\text{erf}\left(\frac{\pi}{\sqrt{2}\sigma_{w}}\right)-\sqrt{2\pi}}{\sigma_{w}\text{erf}\left(\frac{\pi}{\sqrt{2}\sigma_{w}}\right)-\sqrt{2\pi}}.$

####  4.3 Self-consistent equations

#####  4.3.1 Placement of sampling points

Self consistent equations ((LABEL:mfeq-selfconsist)) are constructed for both
rate-based and spiking neurons by using the low-dimensional parametrization
(LABEL:g-gen-gaussian) described in the main text. As mentioned there, we
choose the top of the firing rate profile $\theta_{1}=0$, as well as the
lowest point $\theta_{4}=\pi$. For the intermediate points $0<\theta_{i}<\pi$
for $i\in\\{2,3\\}$, we sample the firing rate profile by inverting the
function $g$ to give a sample at a desired height
$h_{i}=g_{0}+a_{i}\left(g_{1}-g_{0}\right)$ with $0<a_{i}<1$. This yields a
relation for the position which depends on the shape parameters $g_{r}$ and
$g_{\sigma}$:

$\theta_{i}=-g_{\sigma}\log\left(a_{i}\right)^{\frac{1}{g_{r}}}.$

For all figures except for 5 and 6, the intermediate points were chosen by
setting $a_{2}=0.2$ and $a_{3}=0.8$, although we show (2 and 3) that the exact
choice of points affects the solutions only slightly.

In 5 we iterate through even numbers $p\geq 4$ of sampling points. As before,
we choose $\theta_{1}=0$, as well as the lowest point $\theta_{p}=\pi$.
Generalizing the placement of $4$ points described above, the remaining $p-2$
points were chosen as $\theta_{k}=\nicefrac{{0.4}}{{k}}$ for
$1<k\leq\nicefrac{{p}}{{2}}$, and $\theta_{k}=1-\nicefrac{{0.4}}{{k}}$ for
$\nicefrac{{p}}{{2}}<k\leq p$. For the optimization of spiking network
parameters shown in 6, we chose $7$ sampling points: we first chose $p=6$
points by the scheme just described, then added the point $\theta_{7}=0.5$.

#####  4.3.2 Derivation of input-output functions for the rate network

For the rate network, we set (LABEL:rate-input) to zero and solve for the
steady-state input $s_{i}$, which yields

$\displaystyle s_{i}$
$\displaystyle=\tau_{s}\sum_{j=0}^{N-1}w_{ij}\nu_{j}=\tau_{s}\sum_{j=0}^{N-1}w_{ij}g(\theta_{j})$
$\displaystyle\approx\frac{\tau_{s}N}{2\pi}\int_{-\pi}^{\pi}d\varphi
w(\theta_{i}-\varphi)g(\varphi).$

Here, we have replaced the activity of neurons in the network by our
parametrization $g$. In the second line we approximated the summation
$\frac{1}{N}\sum_{j=0}^{N-1}$ by the integral
$\frac{1}{2\pi}\int_{-\pi}^{\pi}d\varphi$ and exploited that the connectivity
$w$ is only dependent on the angular distance between the neuron $i$ (at
position $\theta_{i}$) and neurons $j$ (at varying positions
$\theta_{j}=\varphi$), to replace $w_{ij}\rightarrow w(\theta_{i}-\varphi)$.
We use this steady state input in (LABEL:rate-neurons) to arrive at
(LABEL:rate-error).

#####  4.3.3 Derivation of input-output functions for the spiking network

In the rate model presented in 4.1 the firing rates are given by (LABEL:rate-
neurons). In the spiking network, we have to approximate the expected firing
rates of neurons. To this end, we first replace the synaptic activation
variables $s^{P}(V,t)$ for $P\in\\{I,E,\text{ext}\\}$ by their expectation
values under Poisson input. For the linear synapses this yields

$\displaystyle\left\langle s^{\text{ext}}(t)\right\rangle_{\text{Poisson}}=$
$\displaystyle\tau_{\text{ext}}\nu_{\text{pre}},$ $\displaystyle\left\langle
s^{\text{I}}(t)\right\rangle_{\text{Poisson}}=$
$\displaystyle\tau_{\text{I}}\nu_{I}.$

The nonlinear synaptic activation of NMDA synapses under stimulation with
Poisson processes at rates $\nu_{j}$ was estimated by simulating Eqs.
(19)-(20) under varying presynaptic firing rates and fitting an interpolating
function to the temporal means of the synaptic activation
$\psi(\nu_{j})\equiv\langle y_{j}\rangle_{t}$. An analytical approximation of
the function $\psi(\nu_{j})$ was stated in [27, p. 80]. We instead chose to
numerically fit this function to simulated data, since for higher firing rates
the analytical approximation tended to over-estimate the synaptic activations.

We then define the expected recurrent excitatory input, assuming presynaptic
Poisson firing, by

$\displaystyle J_{i}$
$\displaystyle\equiv\frac{1}{N_{E}}\sum_{i=0}^{N_{E}-1}w_{ij}^{E}\psi(\nu_{j}).$
(22)

Following [27], we linearize the voltage dependence (LABEL:Mg-NMDA) at the
mean voltage $\left\langle V\right\rangle$ and reduce the differential
equations of (LABEL:base_deqs) to dimensionless form. The resulting
expressions depend only on the mean firing rates and mean voltages of
excitatory and inhibitory neurons (see A of the Appendix for the detailed
expressions and derivations):

$\displaystyle\tau_{i}\dot{V}_{i}$ $\displaystyle=$
$\displaystyle-\left(V_{i}-V_{L}\right)+\mu_{i}+\sigma_{i}\sqrt{\tau_{i}}\eta_{i}(t)$
(23) $\displaystyle\mu_{i}$ $\displaystyle=$
$\displaystyle\mu_{i}\left(J_{i},\nu_{I},\nu_{\text{ext}},\left\langle
V_{i}\right\rangle\right)$ $\displaystyle\sigma_{i}$ $\displaystyle=$
$\displaystyle\frac{g_{\text{ext}}}{C_{m}}\left(\left\langle V\right\rangle-
V_{E}\right)\tau_{\text{ext}}\sqrt{\tau_{i}N_{\text{ext}}\nu_{\text{ext}}}.$
$\displaystyle\tau_{i}$ $\displaystyle=$
$\displaystyle\tau_{i}\left(J_{i},\nu_{I},\nu_{\text{ext}},\left\langle
V_{i}\right\rangle\right)$ $\displaystyle\left\langle\eta_{i}(t)\right\rangle$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\left\langle\eta_{i}(t)\eta_{i}(t^{\prime})\right\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{\tau_{\text{ext}}}\exp(-\frac{\left|t-t^{\prime}\right|}{\tau_{\text{ext}}})$

Here, $\mu_{i}$ is the bias of the membrane potential due to synaptic inputs,
and $\sigma_{i}$ measures the scale of fluctuations in the membrane potential
due to random spike arrival approximated by the Gaussian process $\eta_{i}$.
Due to active synaptic conductances, the effective membrane time constant
$\tau_{i}$ is decreased from the intrinsic membrane time-constant
$\nicefrac{{C_{m}}}{{g_{L}}}$ – its value thus depends on all presynaptic
firing rates (and the mean voltage, see A of the Appendix).

The prediction $F$ of the mean firing rates and $\left\langle
V_{i}\right\rangle$ of mean voltages of populations of neurons governed by
this type of differential equation can be well approximated by [27] (see also
the published corrections in [47]):

$\displaystyle\phi\left[\mu_{i},\sigma_{i},\tau_{i}\right]$ $\displaystyle=$
$\displaystyle\left(\tau_{\text{ref}}+\sqrt{\pi}\tau_{i}\int_{\beta(\mu_{i},\sigma_{i})}^{\alpha(\mu_{i},\sigma_{i})}du\exp(u^{2})\left[1+\text{erf}\left(u\right)\right]\right)^{-1},$
(24) $\displaystyle\alpha(\mu_{i},\sigma_{i})$ $\displaystyle=$
$\displaystyle\frac{V_{\text{reset}}-V_{L}-\mu_{i}}{\sigma_{i}}\left(1+\frac{\tau_{\text{ext}}}{2\tau_{i}}\right)+1.03\sqrt{\frac{\tau_{\text{ext}}}{\tau_{i}}}-\frac{\tau_{\text{ext}}}{\tau_{i}},$
(25) $\displaystyle\beta(\mu_{i},\sigma_{i})$ $\displaystyle=$
$\displaystyle\frac{V_{\text{reset}}-V_{L}-\mu_{i}}{\sigma_{i}},$ (26)
$\displaystyle\left\langle V_{i}\right\rangle$ $\displaystyle=$
$\displaystyle\mu_{i}+V_{L}-\left(V_{\text{thr}}-V_{\text{reset}}\right)\phi\left[\mu_{i},\sigma_{i},\tau_{i}\right]\tau_{i}.$
(27)

As in the rate model, we first replace the network activity $\nu_{j}$ on the
right hand side of (LABEL:total-input) by our parametrization $g(\theta_{j})$.
We then approximate the summation $\frac{1}{N_{E}}\sum_{j=0}^{N_{E}-1}$ with
an integral $\frac{1}{2\pi}\int_{-\pi}^{\pi}d\varphi$, and replace the
connectivity by its continuous equivalent $w_{ij}\rightarrow
w(\theta_{i}-\varphi)$ to arrive at:

$\displaystyle J_{i}$
$\displaystyle\approx\frac{1}{2\pi}\int_{-\pi}^{\pi}d\varphi
w^{E}\left(\theta_{i}-\varphi\right)\psi(g\left(\varphi\right))$
$\displaystyle\equiv\frac{1}{2\pi}\text{input}_{\theta_{i}}\left[g\right].$

We then substitute this relation in Eqs. (24) and (27) to arrive at

$\displaystyle\mu_{i}$
$\displaystyle=\mu_{i}\left(\text{input}_{\theta_{i}}\left[g\right],\nu_{I},\nu_{\text{ext}},\left\langle
V_{i}\right\rangle\right),$ $\displaystyle\tau_{i}$
$\displaystyle=\tau_{i}\left(\text{input}_{\theta_{i}}\left[g\right],\nu_{I},\nu_{\text{ext}},\left\langle
V_{i}\right\rangle\right),$ $\displaystyle
F(\text{input}_{\theta_{i}}\left[g\right],\nu_{I},\left\langle
V_{i}\right\rangle)$
$\displaystyle\equiv\phi\left[\mu_{i},\sigma_{i},\tau_{i}\right],$
$\displaystyle G(\text{input}_{\theta_{i}}\left[g\right],\nu_{I},\left\langle
V_{i}\right\rangle)$
$\displaystyle\equiv\mu_{i}+V_{L}-\left(V_{\text{thr}}-V_{\text{reset}}\right)g_{i}(\theta_{i})\tau_{i},$
(28)

which defines Eqs. (5) and (6) of the main text.

#####  4.3.4 Optimization of self-consistent equations

For each point $\theta_{i}$ that we choose to sample from the excitatory
population, the theory of 4.3.3 yields $2$ constraining Equations (24) and
(27). The inhibitory population, being homogeneous and unstructured, yields
$2$ equations, for the 2 free variables $\nu_{I}$ and $\left\langle
V_{i}\right\rangle$. Since we choose a low-dimensional parametrization for the
excitatory population, the number of free variables increases only by $1$ (the
mean voltage $\left\langle V_{i}\right\rangle$) for each point $\theta_{i}$
that we choose to evaluate, while yielding the same $2$ constraining
equations. This allows us to choose at minimum $4$ evaluation points to
constrain the $4$ free parameters of the parametrization (see 3 for a
listing).

The errors $\text{Err}_{i}$ (and $\text{Err}_{I}$, for spiking networks)
between firing rate predictions and the firing rate parametrization are
numerically minimized using methods provided in the _Scipy_ package [48]. In
particular, if the dimension of the error function matches the number of
parameters, we are able to use the efficient _optimize.root_ solver (_Root_ in
the main text), which applies a modified version of the Powell hybrid method
[49], but does not provide constraints on valid parameter regions. Here, we
implemented artificial constraints by returning a high error for dimensions
that leave the bounded region. The same optimization results were achieved by
using the slower _optimize.minimize_ method, which allows optimization (of the
sum of squared errors $\text{SSE}=\sum_{i}\text{Err}_{i}^{2}$, or
$\text{SSE}=\sum_{i}\text{Err}_{i}^{2}+\text{Err}_{I}^{2}$ for spiking
networks) in constrained parameter regions via the L-BFGS-B [50] and __ SLSQP
[51]. For spiking networks, we normalized firing rate errors by the firing
rate $\nu_{\text{max}}=100Hz$ and voltage differences by the voltage range
$V_{\text{thres}}-V_{\text{reset}}$, to ensure comparable contributions to the
SSE for variables with different dimensions.

For the optimization results of 6 we chose $7$ sampling points (see 4.3.1 for
details), which yielded $16$ errors, including those of the inhibitory
population. These were used to optimize $17$ free parameters using the SLSQP
algorithm (see 3 for a listing). We also tried using $8$ sampling points,
which brings both the number of equations and free parameters up to $20$ –
this yielded similar results at increased processing time (the possibly faster
_Root_ solver failed to converge most of the time).

Wall clock times for error functions in 5A were measured on a single core of a
MacBook Pro with 2,6 GHz Intel Core i5 processor, using the Python benchmark
_timeit.timeit_ (minimum wall clock time of $100$ repetitions). We first
measured average time for evaluation of a single error $\text{Err}_{i}$ of
(LABEL:spike-selfcons-err), which evaluated to $t_{E}=4.59ms$ (100 repetitions
of 10 executions). For a single evaluation of the inhibitory error
$\text{Err}_{I}$ of (LABEL:spike-selfcons-err-I) we found $t_{I}=0.98ms$ (the
numerical integration performed in the calculation of
$\text{input}_{\text{I}}=\frac{1}{2\pi}\int
d\varphi\psi\left(g\left(\varphi\right)\right)$ is faster). The wall clock
time $T$ (see 5A, right axis) for a given number $n$ of error vector
evaluations on $p$ points was then calculated by $T=n\left(p\cdot
t_{E}+t_{I}\right)$.

####  4.4 Spiking simulations

All network simulations where performed in the NEST simulator [44] using
fourth-order Runge-Kutta integration as implemented in the GSL package [52].
For the simulation results shown in 4 and 6B, networks underwent a transient
initial period of $t_{\text{initial}}$. Neurons centered at a position of
$\theta=0$ then received a short and strong excitatory input mediated by
additional Poisson firing onto AMPA receptors ($500ms,2\text{kHz}$) with
connections scaled down by a factor of $g_{\text{signal}}=0.5$. The external
input ceased at $t=t_{\text{off}}$.

Simulations were run until $t=t_{\text{max }}$ and spikes were recorded and
converted to firing rates by spike counts in a $75ms$ window shifted at a time
resolution of $1ms$. For every time step, the firing rates across the whole
population were then rectified (by measuring the phase of the first spatial
Fourier coefficient and setting it to $\theta=0$ by rotation of the angular
space) to center the bump of activity around the position $\varphi=0$. The
resulting centered firing rates were then sampled at an interval of $60ms$ in
the interval $\left[t_{\text{off}}+500ms,t_{\text{max}}\right]$ for 5
repetitions of the network simulation with the same microscopic parameters. In
4B times were: $t_{\text{initial}}=1.8s,t_{\text{off }}=2s,t_{\text{max}}=5s$.
In 4D,E we chose: $t_{\text{initial}}=0.5s,t_{\text{off
}}=1s,t_{\text{max}}=3s$. For 6B: $t_{\text{initial}}=0.5s,t_{\text{off
}}=1s,t_{\text{max}}=5s$.

###  5 Acknowledgements

### ?refname?

  * [1] Alexander Kevin Seeholzer.  Continuous Attractor Working Memory and Provenance of Channel Models.  Doctoral thesis, EPFL, Lausanne, 2017.  doi:10.5075/epfl-thesis-7845.
  * [2] Donald Olding Hebb.  The Organization of Behavior: A Neuropsychological Approach.  John Wiley & Sons, New York, NY, USA, 1949.
  * [3] S. J. Martin, P. D. Grimwood, and R. G. Morris.  Synaptic plasticity and memory: An evaluation of the hypothesis.  Annual Review of Neuroscience, 23:649–711, 2000.  doi:10.1146/annurev.neuro.23.1.649.
  * [4] Tomonori Takeuchi, Adrian J. Duszkiewicz, and Richard G. M. Morris.  The synaptic plasticity and memory hypothesis: Encoding, storage and persistence.  Philosophical Transactions of the Royal Society B: Biological Sciences, 369(1633), January 2014.  doi:10.1098/rstb.2013.0288.
  * [5] P. S Goldman-Rakic.  Cellular basis of working memory.  Neuron, 14(3):477–485, March 1995.  doi:10.1016/0896-6273(95)90304-6.
  * [6] Clayton E. Curtis and Mark D’Esposito.  Persistent activity in the prefrontal cortex during working memory.  Trends in Cognitive Sciences, 7(9):415–423, September 2003.  doi:10.1016/S1364-6613(03)00197-9.
  * [7] Rishidev Chaudhuri and Ila Fiete.  Computational principles of memory.  Nature Neuroscience, 19(3):394–403, March 2016.  doi:10.1038/nn.4237.
  * [8] H R Wilson and J D Cowan.  A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue.  Kybernetik, 13(2):55–80, September 1973.
  * [9] Shun-ichi Amari.  Dynamics of Pattern Formation in Lateral-Inhibition Type Neural Fields.  Biological cybernetics, 87:77–87, 1977.
  * [10] R Ben-Yishai, R L Bar-Or, and H Sompolinsky.  Theory of orientation tuning in visual cortex.  Proceedings of the National Academy of Sciences, 92(9):3844–8, April 1995.
  * [11] S Funahashi, C J Bruce, and P S Goldman-Rakic.  Mnemonic coding of visual space in the monkey’s dorsolateral prefrontal cortex.  Journal of Neurophysiology, pages 331–349, 1989.
  * [12] K. Zhang.  Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: A theory.  Journal of Neuroscience, 16(6):2112–2126, March 1996.
  * [13] Albert Compte, Nicolas Brunel, P.S. Goldman-Rakic, and X.J. Wang.  Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model.  Cerebral Cortex, 10:910–923, 2000.
  * [14] Vladimir Itskov, Carina Curto, Eva Pastalkova, and György Buzsáki.  Cell assembly sequences arising from spike threshold adaptation keep track of time in the hippocampus.  The Journal of Neuroscience, 31(8):2828–34, February 2011.  doi:10.1523/JNEUROSCI.3773-10.2011.
  * [15] KiJung Yoon, Michael A. Buice, Caswell Barry, Robin Hayman, Neil Burgess, and Ila R. Fiete.  Specific evidence of low-dimensional continuous attractor dynamics in grid cells.  Nature Neuroscience, 16(8):1077–1084, August 2013.  doi:10.1038/nn.3450.
  * [16] Paul C Bressloff and Jack D Cowan.  A spherical model for orientation and spatial-frequency tuning in a cortical hypercolumn.  Phil. Trans. R. Soc. Lond. B, 358(1438):1643–67, October 2003.  doi:10.1098/rstb.2002.1109.
  * [17] C. Fung, K. Wong, and Si Wu.  A moving bump in a continuous manifold: A comprehensive study of the tracking dynamics of continuous attractor neural networks.  Neural Computation, 22(3):752–792, 2010.
  * [18] V Itskov, D Hansel, and Misha Tsodyks.  Short-Term Facilitation may Stabilize Parametric Working Memory Trace.  Frontiers in Computational Neuroscience, 5(October):40–40, January 2011.  doi:10.3389/fncom.2011.00040.
  * [19] Carlo R. Laing.  Derivation of a neural field model from a network of theta neurons.  Physical Review E, 90(1):010901, 2014.
  * [20] Robert Rosenbaum and Brent Doiron.  Balanced Networks of Spiking Neurons with Spatially Dependent Recurrent Connections.  Physical Review X, 4(2), May 2014.  doi:10.1103/PhysRevX.4.021039.
  * [21] Alfonso Renart, Pengcheng Song, and Xiao-Jing Wang.  Robust spatial working memory through homeostatic synaptic scaling in heterogeneous cortical networks.  Neuron, 38(3):473–85, May 2003.
  * [22] Romain Brette and Wulfram Gerstner.  Adaptive exponential integrate-and-fire model as an effective description of neuronal activity.  Journal of Neurophysiology, 94(5):3637–42, November 2005.  doi:10.1152/jn.00686.2005.
  * [23] James P. Roach, Eshel Ben-Jacob, Leonard M. Sander, and Michal R. Zochowski.  Formation and Dynamics of Waves in a Cortical Model of Cholinergic Modulation.  PLoS Comput Biol, 11(8):e1004449, August 2015.  doi:10.1371/journal.pcbi.1004449.
  * [24] Robert S Zucker and Wade G Regehr.  Short-term synaptic plasticity.  Annual Review of Physiology, 64:355–405, January 2002.  doi:10.1146/annurev.physiol.64.092501.114547.
  * [25] a Destexhe, ZF Mainen, and TJ Sejnowski.  Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism.  Journal of Computational Neuroscience, 1(3):195–230, August 1994\.
  * [26] M Spiridon and W Gerstner.  Effect of lateral connections on the accuracy of the population code for a network of spiking neurons.  Network, 12(4):409–21, November 2001.
  * [27] Nicolas Brunel and XJ Wang.  Effects of Neuromodulation in a Cortical Network Model of Object Working Memory Dominated by Recurrent Inhibition.  Journal of Computational Neuroscience, 11:63–85, 2001.
  * [28] Daniel J Amit and Nicolas Brunel.  Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex.  Cerebral cortex, 7:237–252, 1997.
  * [29] Nicolas Brunel and V Hakim.  Fast global oscillations in networks of integrate-and-fire neurons with low firing rates.  Neural computation, 11(7):1621–71, October 1999.
  * [30] X J Wang.  Synaptic basis of cortical persistent activity: The importance of NMDA receptors to working memory.  The Journal of Neuroscience, 19(21):9587–603, November 1999.
  * [31] John D Murray, Alan Anticevic, Mark Gancsos, Megan Ichinose, Philip R Corlett, John H Krystal, and Xiao-Jing Wang.  Linking Microcircuit Dysfunction to Cognitive Impairment: Effects of Disinhibition Associated with Schizophrenia in a Cortical Working Memory Model.  Cerebral Cortex, November 2012.  doi:10.1093/cercor/bhs370.
  * [32] Ziqiang Wei, Xiao-Jing Wang, and Da-Hui Wang.  From distributed resources to limited slots in multiple-item working memory: A spiking network model with normalization.  Journal of Neuroscience, 32(33):11228–40, August 2012.  doi:10.1523/JNEUROSCI.0735-12.2012.
  * [33] Jacinto Pereira and Xiao-Jing Wang.  A Tradeoff Between Accuracy and Flexibility in a Working Memory Circuit Endowed with Slow Feedback Mechanisms.  Cerebral Cortex, 25(10):3586, 2015.  doi:10.1093/cercor/bhu202.
  * [34] Rita Almeida, João Barbosa, and Albert Compte.  Neural circuit basis of visuo-spatial working memory precision: A computational and behavioral study.  Journal of Neurophysiology, 114(3):1806–1818, September 2015.  doi:10.1152/jn.00362.2015.
  * [35] Steven Strogatz.  Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering.  Westview Press, 1st edition, 2000.
  * [36] Oscar Marin.  Interneuron dysfunction in psychiatric disorders.  Nature Reviews Neuroscience, 13(2):107–120, February 2012.  doi:10.1038/nrn3155.
  * [37] David Hansel and German Mato.  Short-term plasticity explains irregular persistent activity in working memory tasks.  The Journal of Neuroscience, 33(1):133–49, January 2013.  doi:10.1523/JNEUROSCI.3455-12.2013.
  * [38] M Tsodyks.  Attractor neural network models of spatial maps in hippocampus.  Hippocampus, 9(4):481–9, January 1999.  doi:10.1002/(SICI)1098-1063(1999)9:4<481::AID-HIPO14>3.0.CO;2-S.
  * [39] A Samsonovich and B L McNaughton.  Path integration and cognitive mapping in a continuous attractor neural network model.  Journal of Neuroscience, 17(15):5900–20, August 1997.
  * [40] Gianluigi Mongillo, O Barak, and Misha Tsodyks.  Synaptic theory of working memory.  Science, 319(5869):1543–1546, March 2008.  doi:10.1126/science.1150769.
  * [41] N Brunel and S Sergi.  Firing frequency of leaky intergrate-and-fire neurons with synaptic current dynamics.  Journal of Theoretical Biology, 195(1):87–95, November 1998.  doi:10.1006/jtbi.1998.0782.
  * [42] Carlo R. Laing, William C. Troy, Boris Gutkin, and G. Bard Ermentrout.  Multiple Bumps in a Neuronal Model of Working Memory.  SIAM Journal on Applied Mathematics, 63(1):62–97, January 2002\.  doi:10.1137/S0036139901389495.
  * [43] Kenneth D Miller and Francesco Fumarola.  Mathematical equivalence of two common forms of firing rate models of neural networks.  Neural Computation, 24(1):25–31, January 2012.  doi:10.1162/NECO_a_00221.
  * [44] Markus Diesmann and Marc-Oliver Gewaltig.  NEST (NEural Simulation Tool).  Scholarpedia, 2(4):1430–1430, 2007.  doi:10.4249/scholarpedia.1430.
  * [45] X J Wang.  Synaptic reverberation underlying mnemonic persistent activity.  Trends in Neurosciences, 24(8):455–63, August 2001.
  * [46] C. E. Jahr and C. F. Stevens.  A quantitative description of NMDA receptor-channel kinetic behavior.  Journal of Neuroscience, 10(6):1830–1837, June 1990.
  * [47] N Brunel and XJ Wang.  Corrections for Effects of Neuromodulation in a Cortical Network Model of Object Working Memory Dominated by Recurrent Inhibition -.  https://galton.uchicago.edu/~nbrunel/pdfs/corrBW01.pdf, 2001.
  * [48] Travis E. Oliphant.  Python for Scientific Computing.  Computing in Science & Engineering, 9(3):10–20, May 2007.  doi:10.1109/MCSE.2007.58.
  * [49] Michael JD Powell.  A hybrid method for nonlinear equations.  Numerical Methods for Nonlinear Algebraic Equations, 7:87–114, 1970\.
  * [50] R. Byrd, P. Lu, J. Nocedal, and C. Zhu.  A Limited Memory Algorithm for Bound Constrained Optimization.  SIAM Journal on Scientific Computing, 16(5):1190–1208, September 1995.  doi:10.1137/0916069.
  * [51] Dieter Kraft.  A software package for sequential quadratic programming.  Forschungsbericht - Deutsche Forschungs- und Versuchsanstalt fur Luft- und Raumfahrt, 1988.
  * [52] Mark Galassi, Jim Davies, James Theiler, Brian Gough, Gerard Jungman, Michael Booth, and Fabrice Rossi.  GNU Scientific Library Reference Manual.  3rd edition, 2009.  ISBN 0954612078.

## Appendix

###  ?appendixname? A Detailed derivation of dimensionless voltage equations

In this section we give details on the derivation of (LABEL:reduced-voltage-
deq) as well as the resulting full expressions. This closely follows [27, pp.
79–81], while keeping a slightly simplified notation.

We first replace all synaptic activations in (LABEL:base_deqs) by their
expected values under Poisson input, which also introduces the expected
recurrent excitatory input $J_{i}$ (cf. (LABEL:total-input)):

$\displaystyle s_{i}^{\text{ext}}(t)$ $\displaystyle\rightarrow
N_{\text{ext}}\tau_{\text{ext}}\nu_{\text{pre}}+\Delta_{S,ext},$
$\displaystyle s_{i}^{\text{I}}(t)$ $\displaystyle\rightarrow
N_{I}\tau_{\text{I}}\nu_{I},$ $\displaystyle s_{i}^{E}(t)$
$\displaystyle\rightarrow\text{Mg}(V_{i})\sum_{i=0}^{N_{E}-1}w_{ij}^{E}\psi(\nu_{j})=\text{Mg}(V_{i})N_{E}J_{i}.$

Here, $\Delta_{S,ext}$ represents fluctuations of around the mean of
$s_{i}^{\text{ext}}$ due to random spike arrival at fast AMPA synapses. Since
the synaptic timescales (GABA, NMDA) of the other synaptic activations are
much longer, these fluctuations can be neglected. We then rearrange
(LABEL:base_deqs) to dimensionless form, which yields:

$\displaystyle\frac{C_{m}}{g_{L}}\dot{V}_{i}$ $\displaystyle=$
$\displaystyle-\left(V_{i}-V_{L}\right)\left[1+T_{I}\nu_{I}+T_{\text{ext}}\nu_{\text{ext}}+\frac{g_{E}}{g_{L}}\text{Mg}(V_{i})N_{E}J_{i}\right]$
$\displaystyle+\left(V_{I}-V_{L}\right)T_{I}\nu_{I}+\left(V_{E}-V_{L}\right)\left[T_{\text{ext}}\nu_{\text{ext}}+\frac{g_{E}}{g_{L}}\text{Mg}(V_{i})N_{E}J_{i}\right]$
$\displaystyle+\frac{g_{\text{ext}}}{g_{L}}\left(V_{i}-V_{E}\right)\Delta_{S,\text{ext}},$

where
$T_{\text{ext}}=N_{\text{ext}}\tau_{\text{ext}}\frac{g_{\text{ext}}}{g_{L}},\,T_{I}=N_{I}\tau_{I}\frac{g_{I}}{g_{L}}$
are effective timescales of external and inhibitory input.

To get rid of the nonlinear voltage dependence of the right hand side through
$\text{Mg}(V_{i})$, we linearize this function (cf. (LABEL:Mg-NMDA)) around
the mean voltage $\left\langle V_{i}\right\rangle$:

$\displaystyle\frac{V_{i}-V_{E}}{1+\gamma\exp(-\beta V_{i})}$ $\displaystyle=$
$\displaystyle\frac{\left\langle V_{i}\right\rangle-
V_{E}}{\rho}+\left(V-\left\langle
V\right\rangle\right)\frac{\rho+\beta\left(\left\langle V_{i}\right\rangle-
V_{E}\right)\left(\rho-1\right)}{\rho^{2}},$

where $\rho=1+\gamma\exp\left(-\beta\left\langle V_{i}\right\rangle\right)$.

After replacing the voltage dependence in the fluctuation term by the mean
voltage, we arrive at

$\displaystyle\frac{C_{m}}{g_{L}}\dot{V}_{i}$ $\displaystyle=$
$\displaystyle-\left(V_{i}-V_{L}\right)\left[1+T_{I}\nu_{I}+T_{\text{ext}}\nu_{\text{ext}}+\left(\rho_{1}+\rho_{2}\right)J_{i}\right]$
$\displaystyle+\left(V_{I}-V_{L}\right)T_{I}\nu_{I}+\left(V_{E}-V_{L}\right)\left[T_{\text{ext}}\nu_{\text{ext}}+\rho_{1}J_{i}\right]$
$\displaystyle+\rho_{2}\left(\left\langle V_{i}\right\rangle-
V_{L}\right)J_{i}+\frac{g_{\text{ext}}}{g_{L}}\left(\left\langle
V_{i}\right\rangle-V_{E}\right)\Delta_{S,\text{ext}}.$ $\displaystyle\rho_{1}$
$\displaystyle=$ $\displaystyle\frac{g_{\text{E}}N_{E}}{g_{L}\rho}$
$\displaystyle\rho_{2}$ $\displaystyle=$
$\displaystyle\beta\frac{g_{\text{E}}N_{E}\left(\left\langle
V_{i}\right\rangle-V_{E}\right)\left(\rho-1\right)}{g_{L}\rho^{2}}$
$\displaystyle\rho$ $\displaystyle=$ $\displaystyle
1+\gamma\exp\left(-\beta\left\langle V_{i}\right\rangle\right).$ (29)

Finally, we replace the fluctuations $\Delta_{S,ext}$ by independent Gaussian
noise processes with zero mean $\left\langle\eta_{i}(t)\right\rangle=0$ and
simpler autocorrelation
$\left\langle\eta_{i}(t)\eta_{j}(t^{\prime})\right\rangle=\frac{1}{\tau_{\text{ext}}}\exp(-\frac{\left|t-t^{\prime}\right|}{\tau_{\text{ext}}})\delta_{ij}$,
to arrive at the full form of (LABEL:reduced-voltage-deq) in the main text:

$\displaystyle\tau_{i}\dot{V}_{i}$ $\displaystyle=$
$\displaystyle-\left(V_{i}-V_{L}\right)+\mu_{i}+\sigma_{i}\sqrt{\tau_{i}}\eta_{i}(t)$
(30) $\displaystyle S_{i}$ $\displaystyle=$ $\displaystyle
1+T_{I}\nu_{I}+T_{\text{ext}}\nu_{\text{ext}}+\left(\rho_{1}+\rho_{2}\right)J_{i}$
$\displaystyle\mu_{i}S_{i}$ $\displaystyle=$
$\displaystyle\left(V_{I}-V_{L}\right)T_{I}\nu_{I}+\left(V_{E}-V_{L}\right)T_{\text{ext}}\nu_{\text{ext}}+$
$\displaystyle\left[\rho_{1}\left(V_{E}-V_{L}\right)+\rho_{2}\left(\left\langle
V\right\rangle-V_{L}\right)\right]J_{i}$ $\displaystyle\sigma_{i}$
$\displaystyle=$ $\displaystyle\frac{g_{\text{ext}}}{C_{m}}\left(\left\langle
V\right\rangle-
V_{E}\right)\tau_{\text{ext}}\sqrt{\tau_{i}N_{\text{ext}}\nu_{\text{ext}}}.$
$\displaystyle\tau_{i}$ $\displaystyle=$
$\displaystyle\frac{C_{m}}{g_{L}S_{i}}$

Reducing the conductance based differential equation (LABEL:base_deqs) of the
main text to the simplified form (LABEL:app-full-voltage-eq), now allows us to
compute the mean firing rate as a functions of the (input-like) bias $\mu_{i}$
and fluctuation term $\sigma_{i}$, according to (LABEL:siegert) of the last
section.

Parameter name | Parameter Symbol | Units | Excitatory neurons | Inhibitory neurons
---|---|---|---|---
Neuron number | $N_{\bullet}$ | 1 | $N_{E}=800$ | $N_{I}=200$
Poisson neuron number | $N_{\text{ext}}$ | 1 | $1000$ | $1000$
Membrane capacitance | $C_{m}$ | pF | $500$ | $200$
Exc. reversal potential | $V_{E}$ | mV | $0$ | $0$
Inh. reversal potential | $V_{I}$ | mV | $-70$ | $-70$
Leak reversal potential | $V_{L}$ | mV | $-70$ | $-70$
After spike reset potential | $V_{\text{res}}$ | mV | $-60$ | $-60$
Spiking threshold | $V_{\text{thr}}$ | mV | $-50$ | $-50$
NMDA parameter 1 | $\beta$ | 1 | $0.062$ | $0.062$
NMDA parameter 2 | $\gamma$ | 1 | $1/3.57$ | $1/3.57$
NMDA rise parameter | $\alpha$ | s | 0.5 | 0.5
External conductance | $g_{ext}$ | nS | $2.08$ | $1.62$
Recurr. inh. conductance | $g_{\text{I}}$ | nS | $1.336$ | $1.024$
Recurr. exc. conductance | $g_{\text{E}}$ | nS | $0.381$ | $0.292$
Leak conductance | $g_{\text{L}}$ | nS | $25$ | $20$
External synaptic timescale | $\tau_{\text{ext}}$ | ms | $2$ | $2$
Recurr. inh. timescale | $\tau_{I}$ | ms | $10$ | $10$
Recurr. exc. timescale | $\tau_{E}$ | ms | $100$ | $100$
Recurr. exc. rise timescale | $\tau_{\text{E,rise}}$ | ms | $2$ | $2$
Membrane time constant $\frac{C_{m}}{g_{L}}$ | $\tau_{m}$ | ms | $20$ | $10$
Refractory period | $\tau_{\text{ref}}$ | ms | $2$ | $1$
Width of distance dep. weights | $\sigma_{w}$ | rad | $\frac{18\deg}{360\deg}\cdot 2\pi\approx 0.31$ | -
Frequency of Poisson neurons | $\nu_{\text{ext}}$ | Hz | $2.4Hz$ | $2.4Hz$

?tablename? 1: Parameters for spiking simulations. Parameter values are modified from [13] and [27]. System | $w_{0}$ | $w_{1}$ | $w_{\sigma}$ | $w_{r}$
---|---|---|---|---
Sys. 0 | $-0.8$ | $2.3$ | $0.9$ | $2.0$
Sys. 1 | $-1.0$ | $10.$ | $0.2$ | $2.0$
Sys. 2 | $-3.0$ | $15.0$ | $0.5$ | $2.0$
?tablename? 2: Connectivity parameters of rate models.

Figure | Fixed | Optimized | # Optimized | Points / Errors
---|---|---|---|---
4 | $\begin{array}[t]{c}w_{+},w_{\sigma}\\\ g_{\text{EE}},g_{\text{IE}},g_{\text{EI}},g_{\text{II}}\end{array}$ | $\begin{array}[t]{c}g_{0},g_{1},g_{\sigma},g_{r}\\\ \left\langle V\right\rangle_{\theta_{1}},\dots,\left\langle V\right\rangle_{\theta_{4}}\\\ \nu_{I},\left\langle V\right\rangle_{I}\end{array}$ | 10 | 4 / 10
6 red | $\begin{array}[t]{c}g_{1}=50Hz\\\ g_{\sigma}=0.6\\\ \\\ \\\ \end{array}$ | $\begin{array}[t]{c}g_{0},g_{r}\\\ w_{+}=4.100,w_{\sigma}=0.1899\\\ g_{\text{EE}}=0.3923,g_{\text{IE}}=0.3958\\\ g_{\text{EI}}=1.1611,g_{\text{II}}=0.9570\\\ \left\langle V\right\rangle_{\theta_{1}},\dots,\left\langle V\right\rangle_{\theta_{7}}\\\ \nu_{I},\left\langle V\right\rangle_{I}\end{array}$ | 17 | 7 / 16
6 blue | $\begin{array}[t]{c}g_{1}=20Hz\\\ g_{\sigma}=1.2\\\ \\\ \\\ \end{array}$ | $\begin{array}[t]{c}g_{0},g_{r}\\\ w_{+}=2.423,w_{\sigma}=0.4750\\\ g_{\text{EE}}=0.1798,g_{\text{IE}}=0.1858\\\ g_{\text{EI}}=0.7882,g_{\text{II}}=0.7632\\\ \left\langle V\right\rangle_{\theta_{1}},\dots,\left\langle V\right\rangle_{\theta_{7}}\\\ \nu_{I},\left\langle V\right\rangle_{I}\end{array}$ | 17 | 7 / 16
6 green | $\begin{array}[t]{c}g_{1}=25Hz\\\ g_{r}=8\\\ \\\ \\\ \end{array}$ | $\begin{array}[t]{c}g_{0},g_{\sigma}\\\ w_{+}=4.4917,w_{\sigma}=0.0909\\\ g_{\text{EE}}=0.4397,g_{\text{IE}}=0.4742\\\ g_{\text{EI}}=1.1933,g_{\text{II}}=1.1948\\\ \left\langle V\right\rangle_{\theta_{1}},\dots,\left\langle V\right\rangle_{\theta_{7}}\\\ \nu_{I},\left\langle V\right\rangle_{I}\end{array}$ | 17 | 7 / 16

?tablename? 3: Parameters optimized in spiking networks. For constraints in 4
see the parameter values in 1. For network parameters in 6, we additionally
give the values obtained by optimization. _Points_ is the number of sampling
points. _Errors_ is the number errors used for optimization, this includes 2
errors for the inhibitory population, in addition to 2 errors per sampling
point.