Introduction
The predominant view held in neuroscience today is that the activity of neurons in the brain and the function of neural circuits can be understood in terms of the computations that underlie perception, motion control, decision making, and cognitive reasoning. However, recent developments in experimental neuroscience have exposed critical holes in our theoretical understanding of how neural dynamics relates to computation. For example, in the prominent existing theories for neural computation, such as energybased attractor networks (Hopfield, 1984; Amit, 1992), population coding (Pouget et al., 2003) and the neural engineering framework (Eliasmith et al., 2012), information is represented in firing rates of neurons, as opposed to precise timing patterns of actionpotentials. The ratecoding assumption (Bialek et al., 1991; Cessac et al., 2010) is hard to reconcile with dynamic phenomena ubiquitously observed in neural recordings, such as precise sequences of actionpotentials (Abeles, 1982; Reinagel and Reid, 2000; Panzeri and Diamond, 2010; Tiesinga et al., 2008), high speed of neural computation during perception (Thorpe et al., 1996), as well as rhythmic activity of neural populations on various temporal and spatial scales (O’Keefe and Recce, 1993; Buzsaki and Draguhn, 2004).
Here, we develop a theory that can elucidate the role of precise spike timing in neural circuits. The theory is based on complexvalued neural networks and a phasetotiming mapping that translates a complex neural state to the timing of a spike. We first introduce a novel model of fixed point attractor networks in complex state space, called threshold phasor associative memory (TPAM) networks. The model dynamics is governed by a Lyapunov function, akin to Hopfield networks (Hopfield, 1982)
, it has high memory capacity and, unlike Hopfield networks, can store continuousvalued data. Second, we show that any TPAM network possesses a corresponding equivalent network with spiking neurons and timedelayed synapses. Our framework can be used to design circuits of spiking neurons to compute robustly with spike times, potentially trailblazing a path towards fully leveraging recent highperformance neuromorphic computing hardware
(Davies et al., 2018). Concurrently, the framework can help connect computations with experimental observations in neuroscience, such as sequential and rhythmic firing dynamics, balance between excitation and inhibition, and synaptic delays.Background
Fixed point attractor models, and more generally, energybased models
(Ackley et al., 1985), have played an extremely important role in neural networks and theoretical neuroscience. The appeal of these models comes from the fact that their dynamics can be described by the descent of an Energy or Lyapunov function, often conceptualized as an “energy landscape”. While energy descent does not describe all of neural computation, it is the basis of neural network algorithms for many important problems, such as denoising/error correction to make computations robust with respect to perturbations, or constrained optimization (Hopfield and Tank, 1985).The landscape of fixed point attractor neural networks
Here, we focus on a prominent class of fixed point attractor networks with Hebbian oneshot learning to store a set of neural activity patterns. To retrieve a pattern, the network is first initialized with a cue – typically a noisy or partial version of a stored pattern. For retrieval, an iterative dynamics successively reduces the initial noise and converges to a clean version of the stored memory. The models can be distinguished along the following dimensions: realvalued versus a complexvalued neural state space, neurons with discretizing versus continuous transfer functions, and the neural firing patterns in the network being sparse versus dense (Fig. 1).
The traditional models of attractor neural networks, which have a symmetric synaptic matrix that guarantees Lyapunov dynamics and fixed points (Cohen and Grossberg, 1983), were inspired by the Ising model of ferromagnetism (Ising, 1925). The neural transfer function is stepshaped (describing the spiking versus silent state of a neuron), and the models can store dense activity patterns with even ratio between active and silent neurons (Little, 1974; Hopfield, 1982). The inclusion of “thermal” noise yielded networks with a sigmoidshaped neural transfer function (Hopfield, 1984; Treves, 1990; Kühn et al., 1991)
, which produces continuous realvalued state vectors that could be conceptually related to firing rates of spiking neurons in the brain. Further, following the insight of Willshaw et al.
(Willshaw et al., 1969), attractor networks for storing sparse activity patterns were proposed, which better matched biology and exhibited greatly enhanced memory capacity (Tsodyks and Feigel’man, 1988; Buhmann et al., 1989; Palm and Sommer, 1992; Schwenker et al., 1996).One drawback of these traditional memory models is that all attractor states are (close to) binary patterns, where each neuron is either almost silent or fully active near saturation. The restriction to binary memories can be overcome by introducing model neurons that can saturate at multiple (more than two) activation levels
(Gross et al., 1985; Meunier et al., 1989; Yedidia, 1989; Bouten and Engel, 1993). This class of models was inspired by the Potts glass model in solid state physics. Another model with multilevel neurons is the socalled complex Hopfield network (Farhat et al., 1985; Cook, 1989; Gerl et al., 1992; Hirose, 1992; Jankowski et al., 1996; DonqLiang and WenJune, 1998; Aoki and Kosugi, 2000; Aizenberg, 2011; Kobayashi, 2017). Here, the model neurons are discretized phasors, and as a result, the states of the model are complex vectors whose components have unity norm and phase angles chosen from a finite, equidistant set. For a discretization number of two, complex Hopfield networks degenerate to the realvalued bipolar Hopfield network (Hopfield, 1982).A unique strength of a complex state space was highlighted by phasor networks (Noest, 1988; Sirat et al., 1989; Hoppensteadt and Izhikevich, 1996). In phasor networks, the neural transfer function is a phasor projection, i. e. each neural state carries the continuous phase value of the postsynaptic sum, but has normalized magnitude. Interestingly, even without phase discretization, phasor networks can store arbitrary continuousvalued phasor patterns. Patterns with arbitrary relative phase angles can be stable because of the ring topology of normalized complex numbers.
In existing phasor networks and complex Hopfield networks, all neurons represent phases at every time step. To provide models of greater computational versatility, here we explore relatively uncharted territory in Figure 1: Attractor networks with complexvalued sparse activity states. These models can also store phasor patterns, in which a fraction of components have zero amplitude that correspond to silent or inactive neurons. Specifically, we introduce and investigate a novel attractor network model called threshold phasor associative memory (TPAM) network. As will be shown, pattern sparsity in TPAM enables high memory capacity – as in realvalued models (Tsodyks and Feigel’man, 1988), and also corresponds to spike timing patterns that are neurobiologically plausible.
Hebbian sequence associative memories
A first foray into temporal neural coding was the development of networks of threshold neurons with Hebbiantype heteroassociative learning that synaptically stores the firstorder Markovian transitions of sequential neural activity patterns
(Amari, 1972; Willwacher, 1976, 1982; Kleinfeld, 1986; Sompolinsky and Kanter, 1986; Amit, 1988; Riedel et al., 1988). When initialized at or near the first pattern of a stored sequence, the parallelupdate and discrete time dynamics of the network will produce the entire sequence, with a pattern transition occurring each time step. These networks have nonsymmetric synaptic matrices and therefore the dynamics are not governed by a Lyapunov function (Cohen and Grossberg, 1983). However, for networks that store cyclic pattern sequences of equal lengths, Herz et al. (Herz, 1991; Herz et al., 1991) have shown that the network dynamics are governed by an energy function, defined in an extended state space in which each entire sequence is represented by a point.Our approach, to describe temporal structure in the complex domain, is related to the extended state space approach. Specifically, we will show that sequence associative networks correspond to attractor networks whose fixed points are phase patterns with equidistantly binned phase values, each phase bin representing a position in the sequence.
Theories of spiking neural networks
The firing rates in a spiking neural network can be described by attractor networks with sigmoid neural transfer function (Hopfield, 1984; Treves, 1990; Kühn et al., 1991). This modeling approach rests on the paradigm of rate coding (Bialek et al., 1991; Cessac et al., 2010), which considers spiking as an inhomogenious Poisson random process, with synaptic inputs controlling only the instantaneous event rate.
Around 1990, the development of theory to describe networks of neurons with explicit spiking mechanisms began. Using phasereturn maps, Mirollo and Strogatz (Mirollo and Strogatz, 1990) found that integrateandfire neurons coupled without synaptic delays synchronize for almost all initial conditions. Such networks with delayless synaptic couplings formed by outerproduct learning rules have been demonstrated to serve as autoassociative memories. Depending on the operation point, retrieval states can either be represented by synchronious spike patterns (Buhmann and Schulten, 1986; Mueller and Herz, 1999; Sommer and Wennekers, 2001), or by asynchroneous rate patterns (Gerstner and van Hemmen, 1992). Further it was shown that for certain types of delayless synaptic matrices, there are stable limit cycle attractors in which every neuron fires once per cycle (Hopfield and Herz, 1995). The convergence to these dynamic attractor states is rapid and can be described by the descent in a Lyapunov function.
Here, we are particularly interested in networks in which spiking neurons are coupled by synapses, each with individual delay time, whose properties have been studied in small networks (Ernst et al., 1995; Nischwitz and Glünder, 1995). Herz (Herz, 1995) derived a Lyapunov function for networks of nonleaky integrateandfire neurons with synaptic delays. As in the delayless case (Hopfield and Herz, 1995), the descent in the Lyapunov function corresponds to the rapid (although not very robust) convergence to particular temporal coding patterns. Somewhat inspired by (Herz, 1995), our novel approach is to map complex fixed point attractor networks to networks of leaky integrateandfire neurons. As we will show, the resulting spiking neural networks employ spike timing codes that are robust with regard to perturbations, exhibit high memory capacity, and share many properties with neuronal circuitry of real brains.
Results
Threshold phasor associative memories
We propose a novel memory model, the threshold phasor associative memory (TPAM), which can store sparse patterns of complex phasors as fixed point attractors. The network dynamics is governed by an energy function, which we derive. Further, simulation experiments show that TPAM has high memory capacity and provides an efficient error correction stage in a neural network for storing images.
Learning and neural dynamics
Phasor neural networks (Noest, 1988; Sirat et al., 1989) were designed to store dense phase angle patterns in which every component represents a phase angle. Similar to autoassociative outer product Hebbian learning (Kohonen, 1972), like in the standard Hopfield network (Hopfield, 1982), phasor networks employ a complexconjugate outerproduct learning rule:
(1) 
where is a matrix of phasor patterns of dimension . A component of one of the stored patterns is given by . The entries along the diagonal of are set to .
During the discretetime neural update in a phasor network, the output of neuron is a unit magnitude complex number . The complex neural states are multiplied by the complex weight matrix to give the postsynaptic dendritic sum:
(2) 
In contrast to the described phasor memory network, the TPAM network is designed for storing patterns in which only a sparse fraction of components have unit magnitude and the rest have zero amplitude, i.e., are inactive. TPAM uses the same learning rule [1] and postsynaptic summation [2] as the original phasor network, but differs in the neural transfer function. The neural transfer function includes a threshold operation on the amplitude of the synaptic sum [2]:
(3) 
with the Heaviside function. If the threshold is met, the output preserves the phase of the sum vector and normalizes the amplitude. Otherwise, the output is zero.
To maintain a given level of network activation, the threshold setting needs to be controlled as a function of the global network activity (Wennekers et al., 1995). Here, we set the threshold proportional to the overall activity:
(4) 
with a scalar between 0 and 1, typically slightly less than 1.
The memory recall in TPAM with parallel update with neurons is demonstrated in Fig. 2. The network has stored sparse random phasor patterns with
and phase values drawn independently from a uniform distribution. The iterative recall is initialized by a partial memory pattern – with some nonzero components set to zero (top panels), and with a superposition of several stored patterns (bottom panels). In both cases, the network dynamics relaxes to one of the stored memories (approximately).
Energy function of TPAM networks
For traditional phasor memory networks (without threshold), Noest (Noest, 1988) showed that the corresponding Lyapunov function is
(5) 
Note, that because [1] results in a Hermitian matrix , [5] is a realvalued function. Further note, that the dynamics in phasor networks is a generalization of phasecoupled systems well studied in physics, such as the Kuramoto model (Kuramoto, 1975), and the XY model (Kosterlitz and Thouless, 1973), and for describing coherent activity in neural networks (Schuster and Wagner, 1990; Sompolinsky et al., 1990; Niebur et al., 1991). Those models are governed by a Lyapunov function of the form [5], but in which is realvalued and symmetric (Van Hemmen and Wreszinski, 1993).
To see how the inclusion of the threshold operation in the TPAM update [3] changes the Lyapunov function, we follow the treatment in (Hopfield, 1984) by extending [5] to describe the dynamics of phasor networks with arbitrary invertible transfer function :
(6) 
The neural transfer function of TPAM, in [3], is not invertible. But it can be approximated by a smooth, invertible function by replacing the Heaviside function in [3] with an invertible function , for example, the logistic function. In the limit of making the approximation tight, i.e., , the corresponding update is given by [3]. For a constant global threshold the Lyapunov function [6) of TPAM is:
(7) 
with a potential barrier function . According to equation [7], a positive constant global threshold [3] has the effect of adding a constraint term, which encourages a lower activity in the network.
For the dynamic threshold control [4], the Lyapunov function for TPAM becomes:
(8) 
According to equation [8], a positve coefficient in the dynamic threshold control, [3] and [4], adds a repulsive selfinteraction between active phasors, thereby reducing the activity in the network.
The derived Lyapunov functions help to clarify the difference between constant and linear threshold control. Consider the case of low memory load. With constant threshold, not only are the individual stored patterns stable fixed points, but also their superpositions will be stable. In contrast, dynamic threshold control introduces competition between active stored memory patterns. The coefficient can be tuned so that only individual patterns are stable (as done here). When lowered, superpositions of two (or more) patterns can become stable, but competition still only allows a limited number of active superpositions. This may be useful behavior for applications outside the scope of this paper.
Information capacity of TPAM networks
To understand the function of TPAM, the impact of its different features on memory performance are studied through simulation experiments. After storing random patterns, we initialize the network to one of the stored patterns with a small amount of noise. The network runs until convergence or for a maximum of 500 iterations. To assess the quality of memory recall we then compare the network state with the errorless stored pattern.
Figure 3A displays on the axes cosine similarity (i.e. correlation) between the output of the memory and the desired target pattern. This normalized metric accounts for both disparity in the phase offset and mismatch in the supports, but does not directly reflect the mutual information between patterns, which also depends on the sparsity level. Fig. 3A compares TPAM with different levels of sparsity (green) to the traditional binary Hopfield network (Hopfield, 1982) (black line), and to the continuous phasor network (Noest, 1988) (orange line). As in the case of the ternary models (Meunier et al., 1989) (see Supplement), the number of patterns that can be recalled with high precision increases significantly with sparsity.
To assess how the memory efficiency of TPAM depends on pattern sparsity, we empirically measure the information in the random phase patterns that can be recalled from the memory (see Supplement). Dividing the recalled information by the number of synapses yields the memory capacity of a network in bits per synapse (Palm and Sommer, 1992; Schwenker et al., 1996). Measuring the memory capacity in bits, rather than by the number of stored patterns (Hopfield, 1982), has the advantage that performances can be compared for memory patterns with different sparsity levels.
The measurements of memory capacity in bits/synapse shows that sparsity greatly increases memory capacity (Fig. 3B) over dense associative memory networks. Interestingly, this result parallels the increase of memory capacity in binary Hopfield networks with pattern sparsity (Tsodyks and Feigel’man, 1988; Buhmann et al., 1989; Palm and Sommer, 1992). This holds up to a limit, however, as the networks with the highest sparsity levels had a slightly decreased maximum memory capacity in the simulation experiments.
Indexing and retrieving data with TPAM networks
One option to store realworld data in TPAM networks is to encode the data in phasor patterns that can be stored in a recurrent TPAM network as described in the last section. A problem with this approach is that data correlations cause interference in the stored patterns, which is a known issue in traditional associative memories with outer product learning rule that reduces the information capacity quite drastically (Amit, 1992).
Here, we explore the ability of TPAM to perform error correction of random indexing patterns within a memory network inspired by the sparse distributed memory (SDM) model (Kanerva, 1988). The original SDM model consists of two feedforward layers of neurons, an indexing stage with random synaptic weights, mapping data points to sparse binary index vectors, and a heteroassociative memory, mapping index vectors back to data points. Our memory architecture deviates from the original SDM model in three regards. First, it uses complexvalued index patterns. Second, synaptic weights in the indexing stage are learned from the data. Third, it consists of an additional third stage, an errorcorrection stage using a recurrent TPAM, that sits between the indexing stage and heteroassociative memory (similar as in (Aoki and Kosugi, 2000); Fig. 4A).
In the following, we denote the data matrix as , where is the dimensionality of the data and is the number of data points. As in previous sections, we denote the matrix of index vectors (randomly chosen sparse phasor patterns) by where is the number of neurons in the TPAM network.
The indexing stage maps incoming data vectors into index patterns. It is a feedforward network of neurons with the synaptic matrix . A simple Hebbian learning rule for heteroassociation is . However, to reduce the level of indexing errors due to inherent data correlations, we use learning that involves the pseudoinverse of the data,
, resulting from the singular value decomposition
. Specifically, the synapses in the indexing stage are formed according to:(9) 
This linear transform performs
pattern separation by amplifying the differences between correlated data vectors. It thereby produces decorrelated projections from the data space to the index space.The output stage of the network is a heteroassociative memory of the data. This is a feedforward network using simple Hebbian learning for mapping an index pattern back into a data vector. To produce realvalued output patterns, the output neural transfer function projects each component to the real part.
In SDM and heteroassociative memories in general, if the indexing or cue patterns are noisy, the quality of the returned data suffers significantly. To improve the retrieval quality in these cases, we store the indexing patterns in TPAM. The TPAM performs error correction on the index patterns produced by the indexing stage before the heteroassociative memory stage.
Empirical comparisons of image storage and retrieval using a simple Hebbian heteroassociative memory, an SDM, and the full network with pattern separation and TPAM for error correction (Fig. 4B; see Supplement for implementation details) were performed with simulation experiments. The simple Hebbian model and the SDM were also extended by incorporating pattern separation in the indexing stage (solid lines in Fig. 4B include pattern separation). We stored image patches of pixels into the networks, and measured the correlation of the retrieved pattern with the true pattern given a noisy input cue. We compute the total information per pixel as (see Supplement). The full network returns a larger amount of information about the stored data than the simpler models. Errors in retrieved TPAM patterns (Fig. 4C) are due to spurious local minima, which are usually superpositions of stored memories. Similarly, errors in SDM are spurious activations of incorrect patterns, leading to readout errors also as superpositions of stored memories. Including the pseudoinverse for indexing (PINV), improves the likelihood of avoiding such superposition errors.
Relating TPAM networks to spiking neural networks
Here, we exploit a natural link between a complex state space and a spike raster through a phasetotiming mapping. To approximate TPAM networks with networks of integrateandfire neurons, we propose biologically plausible mechanisms for the key computations: complex synaptic multiplication, summation of complex postsynaptic signals [2], and the neural transfer function with dynamic threshold [3].
Phasetotiming mapping
The complex state vector of TPAM (Fig. 5A) can be uniquely mapped to a spike pattern in a population of neurons (Fig. 5B) through a phasetotiming mapping. The phase angle of a component is represented by the timing of a spike within an interval , where the times between and represent the phase angles between and . A stable fixed point of a complex TPAM state corresponds to a limitcycle of precisely timed spiking activity, where neurons fire periodically with a rate of or are silent. The cycle period can be chosen arbitrarily, and is in the following simulation experiments.
Note that the phasetotime mapping is not bijective. To map back from spike rasters to complex vectors, one needs to specify the offset that divides the time axis into intervals of length . Different choices of offsets lead to different sequences of complex vectors. At fixed points, the resulting vector sequences just differ by a global phase shift. Away from fixed points, however, different offsets can also lead to vector sequences with different amplitude structure.
Complex synaptic multiplication
Via the phasetotiming mapping, it is also possible to translate the synaptic phasor interactions of a given TPAM to synapses between spiking neurons. Specifically, a phase shift of a complex TPAM synapse translates into a synapse with a specific conduction delay. The synaptic multiplication then corresponds to the summation of spiketime and synaptic delay. Connecting neurons with a deterministic spiking mechanism (such as integrateandfire) with this synaptic structure results in a network in which neural updates are eventdriven and evolve continuously. To make all interactions causal, a cycle time can be added to synaptic delays that would otherwise be negative.
In the TPAM with discrete time dynamics, the new state of a neuron is exclusively a function of the previous time state. Thus, for a given phasetotiming mapping, spike responses in one time window should only depend on spikes within the immediately preceding time window. For any fixed point, multiples of can be added or subtracted to the synaptic delays to match the discrete time dynamics (Fig. 5C; dashed). However, this multiple of depends on the particular fixed point, and the mapping to discrete time dynamics cannot be simultaneously guaranteed for multiple fixed points stored in the network.
Here, we selectively increment delays by in order to wrap the synaptic delay values to a range between and . This choice maximizes interactions between subsequent time steps, but it cannot fully eliminate interactions within the same timestep or across two timesteps (Fig. 5C, D; solid). Yet, at any fixed point this choice of synaptic delays produces exactly the same postsynaptic potentials as the synaptic delays that correspond to discrete time dynamics. Thus, even though the spiking network does not strictly obey a discrete time dynamics, at each equilibrium state an equivalent network exists that satisfies discrete time dynamics and produces exactly the same postsynaptic inputs to all neurons.
Postsynaptic complex vector summation
The complex sum of TPAM [2] corresponds to the addition of sine waves with different phases and amplitudes that oscillate with period in the time domain. To perform this computation with a spiking neural network, the effect of each spike should produce a phaselocked oscillatory zeromean postsynaptic current. In our model of standard integrateandfire neurons (Goodman and Brette, 2009; Eliasmith and Anderson, 2003; Izhikevich, 2007), the postsynaptic oscillation caused by each spike is generated by a monosynaptic excitatory current, as well as inhibitory current routed through interneurons. The actionpotential and the opening of synaptic channels act on a fast timescale, which allows these kinetic features to be simplified as instant jumps in the dynamic state variables. When a neuron’s membrane potential reaches threshold , it fires a spike and the membrane potential is set to the reset potential . After the synaptic delay, the spike causes postsynaptic channels to open. This is modeled as a jump in a synaptic state variable that injects current proportional to the synaptic magnitude , which then decays exponentially as the channels close (see Supplement for details).
The time constants determining the decays of neural and synaptic variables are tuned proportional to the cycle time to create the oscillatory current. The time constant of the membranes of excitatory neurons is , and for interneurons. The time constant for inhibitory and excitatory synapses is and , respectively. With these settings, the total postsynaptic current elicited by a spike forms (approximately) a sine wave oscillation with cycle period of (Fig. 6A).
The excitatory synapses have individual synaptic delays determined by the phases of the recurrent weight matrix , whereas the delays to inhibitory neurons are all the same. Because inhibitory neurons operate linearly and their synaptic delays are nonspecific, a single inhibitory neuron (or arbitrary sized population) can serve to route inhibition for multiple complex synaptic connections. The delays of the excitatory synapses determine the phase of the oscillatory currents. Depending on presynaptic spike times and synaptic delays, the currents can sum either decoherently or coherently (Fig. 6
B, C), estimating the complex sum. Altogether, the complex dot product is approximated by the spiking network.
Neural transfer function and threshold control
The inhibitory neurons serve a second purpose besides shaping the postsynaptic oscillatory currents. They also account for the global normalization needed for the threshold strategy [3], which keeps the activity sparse. The inhibitory population integrates the activity from the excitatory neurons and routes inhibition back into the population in proportion to the overall activity. The magnitude of this feedback inhibition can be tuned to limit the activity of the population, and it implements the threshold strategy needed for TPAM, e.g. [4]. The gain of the inhibitory feedback is modulated by several mechanisms: the gain of the EtoI synapses, the gain of the ItoE synapses, the time constants of the synapses, the number of inhibitory neurons, and the gain of the inhibitory neural transfer function. Each of these gain factors can be analytically accounted, with a linear approximation being useful to understand the gain of the spiking neural transfer function (see Supplement).
The deterministic dynamics of the integrateandfire neuron will cause the neuron to spike whenever the current drives the voltage above threshold. If the magnitude of the input current oscillation is large enough, then the neuron will fire at a consistent phase. For the excitatory neurons, a refractory period slightly over half the cycle time (i.e. ) acts as the Heaviside function on the magnitude. This implements the phasorprojection of TPAM by limiting the neuron to one spike per cycle, while preserving the phase through the spike timing [3]. The parameter value sets the threshold of whether the neuron will fire or not [4].
Simulation experiments with spiking networks
For storing a small collection of RGB images, a network architecture with spiking neurons was implemented as depicted in Fig. 4A. The indexing stage transforms the realvalued input image into a complex vector from the encoding matrix (see Supplement). The complex vector is mapped into a timed sequence of input spikes (Fig. 7A), which is the initial input to a spiking TPAM network. The spiking TPAM network is decoded with a Hebbian heteroassociative memory. This output stage uses the same synaptic mechanisms to implement the complex dot product as described for TPAM. However, the readout neurons respond proportionally to the input magnitude (i.e. no refractory period, see Supplement).
In the simulation experiment, the network is cued with several overlapping patterns and noise. After initialization through the indexing stage, the spiking activity in the TPAM network quickly settles into an oscillatory pattern (Fig. 7B), which corresponds to the index of one of the stored patterns (Fig. 7C). Interestingly, the global oscillation of the network is generated intrinsically – it does not require any external driving mechanism. The output of the heteroassocitive memory stage during the convergence of the TPAM dynamics (Fig. 7D) shows how the dynamics rapidly settles to one of the stored patterns superposed in the input, outcompeting the other two.
Examples of the dynamics of individual neurons in the network are illustrated. A neuron participating in the stable oscillation receives coherent volleys of excitatory input and has a large oscillation in its membrane potential (Fig. 7E). Another neuron shown is initially excited, but then becomes silent in the stable attractor state (Fig. 7F). In this neuron excitation becomes decoherent, and as a consequence the voltage is never large enough to overcome threshold.
Finally, the capacity of the spiking network is examined in simulation experiments (Fig. 7G). The spiking network is robust without much parameter tuning and not perfectly optimized, but still retrieval performance based on the number of stored patterns (Fig. 7G, black x’s) easily exceeds the performance of a traditional bipolar Hopfield network (Fig. 7G, black line). The performance curve of the spiking model is nicely described by the performance curve of the complex TPAM with similar settings. However, the spiking network does not reach the full performance of the complex TPAM with optimized threshold. The sparsity and threshold parameters of the spiking network were not fully optimized in our experiments. But also, some of the approximations in the spiking network add noise, which prevents it from reaching the full capacity of the ideal complex model. Nonetheless, these experiments show that the spiking model does behave in a manner that can be captured by the complex TPAM model.
Sequence associative memories and complex attractor networks
Last, we investigate how complex fixed point attractor networks can help to understand sequence associative memories: simple networks with binary threshold neurons and parallel, timediscrete update dynamics for storing sequences of patterns of fixed length (see Background). Consider the storage of a closed sequences or limit cycles of fixed length : , with . In case of storing multiple sequences an index is added to label the different sequences: . The learning in these models is also described by a Hebbian outerproduct learning scheme (Amari, 1972)
. Here we use a combination of Hebbian and antiHebbian learning to produce a skewsymmetric interaction matrix:
(10) 
Since the matrix is skewsymmetric, there is no Lyapunov function describing the limit cycle dynamics in the network. However, we can use the spectral properties of the weights to construct an equivalent fixed point attractor network.
Consider [10] for the simple example with , and the stored patterns being the cardinal basis vectors of
. One of the complex eigenvectors of
is , which is the (equidistant) phasor pattern that represents the entire stored limit cycle in complex space. One can now form a complex matrix that posessesas a fixed point, i.e., has eigenvalue of one, simply by dividing
by the eigenvalue associated with , which is :(11) 
Since the eigenvalues of any skewsymmetric matrix have zero real part (Eves, 1980) the interaction matrix is always Hermitian in general. Thus, the described construction is a recipe to translate sequence memory networks into complex neural networks governed by a Lyapunov dynamics. In the resulting networks, the synaptic matrix is , the neural nonlinearity is , and the Lyapunov function is [5].
One could now suspect that storing the pattern in a phasor network (Noest, 1988) would result in the same network interaction matrix . However, this is not the case. The weight matrix resulting from learning the phase vector with the conjugate outerproduct learning rule [1] is:
(12) 
The phase vector is again an eigenvector of .
The take away from this example are the following points:
(1.) Outerproduct sequence memories with skewsymmetric weights can be mapped to continuousvalued phasor networks with Hermitian weights , whose dynamics is described by the Lyapunov function [5]. A similar idea of deriving a Lyapunov function for periodic sequences of binary patterns was proposed in (Herz, 1991; Herz et al., 1991). Specifically, these authors proposed to embed sequence trajectories in a realvalued space, consisting of copies of the original state space.
(2.) Continuousvalued phasor networks derived from sequence memories with outerproduct rule [10] are different from phasor networks using the conjugate complex outerproduct learning rule [1], such as TPAM networks.
(3.) Phasor networks derived from outerproduct sequence memories have two severe restrictions. First, since the synaptic weights are imaginary without real part, they only can rotate the presynaptic signals by 90 degrees. Second, the firstorder Markov property of sequence memory networks translates into phasor networks with interactions only between direct phase angle neighbors.
(4.) Phasor networks derived from sequence memories can, by construction, only store phase patterns whose phase angles are equidistantly dividing . Because of symmetry reasons, such equidistant phase patterns can be stabilized in spite of the restrictions described in the previous point.
Discussion
We present a theory framework for temporal coding with spiking neurons that is based on fixed point attractor dynamics in complex state space. Our results describe a new type of complex attractor neural network, and how it can be used to design robust computations in networks of spiking neurons with delay synapses.
Threshold phasor networks
A novel type of complex attractor network is introduced, threshold phasor associative memory (TPAM). TPAM inherits from previous complex fixed point attractor networks the Lyapunov dynamics, and the capability to store arbitrary continuousvalued phasor patterns. The neural threshold mechanism added in TPAM yields advantages over existing phasor memory networks (Noest, 1988). First, it significantly increases the memory capacity, the amount of information that can be stored per synapse, as we show in simulation experiments. Second, the amplitude threshold operation in complex space prevents neurons to represent phase angles whose postsynaptic phase concentration is diffuse. Therefore, all phase angles in a TPAM state vector carry high information about the postsynaptic inputs. Third, TPAM networks describe physiologically realistic sparse firing patterns.
The ability to store continuous phasor patterns as attractor states permits the processing of analog data in TPAM without preceding data discretization. However, data correlations still causes problems in TPAM networks as in other associative memory models. For the storage of correlated sensor data, such as images, we propose a threelayer network extending from previous memory indexing models, such as sparse distributed memory (Kanerva, 1988)
. We demonstrate in an image retrieval task that a network consisting of a pattern separation stage, a TPAM network for error correction, and a heteroassociative memory has improved performance, compared to earlier models.
Mapping complex attractors into periodic temporal codes
Previous models of spiketiming based computation often proved brittle, which even served as an argument for rendering spiketiming codes as irrelevant for neuroscience (London et al., 2010). We employ a straightforward phasetotiming mapping to translate the complex fixed point attractors in TPAM to periodic patterns of precisely timed spiking activity. Although the timecontinuous nature of spike interactions through delay synapses and the timediscrete dynamics in TPAM are not equivalent in the entire state space, they are so at the fixed points. The correspondence between complex fixed point states and periodic spiking patterns enables robust computation with spiketiming codes, as we demonstrate with integrateandfire neural networks. Such periodic spiking patterns are (superficially) similar to and compatible with spiketiming based codes proposed previously. For example, Hopfield (Hopfield, 1995) proposed a network for transducing ratecoded signals into periodic spiketiming patterns, which encodes the realvalued input into the phase of a complex number. Other examples are synfire braids (Bienenstock, 1995) and polychronization (Izhikevich, 2006) in networks with synaptic delays and neurons that detect synchrony.
While TPAM provides valuable insights into network computations using the timing of single spikes, a simple extension is able to also capture aspects of rate coding. By relaxing the constraint on binary magnitudes in TPAM, complex values with variable magnitudes can be represented, which correspond to variable firing rates. For instance, if the neural response magnitude is a saturating sigmoid function, attractors would lie near saturation points, as in the realvalued Hopfield model
(Hopfield, 1984). Yet, phase values can still be arbitrary. The corresponding spike patterns would still be oscillatory, with bursts at a particular phase followed by durations of silence in the antiphase. This type of burst spiking is indeed seen in hippocampal place cells that are locked to the theta rhythm (O’Keefe and Recce, 1993).Further, it is shown how Hebbian sequence memories in discrete time (Amari, 1972) can be described by fixed point dynamics in the complex domain. This description is similar to the construction of a Lyapunov function for sequence memories in an enhanced (realvalued) state space (Herz et al., 1991). Interestingly, the complex attractor networks corresponding to sequence memories [11] are different from phasor networks with conjugate outerproduct Hebbian learning [1].
Complex algebra with spiking neural circuits
We describe how the postsynaptic computation in TPAM neurons can be implemented by circuits of excitatory and inhibitory leaky integrateandfire neurons. In essence, the required complex dot product can be achieved by shaping the postsynaptic currents caused by individual spikes into (approximately) a sinusoid. The shaping is achieved through a combination of dendritic filtering, synaptic delays and adding inhibitory neurons that balance excitation. The required mechanisms are in accordance with the following observations in neuroscience: Periodic firing and synchrony: Actionpotentials in the brain are often periodic, synchronized with intrinsic rhythms visible in local field potentials (Buzsaki and Draguhn, 2004; Riehle et al., 1997; Lynch et al., 2016). Delayed synaptic transmission: Due to variability in axon length and myelination, the distribution in measured delay times in monosynaptic transmission is broad, easily spanning the full cycle length of gamma and even theta oscillations (Swadlow, 2000). Balance between excitation and inhibition: Excitation/inhibition balance is widely seen throughout cortex (Mariño et al., 2005; Haider et al., 2006; Atallah and Scanziani, 2009), and inhibitory feedback onto pyramidal cells is a major feature of the canonical cortical microcircuit (Douglas et al., 1989; Isaacson and Scanziani, 2011).
Complex algebra with spiking neurons is relevant for other types of complex neural networks as well, such as working memory networks based on the principles of reservoir computing (Frady et al., 2018). Because spiking neurons represent both the real and imaginary parts, complex spiking reservoir networks would have twice the memory capacity per neuron than ratebased models. Further, TPAMlike networks with modified learning rules could potentially be used to construct spiketiming implementations of lineattractors–useful for understanding place coding in hippocampus (Welinder et al., 2008; Campbell et al., 2018), or for modeling/predicting dynamical systems (Eliasmith et al., 2012; Denève et al., 2017).
The conjugate outerproduct learning rule in TPAM requires tunable synaptic delays, as suggested in previous models (Hüning et al., 1998), but without compelling biological evidence. Alternatively, the wide distribution of fixed synaptic delays in brain tissue could enable biologically observed spiketimingdependent plasticity (STDP) (Szatmáry and Izhikevich, 2010) to shape synaptic connectivity compatible with the conjugate outerproduct learning rule. An interesting subject of future research is the effect of STDP mechanisms in the presence of oscillations and synaptic delays, and how it relates to outerproduct Hebbian learning rules.
The presented theory also has impact on neuromorphic computing. Several groups have been exploring and engineering spiking neural networks as a new paradigm for computation, such as Braindrop (Neckar et al., 2019) and IBM’s True North (Merolla et al., 2014). Recently, Intel announced the neuromorphic chip Loihi (Davies et al., 2018), with features such as individual synaptic delays and onchip learning. Our theory offers a principled way of “programming” spikingneuron hardware, leveraging the speed of temporal codes and providing straightforward connections to complex matrix algebra and Lyapunov dynamics.
Acknowledgements
This work was supported by the Intel Corporation (ISRA on Neuromorphic Architectures for Mainstream Computing), NSF award 1718991, and Berkeley DeepDrive. The authors would like to thank Pentti Kanerva and members of the Redwood Center for valuable feedback.
Supplemental Methods
Threshold phasor associative memories
Derivation of the Lyapunov function for TPAM
To derive a Lyapunov function for TPAM, we first consider Hopfield networks with realvalued neurons. Hopfield (Hopfield, 1984) showed that for an invertible neural transfer function , the corresponding Lyapunov function is:
(13) 
The first term is the energy function of the Hopfield network with binary neurons (Hopfield, 1982). The derivative of the second term in [13] is just . Thus, setting the derivative with respect to to zero leads to coordinatewise update equations with transfer function : .
To understand the behavior of binary neurons with nonzero thresholds, we consider the neural transfer function , where is the constant, individual threshold of neuron . Although the transfer function is not invertible, it can be approximated by an invertible function, for example, the logistic function. In the limit of making the approximation tight, i.e., , the Lyapunov function [13] becomes
(14) 
Thus, the second term in [13] decomposes into a barrier and a bias term. The potential barrier term is , which blocks components to assume any values outside the interval . The other is a bias term, well known from the Ising model. Note that both additional terms in [14] vanish for governing the updates of binary threshold neurons with , leading to the energy function of the orginial Hopfield model (Hopfield, 1982). The barrier term is unnecessary because enforced implicitly by any binaryvalued neural update.
If all neurons in the network experience the same global threshold, , and the barrier term enforces all state components to lie within the interval, the bias term becomes an constraint:
(15) 
For the term penalizes states with nonzero components, encouraging sparser states with fewer active neurons.
Further, if the threshold control is dynamic, with the threshold a linear function of the network activity, , the Lyapunov function becomes:
(16) 
Thus, a linear dynamic threshold control corresponds to a antiferromagnetic term added to the network interactions, which can be easily modeled by inhibition.
The Lyapunov functions derived above can be easily be generalized to TPAM. In analogy to [13], the energy function of the phasor neural network (Noest, 1988) for arbitrary invertible transfer function can be extended to:
(17) 
The neural transfer function of TPAM, can be approximated by an invertible function by replacing the Heaviside function in the neural transfer function with an invertible function, such as the logistic function. In the limit of making the approximation tight, i.e., , and with constant threshold the Lyapunov function [6] becomes
(18) 
the original energy function of phasor neural networks (Noest, 1988) with, again, two additional terms. The barrier function in this case is , the bias term is again a term, which encourages networks with lower activity.
In analogy with [16], for linear dynamic threshold control in the threshold strategy, the Lyapunov function of TPAM becomes:
(19) 
Thus, the transition from a constant threshold setting to a dynamic threshold control replaces the constraint by a self interaction term.
Comparison of TPAM to other associative memory models
We performed capacity experiments for previous models of associative memory to compare with TPAM. In 8A, we compare a phasor memory network with continuous phase variables (without threshold) to complex Hopfield networks with discretized phase representations, in which the full phase circle is equidistantly divided into bins (Aoki and Kosugi, 2000). For , this becomes the traditional bipolar Hopfield network. For all models the similarity is high at small , then falls off rapidly. In models with larger numbers of bins, the dropoff of the similarity starts at smaller numbers of stored patterns.
Fig. 8B compares threshold phasor networks with binary phase discretization. The black line is the Hopfield model without threshold storing dense bipolar patterns (same as the black line in subfigure A). The blue lines correspond to models with threshold storing sparse ternary patterns containing components (the lighter the color, the sparser the patterns). For ternary attractor networks, the capacity decreases with moderate sparsity values (dark blue lines below the black line). As the sparsity increases further, the performance supersedes the bipolar Hopfield network (light blue lines). Thus, at sufficient levels of sparsity, the thresholded models can store significantly more patterns than the standard Hopfield model.
Information theory for TPAM capacity
To measure the total information in sparse phasor vectors, we treat the information in phases and amplitudes separately.
The amplitude structure is binary and one has just to estimate the two error types, false positives and misses in the simulation experiments. The information needed to correct the errors of a sparse binary vector is given by:
(20) 
where is the Shannon information, and .
The phase information is estimated from the statistics of phasor vectors observed in the simulation experiments. A von Mises distribution is fit to the difference between the true phasor and retrieved phasor variable, yielding the phase concentration parameter
which is inversely proportional to the variance. In the highfidelity regime, when the network has only a few patterns stored, the von Mises concentration parameter
can be approximated with a Gaussian fit. We measured empirically from simulation experiments. The entropy of the von Mises distribution based on is:(21) 
and the information of a phasor variable is:
(22) 
The total information for a single item stored in memory is then the information for each phasor variable plus the information of the sparse binary vector:
(23) 
Given that the vectors are i.i.d., each item stored yields the same information. The total information is then given by the information per item and the total number of items stored, normalized by the number of synapses:
(24) 
Heteroassociative capacity experiments
We compared three models of heteroassociative memories: a simple model with Hebbian learning, the Sparse Distributed Memory, and a novel indexing method with TPAM as an autoassociative cleanup memory. We compared designs that have random indexing to designs that also include pattern separation in the indexing stage (dashed lines versus solid lines, respectively, in Fig. 4B).
For random indexing, the encoding matrix is random. For the Hebbian learning model and SDM, the encoding matrix was a random sparse binary matrix. For TPAM the encoding matrix was sparse random phasors with . The retrieval of a data vector with a given phasor pattern is , with the decoding matrix.
The decoding matrix for each method is computed from the activity of the hidden layer. For the Hebbian learning model, , where is the cue pattern. The decoding matrix is then . The SDM has a nonlinear threshold function in the hidden layer that can provide some cleanup, . In this case, the function takes the largest values of the input vector and sets them to 1, with all others set to 0. The decoding matrix is similarly .
For the TPAM with random indexing, a codebook of random phasors is chosen . The encoding matrix is . The index patterns are stored in the autoassociative recurrent matrix . The index vector can be decoded with the simple heteroassociative conjugate outerproduct learning rule
To include pattern separation, the indexing matrices are altered. Rather than being purely random, they include both a random part (the random index pattern) and the pseudoinverse of the patterns, . We also replaced the original SDM encoding matrix with a pattern separation matrix to understand the effects of pattern separation and cleanup on memory retrieval performance. For Hebbian learning and SDM models, the matrix is a random sparse binary matrix.
The index patterns are random, but in highdimensions this means that they are approximately orthogonal. We thus do not need to worry too much about correlations in the index vectors (however, performance can be increased by ensuring that the index patterns are exactly orthogonal).
If the cue pattern is a perfect index pattern, i.e. , then the postsynaptic output can be written as:
(25) 
where is the th cardinal basis vector of . The RHS of [25] contains the signal and the noise term. The noise term is zero if the index vectors are exactly orthogonal. Random index patterns will have Gaussian interference noise, and readout will have a signaltonoise ratio of approximately (Frady et al., 2018). The crosstalk noise is complexvalued and will prevent [25] to produce a purely realvalued pattern.
Relating TPAM networks to spiking neural networks
Details of spiking model
Simulations were carried out using the Brian 2 simulator for spiking neural networks (Goodman and Brette, 2009). The dynamics of the integrateandfire neurons is given by:
(26) 
where is the capacitance, is the leak conductance and is the resting potential. When the voltage of the neuron exceeds threshold then the neuron fires a spike and is reset to the reset potential .
Each presynaptic spike causes a synaptic state variable to increment. The synaptic state variable instantly jumps up after the delay time , and decays exponentially:
(27) 
where indicates the timing of the presynaptic spikes.
Each synapse contributes current to the postsynaptic cell proportional to its statevariable :
(28) 
The excitatory currents are balanced by inhibitory currents, which are routed through a population of inhibitory neurons. The timeconstants of the synapses are set to for inhibitory and for excitatory, and the current recombines to create an oscillation with cycle period of . The filtering of the neural capacitance also contributes to shaping the synaptic inputs into a sine wave. The neural timeconstant is for excitatory neurons and for inhibitory neurons. These are relatively small and their contributions to the overall dynamics can be generally ignored. The timeconstants can be tuned to better approximate a perfect oscillation.
The gain of the inhibition can be controlled through several mechanisms. Scaling the synaptic weights, which can be the EtoI or ItoE synapses, the size of the inhibitory population, and the timeconstants of the synaptic connection each contribute a proportional factor to the overall gain of the inhibition. The final factor is accounted by the gain of the neural transfer function, which can be approximated as follows.
The integrateandfire neuron has an analytically defined approximation of the neural transfer function (for fixed/slow input), which we utilize to set the parameters of the spiking model. Based on a constant current into the neuron, one can compute the time it takes to integrate from the reset potential to the threshold potential:
(29) 
where , , and are parameters of the integrateandfire neuron model. The ‘current’ is the input value, directly related to the dendritic sum variable in the normative TPAM network.
The instantaneous firing rate () is used to map the neural transfer function to the spiking neurons. It is the inverse of the spiking period: . If there is no refractory period (), then the will asymptotically converge to a straight line for large input currents , with
(30) 
With a refractory period, the will have similar properties, but will saturate at .
The refractory period is used to change the nonlinear transfer function of the neuron. For neurons as part of a TPAM network, with a phasorprojection as the nonlinearity, a large refractory period that limits the neurons to one spike per cycle can mimic the phasorprojection. Typically, we use something like .
For the inhibitory population and the readout neurons, which act linearly with input current, we set the parameters so that the is a linear function of the input, by using zero (or very small) refractory period. The gain of the neural transfer function is then estimated by the slope of the IFR [30].
The parameters of the network are chosen to approximate biological parameters. The capacitance of the neurons does affect the dynamics of the network, and can itself act as a slight fixed delay. However, any such delay due to timeconstant could be compensated for by reducing the synaptic conduction delays. The jump discontinuity in the synaptic dynamics could be modeled with a risetime, which can be better tuned to produce a more ideal postsynaptic currents. With careful consideration, the more detailed parameters of the spiking model could be better accounted for. However, the parameters do not have to be perfect to get a model working.
Spiking capacity experiments
To be sure that the spiking model approximations did not catastrophically interfere with the attractor dynamics, we ran similar capacity experiments on the spiking model as the normative model, but at a much smaller scale. We chose a fixed parameter set of neurons with active and stored random phasor patterns in the network. These parameters were not the optimized parameters.
The spiking network was initialized to one of the stored patterns and iterated for 5 seconds, which is about 25 cycles. We measured the spike timings in the last second of the simulation and computed the exact frequency of the spiking oscillation. This was then translated into a complex vector. The complex phase of each element in this vector was then rotated to best align with the target pattern. The similarity is then the normalized dot product between the aligned vector and the target pattern.
The details of the spiking model currently prevent an exact parameter match between the normative and the spiking model, but we do see that a consistent spiking model follows the same capacity trajectory as a normative model. Without much parameter optimization, we can build a spiking attractor network that can easily store more patterns than a traditional Hopfield model.
References
 Abeles (1982) Abeles, M. (1982). Local Cortical Circuits: An Electrophysiological Study. Springer, Berlin.

Ackley et al. (1985)
Ackley, D. H., Hinton, G. E., and Sejnowski, T. J. (1985).
A learning algorithm for boltzmann machines.
Cognitive Science, 9(1):147–169.  Aizenberg (2011) Aizenberg, I. (2011). Why we need complexvalued neural networks? In ComplexValued Neural Networks with MultiValued Neurons, pages 1–53. Springer.
 Amari (1972) Amari, S.I. (1972). Learning patterns and pattern sequences by selforganizing nets of threshold elements. IEEE Transactions on Computers, C21(11):1197–1206.
 Amit (1988) Amit, D. J. (1988). Neural networks counting chimes. PNAS, 85(7):2141–2145.
 Amit (1992) Amit, D. J. (1992). Modeling brain function: The world of attractor neural networks. Cambridge university press.
 Aoki and Kosugi (2000) Aoki, H. and Kosugi, Y. (2000). An Image Storage System Using ComplexValued Associative Memories. IEEE, pages 626–629.
 Atallah and Scanziani (2009) Atallah, B. V. and Scanziani, M. (2009). Instantaneous Modulation of Gamma Oscillation Frequency by Balancing Excitation with Inhibition. Neuron, 62(4):566–577.
 Bialek et al. (1991) Bialek, W., Rieke, F., de Ruyter van Steveninck, R. R., and Warland, D. (1991). Reading a Neural Code. Science, 252.
 Bienenstock (1995) Bienenstock, E. (1995). A model of neocortex. Network: Computation in Neural Systems, 6(2):179–224.
 Bouten and Engel (1993) Bouten, M. and Engel, A. (1993). Basin of attraction in networks of multistate neurons. Physical Review E, 47(2):1397.
 Buhmann et al. (1989) Buhmann, J., Divko, R., and Schulten, K. (1989). Associative memory with high information content. Physical Review A, 39(5):2689.
 Buhmann and Schulten (1986) Buhmann, J. and Schulten, K. (1986). Associative recognition and storage in a model network of physiological neurons. Biological cybernetics, 54(45):319–335.
 Buzsaki and Draguhn (2004) Buzsaki, G. and Draguhn, A. (2004). Neuronal Oscillations in Cortical Networks. Science, 304(5679):1926–1929.
 Campbell et al. (2018) Campbell, M. G., Ocko, S. A., Mallory, C. S., Low, I. I., Ganguli, S., and Giocomo, L. M. (2018). Principles governing the integration of landmark and selfmotion cues in entorhinal cortical codes for navigation. Nature Neuroscience, 21(8):1096–1106.
 Cessac et al. (2010) Cessac, B., PaugamMoisy, H., and Viéville, T. (2010). Overview of facts and issues about neural coding by spikes. Journal of Physiology Paris, 104(12):5–18.
 Cohen and Grossberg (1983) Cohen, M. A. and Grossberg, S. (1983). Absolute stability of global pattern formation and parallel memory storage by competitive neural networks. IEEE transactions on systems, man, and cybernetics, (5):815–826.
 Cook (1989) Cook, J. (1989). The meanfield theory of a qstate neural network model. Journal of Physics A: Mathematical and General, 22(12):2057.
 Davies et al. (2018) Davies, M., Srinivasa, N., Lin, T. H., Chinya, G., Cao, Y., Choday, S. H., Dimou, G., Joshi, P., Imam, N., Jain, S., Liao, Y., Lin, C. K., Lines, A., Liu, R., Mathaikutty, D., McCoy, S., Paul, A., Tse, J., Venkataramanan, G., Weng, Y. H., Wild, A., Yang, Y., and Wang, H. (2018). Loihi: A Neuromorphic Manycore Processor with OnChip Learning. IEEE Micro, 38(1):82–99.
 Denève et al. (2017) Denève, S., Alemi, A., and Bourdoukan, R. (2017). The Brain as an Efficient and Robust Adaptive Learner. Neuron, 94(5):969–977.
 DonqLiang and WenJune (1998) DonqLiang, L. and WenJune, W. (1998). A multivalued bidirectional associative memory operating on a complex domain. Neural Networks, 11(9):1623–1635.
 Douglas et al. (1989) Douglas, R. J., Martin, K. A., and Whitteridge, D. (1989). A Canonical Microcircuit for Neocortex. Neural Computation, 1(4):480–488.
 Eliasmith and Anderson (2003) Eliasmith, C. and Anderson, C. H. (2003). Neural Engineering: Computation, Representation, and Dynamics in Neurobiological Systems. MIT Press, Cambridge, MA.
 Eliasmith et al. (2012) Eliasmith, C., Stewart, T. C., Choo, X., Bekolay, T., DeWolf, T., Tang, Y., and Rasmussen, D. (2012). A LargeScale Model of the Functioning Brain. Science, 338(6111):1202–1205.
 Ernst et al. (1995) Ernst, U., Pawelzik, K., and Geisel, T. (1995). Synchronization induced by temporal delays in pulsecoupled oscillators. Physical review letters, 74(9):1570.
 Eves (1980) Eves, H. W. (1980). Elementary matrix theory. Courier Corporation.
 Farhat et al. (1985) Farhat, N. H., Psaltis, D., Prata, A., and Paek, E. (1985). Optical implementation of the hopfield model. Applied optics, 24(10):1469–1475.

Frady et al. (2018)
Frady, E. P., Kleyko, D., and Sommer, F. T. (2018).
A theory of sequence indexing and working memory in recurrent neural networks.
Neural Computation, 30(6):1449–1513.  Gerl et al. (1992) Gerl, F., Bauer, K., and Krey, U. (1992). Learning withqstate clock neurons. Zeitschrift für Physik B Condensed Matter, 88(3):339–347.
 Gerstner and van Hemmen (1992) Gerstner, W. and van Hemmen, J. L. (1992). Associative memory in a network of ‘spiking’neurons. Network: Computation in Neural Systems, 3(2):139–164.
 Goodman and Brette (2009) Goodman, D. F. M. and Brette, R. (2009). The brian simulator. Frontiers in neuroscience, 3(2):192–7.
 Gross et al. (1985) Gross, D., Kanter, I., and Sompolinsky, H. (1985). Meanfield theory of the potts glass. Physical review letters, 55(3):304.
 Haider et al. (2006) Haider, B., Duque, A., Hasenstaub, A., and McCormick, D. A. (2006). Neocortical Network Activity In Vivo Is Generated through a Dynamic Balance of Excitation and Inhibition. Journal of Neuroscience, 26(17):4535–4545.
 Herz (1991) Herz, A. V. (1991). Global analysis of parallel analog networks with retarded feedback. Physical review A, 44(2):1415.
 Herz et al. (1991) Herz, A. V., Li, Z., and Van Hemmen, J. (1991). Statistical mechanics of temporal association in neural networks with transmission delays. Physical review letters, 66(10):1370.
 Herz (1995) Herz, A. V. M. (1995). Do signal delays change rapid phase locking of pulsecoupled oscillators. Preprint.
 Hirose (1992) Hirose, A. (1992). Dynamics of Fully ComplexValued Neural Networks. Electronics Letters, 28(16):1492–1494.
 Hopfield (1982) Hopfield, J. J. (1982). Neural networks and physical systems with emergent collective computational abilities. PNAS, 79(8):2554–2558.
 Hopfield (1984) Hopfield, J. J. (1984). Neurons with graded response have collective computational properties like those of twostate neurons. PNAS, 81(10):3088–3092.
 Hopfield (1995) Hopfield, J. J. (1995). Pattern recognition computation using action potential timing for stimulus representation. Nature, 376(6535):33–36.
 Hopfield and Herz (1995) Hopfield, J. J. and Herz, A. V. (1995). Rapid local synchronization of action potentials: Toward computation with coupled integrateandfire neurons. PNAS, 92(15):6655–6662.
 Hopfield and Tank (1985) Hopfield, J. J. and Tank, D. W. (1985). "Neural" computation of decisions in optimization problems. Biological Cybernetics, 52(3):141–152.
 Hoppensteadt and Izhikevich (1996) Hoppensteadt, F. C. and Izhikevich, E. M. (1996). Synaptic organizations and dynamical properties of weakly connected neural oscillators: I. Analysis of a canonical model. Biological Cybernetics, 75(2):117–127.
 Hüning et al. (1998) Hüning, H., Glünder, H., and Palm, G. (1998). Synaptic delay learning in pulsecoupled neurons. Neural computation, 10(3):555–565.
 Isaacson and Scanziani (2011) Isaacson, J. S. and Scanziani, M. (2011). How inhibition shapes cortical activity. Neuron, 72(2):231–243.
 Ising (1925) Ising, E. (1925). Beitrag zur theorie des ferromagnetismus. Zeitschrift für Physik, 31(1):253–258.
 Izhikevich (2006) Izhikevich, E. M. (2006). Polychronization: computation with spikes. Neural computation, 18(2):245–82.
 Izhikevich (2007) Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. The MIT Press, Cambridge, MA.
 Jankowski et al. (1996) Jankowski, S., Lozowski, A., and Zurada, J. M. (1996). Complexvalued multistate neural associative memory. IEEE Transactions on Neural Networks, 7(6):1491–1496.
 Kanerva (1988) Kanerva, P. (1988). Sparse distributed memory. MIT press.
 Kleinfeld (1986) Kleinfeld, D. (1986). Sequential state generation by model neural networks. Biophysics, 83(December):9469–9473.
 Kobayashi (2017) Kobayashi, M. (2017). Fast Recall for ComplexValued Hopfield Neural Networks with Projection Rules. Comoputational Intelligence and Neuroscience, 2017:1–6.
 Kohonen (1972) Kohonen, T. (1972). Correlation matrix memories. IEEE transactions on computers, 100(4):353–359.

Kosterlitz and
Thouless (1973)
Kosterlitz, J. M. and Thouless, D. J. (1973).
Ordering, metastability and phase transitions in twodimensional systems.
Journal of Physics C: Solid State Physics, 6(7):1181.  Kühn et al. (1991) Kühn, R., Bös, S., and van Hemmen, J. L. (1991). Statistical mechanics for networks of gradedresponse neurons. Physical Review A, 43(4):2084.
 Kuramoto (1975) Kuramoto, Y. (1975). Selfentrainment of a population of coupled nonlinear oscillators. In International symposium on mathematical problems in theoretical physics, pages 420–422. Springer.
 Little (1974) Little, W. A. (1974). The existence of persistent states in the brain. In From HighTemperature Superconductivity to Microminiature Refrigeration, pages 145–164. Springer.
 London et al. (2010) London, M., Roth, A., Beeren, L., Häusser, M., and Latham, P. E. (2010). Sensitivity to perturbations in vivo implies high noise and suggests rate coding in cortex. Nature, 466(7302):123.
 Lynch et al. (2016) Lynch, G. F., Okubo, T. S., Hanuschkin, A., Hahnloser, R. H., and Fee, M. S. (2016). Rhythmic ContinuousTime Coding in the Songbird Analog of Vocal Motor Cortex. Neuron, 90(4):877–892.
 Mariño et al. (2005) Mariño, J., Schummers, J., Lyon, D. C., Schwabe, L., Beck, O., Wiesing, P., Obermayer, K., and Sur, M. (2005). Invariant computations in local cortical networks with balanced excitation and inhibition. Nature Neuroscience, 8(2):194–201.
 Merolla et al. (2014) Merolla, P. A., Arthur, J. V., AlvarezIcaza, R., Cassidy, A. S., Sawada, J., Akopyan, F., Jackson, B. L., Imam, N., Guo, C., Nakamura, Y., Brezzo, B., Vo, I., Esser, S. K., Appuswamy, R., Taba, B., Amir, A., Flickner, M. D., Risk, W. P., Manohar, R., and Modha, D. S. (2014). A million spikingneuron integrated circuit with a scalable communication network and interface. Science, 345(6197):668–673.
 Meunier et al. (1989) Meunier, C., Hansel, D., and Verga, A. (1989). Information processing in threestate neural networks. Journal of Statistical Physics, 55(56):859–901.
 Mirollo and Strogatz (1990) Mirollo, R. E. and Strogatz, S. H. (1990). Synchronization of pulsecoupled biological oscillators. SIAM Journal on Applied Mathematics, 50(6):1645–1662.
 Mueller and Herz (1999) Mueller, R. and Herz, A. (1999). Contentaddressable memory with spiking neurons. Physical Review E, 59(3):3330.
 Neckar et al. (2019) Neckar, A., Fok, S., Benjamin, B. V., Stewart, T. C., Oza, N. N., Voelker, A. R., Eliasmith, C., Manohar, R., and Boahen, K. (2019). Braindrop: A mixedsignal neuromorphic architecture with a dynamical systemsbased programming model. Proceedings of the IEEE, 107(1):144–164.
 Niebur et al. (1991) Niebur, E., Schuster, H. G., Kammen, D. M., and Koch, C. (1991). Oscillatorphase coupling for different twodimensional network connectivities. Physical Review A, 44(10):6895.
 Nischwitz and Glünder (1995) Nischwitz, A. and Glünder, H. (1995). Local lateral inhibition: a key to spike synchronization? Biological Cybernetics, 73(5):389–400.
 Noest (1988) Noest, A. J. (1988). Phasor neural networks. In Neural information processing systems, pages 584–591.
 O’Keefe and Recce (1993) O’Keefe, J. and Recce, M. L. (1993). Phase relationship between hippocampal place units and the hippocampal theta rhythm. Hippocampus, 3(3):317–330.
 Palm and Sommer (1992) Palm, G. and Sommer, F. T. (1992). Information capacity in recurrent mcculloch–pitts networks with sparsely coded memory states. Network: Computation in Neural Systems, 3(2):177–186.
 Panzeri and Diamond (2010) Panzeri, S. and Diamond, M. E. (2010). Information carried by population spike times in the whisker sensory cortex can be decoded without knowledge of stimulus time. Frontiers in Synaptic Neuroscience, 2(JUN):1–14.
 Pouget et al. (2003) Pouget, A., Dayan, P., and Zemel, R. S. (2003). Inference and computation with population codes. Annual review of neuroscience, 26:381–410.
 Reinagel and Reid (2000) Reinagel, P. and Reid, R. C. (2000). Temporal coding of visual information in the thalamus. Journal of Neuroscience, 20(14):5392–400.
 Riedel et al. (1988) Riedel, U., Kühn, R., and Van Hemmen, J. (1988). Temporal sequences and chaos in neural nets. Physical review A, 38(2):1105.
 Riehle et al. (1997) Riehle, A., Grun, S., Diesmann, M., Aertsen, A., Grün, S., Diesmann, M., Aertsen, A., Grun, S., Diesmann, M., and Aertsen, A. (1997). Spike Synchronization and Rate Modulation Differentially Involved in Motor Cortical Function. Science, 278(5345):1950–1953.
 Schuster and Wagner (1990) Schuster, H. and Wagner, P. (1990). A model for neuronal oscillations in the visual cortex. Biological cybernetics, 64(1):77–82.
 Schwenker et al. (1996) Schwenker, F., Sommer, F. T., and Palm, G. (1996). Iterative retrieval of sparsely coded associative memory patterns. Neural Networks, 9(3):445–455.
 Sirat et al. (1989) Sirat, G. Y., Maruani, A. D., and Chevallier, R. C. (1989). Grey level neural networks. Applied optics, 28(3):414–415.
 Sommer and Wennekers (2001) Sommer, F. T. and Wennekers, T. (2001). Associative memory in networks of spiking neurons. Neural networks, 14(67):825–834.
 Sompolinsky et al. (1990) Sompolinsky, H., Golomb, D., and Kleinfeld, D. (1990). Global processing of visual stimuli in a neural network of coupled oscillators. PNAS, 87(18):7200–7204.
 Sompolinsky and Kanter (1986) Sompolinsky, H. and Kanter, I. (1986). Temporal association in asymmetric neural networks. Physical Review Letters, 57(22):2861–2864.
 Swadlow (2000) Swadlow, H. A. (2000). Information flow along neocortical axons. In Miller, R., editor, Time and the Brain. Harwood Academic Publishers.
 Szatmáry and Izhikevich (2010) Szatmáry, B. and Izhikevich, E. M. (2010). Spiketiming theory of working memory. PLoS Computational Biology, 6(8).
 Thorpe et al. (1996) Thorpe, S., Fize, D., and Marlot, C. (1996). Speed of processing in the human visual system.
 Tiesinga et al. (2008) Tiesinga, P., Fellous, J.M., and Sejnowski, T. J. (2008). Regulation of spike timing in visual cortical circuits. Nature Reviews Neuroscience, 9(2):97–107.
 Treves (1990) Treves, A. (1990). Gradedresponse neurons and information encodings in autoassociative memories. Physical Review A, 42(4):2418.
 Tsodyks and Feigel’man (1988) Tsodyks, M. V. and Feigel’man, M. V. (1988). The enhanced storage capacity in neural networks with low activity level. EPL (Europhysics Letters), 6(2):101.
 Van Hemmen and Wreszinski (1993) Van Hemmen, J. L. and Wreszinski, W. F. (1993). Lyapunov function for the kuramoto model of nonlinearly coupled oscillators. Journal of Statistical Physics, 72(12):145–166.
 Welinder et al. (2008) Welinder, P. E., Burak, Y., and Fiete, I. R. (2008). Grid cells: The position code, neural network models of activity, and the problem of learning. Hippocampus, 18(12):1283–1300.
 Wennekers et al. (1995) Wennekers, T., Sommer, F. T., and Palm, G. (1995). Iterative retrieval in associative memories by threshold control of different neural models. Supercomputing in brain research: From tomography to neural networks, pages 301–319.
 Willshaw et al. (1969) Willshaw, D. J., Buneman, O. P., and LonguetHiggins, H. C. (1969). Nonholographic associative memory. Nature.
 Willwacher (1976) Willwacher, G. (1976). Capabilities of an associative storage system compared with the function of the brain. Biological Cybernetics, 24:181–198.
 Willwacher (1982) Willwacher, G. (1982). Storage of a Temporal Pattern Sequence in a Network. Biological Cybernetics, 43:115–126.
 Yedidia (1989) Yedidia, J. (1989). Neural networks that use threestate neurons. Journal of Physics A: Mathematical and General, 22(12):2265.
Comments
There are no comments yet.