PMCCPMCCPMCC

Conseils de recherche
Les critères de recherche 

Avancée

 
Logo of jmathneurSpringerOpen.comSubmit OnlineRegisterThis journalThis article
 
J Math Neurosci. 2017; 7: 9.
Published online 2017 August 25. doi:  10.1186/s13408-017-0051-7
PMCID: PMC5572897

An Analysis of Waves Underlying Grid Cell Firing in the Medial Enthorinal Cortex

Abstract

Layer II stellate cells in the medial enthorinal cortex (MEC) express hyperpolarisation-activated cyclic-nucleotide-gated (HCN) channels that allow for rebound spiking via an Ih current in response to hyperpolarising synaptic input. A computational modelling study by Hasselmo (Philos. Trans. R. Soc. Lond. B, Biol. Sci. 369:20120523, 2013) showed that an inhibitory network of such cells can support periodic travelling waves with a period that is controlled by the dynamics of the Ih current. Hasselmo has suggested that these waves can underlie the generation of grid cells, and that the known difference in Ih resonance frequency along the dorsal to ventral axis can explain the observed size and spacing between grid cell firing fields. Here we develop a biophysical spiking model within a framework that allows for analytical tractability. We combine the simplicity of integrate-and-fire neurons with a piecewise linear caricature of the gating dynamics for HCN channels to develop a spiking neural field model of MEC. Using techniques primarily drawn from the field of nonsmooth dynamical systems we show how to construct periodic travelling waves, and in particular the dispersion curve that determines how wave speed varies as a function of period. This exhibits a wide range of long wavelength solutions, reinforcing the idea that rebound spiking is a candidate mechanism for generating grid cell firing patterns. Importantly we develop a wave stability analysis to show how the maximum allowed period is controlled by the dynamical properties of the Ih current. Our theoretical work is validated by numerical simulations of the spiking model in both one and two dimensions.

Electronic Supplementary Material

The online version of this article (doi:10.1186/s13408-017-0051-7) contains supplementary material.

Keywords: Grid cell, Medial enthorinal cortex, h-current, Rebound spiking, Integrate-and-fire neural field model, Nonsmooth dynamics, Travelling wave, Evans function

Introduction

The ability to remember specific events occurring at specific places and times plays a major role in our everyday life. The question of how such memories are established remains an active area of research, but several key facts are now known. In particular, the 2004 discovery of grid cells in the medial enthorinal cortex (MEC) by Fyhn et al. [2], supports the notion of a cognitive map for navigation. This is a mental representation whereby individuals can acquire, code, store, recall, and decode information about relative spatial locations in their environment. The concept was introduced by Tolman in 1948 [3], with the first neural correlate being identified as the place cell system in the hippocampus [4]. Place cells are found in the hippocampus and fire selectively to spatial locations, thereby forming a place field whose properties change from one environment to another. More recently, a second class of cells was identified that fire at the nodes of a hexagonal lattice tiling the surface of the environment covered by the animal—these are termed grid cells [5]. As an animal approaches the centre of a grid cell firing field, the spiking output of grid cell will increase in frequency. The grid field size and spacing increases from dorsal to ventral positions along the MEC and is independent of the animal’s speed and direction (even in the absence of visual input) and independent of the arena size. In rats, the grid field spacing can range from roughly 30 cm up to several meters [6]. Other grid cell properties include firing field patterns that manifest instantly in novel environments and maintain alignment with visual landmarks. Furthermore, neighbouring grid cells have firing fields with different spatial phases, whilst grid cells with a common spacing also have a common orientation (overturning an original suggestion that they have different orientations) [7].

May-Britt Moser and Edvard Moser shared the 2014 Nobel Prize in Physiology or Medicine with John O’Keefe for their discoveries of cells (grid and place cells, respectively) that subserve the brain’s internal positioning system. From a modelling perspective grid cells have attracted a lot of attention, due in part to their relatively recent and unexpected discovery, but also due to the very geometric firing patterns that they generate. There are now three main competing mathematical models for the generation of grid-like firing patterns: oscillatory interference models, continuous attractor network models, and “self-organised” models—see Giocomo et al. [8] and Schmidt-Hieber and Häusser [9] for excellent reviews. The first class of model uses interference patterns generated by multiple oscillations to explain grid formation [10]. They have been especially fruitful in addressing the theta rhythmic firing of grid cells (5–12 Hz) and their phase precession. Here spikes occur at successively earlier phases of the theta rhythm during a grid field traversal, suggestive of a spike-timing code [11]. The second class uses activity in local networks with specific connectivity to generate the grid pattern and its spacing. In this category, the models can be further sub-divided into those that utilise spatial pattern formation across the whole tissue (possibly arising via a Turing instability), such as in the work of Fuhs and Touretzky [12] and Burak and Fiete [13], and those that rely only upon spatially localised pattern states (or bumps) in models with (twisted) toroidal connectivity as described by McNaughton et al. [14] and Guanella et al. [15]. The third class proposes that grid cells are formed by a self-organised learning process that borrows elements from both former classes [1618]. Recent experiments revealing the in vivo intracellular signatures of grid cells, the primarily inhibitory recurrent connectivity of grid cells, and the topographic organisation of grid cells within anatomically overlapping modules of multiple spatial scales along the dorsoventral axis of MEC provide strong constraints and challenges to all three classes of models [1820]. This has led to a variety of new models, each with a focus on one or more aspects of biophysical reality that might underlie the functionality of grid cell response. For example Couey et al. [21] have shown that a continuous attractor network with pure inhibition can support grid cell firing, with the caveat that there is sufficient excitatory input to the MEC, supposedly from hippocampus, to cause principal cells to fire. However, recent optogenetic and electrophysiological experiments have challenged this simple description [22], highlighting the importance of intrinsic nonlinear ionic currents and their distribution amongst the main cell types in MEC.

Stellate and pyramidal cells constitute the principal neurons in layer II of medial enthorinal cortex (MEC II) that exhibit grid cell firing. The former comprise approximately 70% of the total MEC II neural population and are believed to represent the majority of the grid cell population. Even before the discovery of grid cells stellate cells were thoroughly studied because of their rapid membrane time constants and resonant behaviour. Indeed, they are well known to support oscillations in the theta frequency range [23, 24]. Interestingly the frequency of these intrinsic oscillations decreases along the dorsal-ventral axis of MEC II [25], suggestive of a role in grid field spacing. The resonant properties of stellate cells have been directly linked to a high density of hyperpolarisation-activated cyclic-nucleotide-gated (HCN) channels [26], underlying the so-called Ih current. The time constant of both the fast and slow component of Ih is significantly faster for dorsal versus ventral stellate cells, providing a potential mechanism for the observed difference in the resonant frequency along the dorsal-ventral axis. However, perhaps of more importance is the fact the Ih current can cause a depolarising rebound spike following a hyperpolarising current injection. Given that stellate cells are mainly interconnected by inhibitory interneurons [21], this means that rebound can play an important role in shaping spatio-temporal network rhythms. The inclusion of important intrinsic biophysical properties into a network has been emphasised by several authors, such as Navratilova et al. [27] regarding the contribution of after-spike potentials of stellate cells to theta phase precession, and perhaps most notably by Hasselmo and colleagues for the inclusion of HCN channels [2831]. This has culminated in a spiking network model of MEC that supports patterns whose periodicity is controlled by a neuronal resonance frequency arising from an Ih current [1]. The model includes many of the features present in the three classes described above, and is able to replicate behaviour from several experiments, including phase precession in response to a phasic medial septum input, theta cycle skipping, and the loss of the spatial periodicity of grid cell firing fields upon a reduction of input from the medial septum. Simulations of the model in one spatial dimension show that the spacing of grid firing fields can be controlled by manipulating the speed of the rebound response. We note that a change that affects the rebound response would also affect the resonant properties of the cell. In contrast to continuous attractor models that rely on the spatial scale of connectivity to control grid spacing, a change in rebound response provides a mechanism of local control via changes in the expression levels of HCN channels. This fascinating observation warrants a deeper mathematical analysis. In this paper we introduce a spiking network model that shares many of the features of the Hasselmo model [1], focussing on the formation of travelling waves that can arise in the absence of (medial septum) input. Importantly our bespoke model is built from piecewise linear and discontinuous elements that allows for an explicit analysis of the periodic waves that arise in an inhibitory network through rebound spiking. In particular our wave stability analysis demonstrates clearly that the maximum allowed period is strongly controlled by the properties of our model Ih current. This gives further credence to the hypothesis that HCN channels can control the properties of tissue level periodic waves that may underpin the spacing of grid cell firing in MEC.

In Sect. 2 we introduce our model of a network of stellate cells in MEC II. The single neuron model is a simple leaky integrate-and-fire (LIF) model with the inclusion of a synaptic current mediated by network firing events, and a model of Ih based on a single gating variable. For ease of mathematical analysis we focus on a continuum description, so that the model may be regarded as a spiking neural field. Simulations of the model in two spatial dimensions are used to highlight the genericity of spike rebound mediated spatio-temporal patterns. To uncover the systematic way in which a cellular Ih current can control the properties of patterns at the tissue level, we focus next on a one-dimensional version of the model without external input. By further developing a piecewise linear (pwl) caricature of the activation dynamics of a HCN channel we show in Sect. 3 how explicitly to construct the dispersion curve for periodic travelling waves. This gives the speed of a wave as a function of its period, and shows the possibility of a wide range of wave periods that could be selected. Next in Sect. 4 we exploit techniques from the analysis of nonsmooth systems to determine the Evans function for wave stability. Importantly an investigation of this function shows that the maximum allowed period can be controlled both by the overall conductance strength of the Ih current as well as the time-scale for activation of HCN channels. Finally in Sect. 5 we discuss natural extensions of our work, and its relevance to further models of grid cell firing.

The Model

The original work of Hasselmo [1] considered both simple resonant models as well as spiking nonlinear integrate-and-fire Izhikevich units to describe MEC stellate cells. The former incorporated a model for Ih using a single gating variable with a linear dynamics whilst the latter were tuned to capture the resonant and rebound spiking properties from experimental data. Synaptic interactions between cells in a discrete one-dimensional network were modelled with a simple voltage threshold process. These models were subsequently studied in more detail in [28], with a further focus on two-dimensional models and travelling waves. Here, we consider a biophysically realistic spiking model in a similar spirit to that of Hasselmo, but within a framework that will allow for a subsequent mathematical analysis. In particular, we consider a spiking model of stellate cells using a generalised LIF model that includes a nonlinear Ih current. Moreover, we opt for a description of synaptic interactions using an event-based scheme for modelling post-synaptic conductances.

We first consider a continuum description defined on the plane and introduce a voltage variable VV(rt), where r ∈ ℝ2 and t ≥ 0. The subthreshold LIF dynamics is given by

CtV(r,t)=glV(r,t)+Ih(r,t)+Isyn(r,t)+Ihd(r,t),
1

with a set of spike times at position r generated according to

Tm(r) = inf{tV(rt) ≥ Vtht ≥ Tm−1(r) + τR},  m ∈ ℤ.
2

Here τR is a refractory time-scale. Upon reaching the threshold Vth the membrane potential is reset to the value Vr < Vth. The infimum operation ensures that a firing time is determined by the first time that the voltage variable (at a fixed position) crosses threshold (from below, remembering the IF reset process) subject to refractoriness. To model the refractory process we hold the voltage variable at the reset value Vr for a duration τR after a firing event. The left-hand side of (1) describes a membrane current with capacitance C. The first term on the right-hand side of (1) represents a leak with a constant conductance gl (and we have set the leakage reversal potential to zero without loss of generality). The terms Ih, Isyn, and Ihd represent currents arising from HCN channels, synaptic input, and head-direction input, respectively. Ih is a slow inward current with a reversal potential Vh that is substantially above resting levels, but which requires hyperpolarisation to become active; that is, the activation curve is monotone decreasing in V. Furthermore, the Ih current does not inactivate, even with prolonged (minutes) hyperpolarisation. Thus it is often modelled with a single gating variable nh such that Ih(rt) = ghnh(rt)(Vh − V(rt)), where

τh(V)tnh=nh,(V)nh.
3

Here the shape for the activation function nh,∞ is the sigmoid:

nh,(V)=11+exp((VV1/2)/k),
4

and fits to experimental data give V1/2 ≈ −10 mV (with respect to rest) and k ≈ 10 [32, 33]. The time constant for activation and deactivation can vary from tens to hundreds of milliseconds [34]; however, here we fix τh = constant for simplicity and ignore any detailed dependence on voltage. To model synaptic interactions we consider a simple effective anatomical model whereby stellate cells interact directly through inhibition. In reality inhibitory interactions are actually mediated by interneurons and there is no direct synaptic coupling between stellate cells. Introducing an overall strength of synaptic conductance gsyn we then write Isyn(rt) = gsynψ(rt), where

ψ(rt) = ∫ΩdrW(rr)E(rt), 
5

where the function W represents anatomical connectivity, Ω ⊆ ℝ2 is the spatial domain, and the function E represents the shape of a post-synaptic response to a train of incoming spikes. We write this in the form

E(r,t)=mZη(tTm(r)),
6

for a given temporal filter shape η with the property that η(t) = 0 for t < 0 (so that interactions are causal). For convenience we will work with normalised responses such that 0dtη(t)=1. As a concrete choice for the function W we shall take a smoothed bump shape W(rr) = w(|r − r|), with

w(x)=w02[tanh(β(σx))+tanh(β(σ+x))],β,σ>0.
7

Here w0 < 0 in accordance with the predominantly inhibitory interactions of MEC, σ controls the spatial scale of interaction, and β the steepness of the surround inhibition function as shown in Fig. 1. The model is completed with the choice of the synaptic filter η, which we shall take to be an α-function of the form

η(t) = α2teαtH(t), 
8

where α−1 is the time-to-peak of the synapse, and H is the Heaviside step function. Note that we work in a regime where the model is excitable, as we do not include any background drive that would be able to make the single neuron model fire in the absence of synaptic coupling or head-direction input.

Fig. 1
Connectivity function w in (7), with σ = 25, β = 0.5 and w0 = −10

An animal’s speed vv(t) and direction of motion ΦΦ(t) generates input to the MEC that is modelled by the head-direction current Ihd. For example, this could be of the form Ihd(rt) = v ⋅ rϕ where vv(cosΦ, sinΦ) and rϕl(cosϕ(r), sinϕ(r)) describes a head-direction preference for the orientation ϕ(r) at position r [13]. Here l is a constant that sets the magnitude of the head-direction current. The hidden assumption here is that head direction matches the direction of motion. However, this is not necessarily true behaviourally, and head direction cells may not code for motion direction [35]. Nonetheless it is a common assumption in most grid cell models, and so we adopt it here too. Fuhs and Touretzky [12] have shown that choosing an anisotropic anatomical connectivity of the particular form W(rr) = w(|r − rϕ − r|) can then induce a spatio-temporal network pattern to flow in accordance with the pattern of head-direction information generated when traversing an environment. For continuous attractor network models that can generate, via a Turing instability, static hexagonal patterns with regions of high activity at the nodes of a triangular lattice, the induced movement of these hot-spots over a given point in the tissue gives rise to grid-like firing patterns. In this case the spacing of the firing field is hard to change, as it is mainly fixed by the spatial scale of the chosen connectivity; however, the mechanism of inducing pattern flow is robust to how the tissue pattern is generated. Thus, given the established success of the Fuhs and Touretzky mechanism we will not focus on this here, and instead turn our attention to the formation of relevant tissue patterns and, in particular, how local control of firing field spacing may be effected.

To investigate the types of solution the spiking neural field model supports, we perform simulations over a two-dimensional square domain. Since the action of the head-direction current is merely to induce a flow of emergent patterns we restrict our attention to the case without such input, i.e. Ihd = 0. See Appendix A for details of our bespoke numerical scheme implemented on a GPU architecture, and Additional file 1 for C++/CUDA code. We observe three general classes of coherent behaviour that take the form of spatially periodic non-travelling structures (though which oscillate in time), travelling periodic waves, and lurching waves. The latter are also generically found in neural systems with rebound currents, such as in models of thalamic slices [3638]. Unlike traditional smoothly propagating waves, which exhibit a stationary profile in a co-moving frame, lurching waves consist of patterns of activity in a localised region of the domain, which after some period of time, decay and an adjacent region of the domain becomes active. These waves appear to ‘lurch’ across the domain. Whilst interesting in their own right, we focus in this article on the analysis of smoothly propagating periodic waves, since these have been suggested by Hasselmo [1] to play a dominant role in the formation of grid-like firing patterns in MEC. We show an example of a non-travelling periodic structure in Fig. 2, at two different time points to illustrate that the pattern is not static, but oscillates in time. An example of a smoothly propagating travelling wave is shown in Fig. 3. The model also supports more exotic spatio-temporal structures, including hexagonal patterns, and for a movie showing a dynamic state with a hexagonal sub-structure see Additional file 2.

Fig. 2
A simulation of spatially periodic non-travelling patterns in a two-dimensional spiking neural field model with an Ih current, solved on a spatial grid of 1000 × 1000 points. Displayed is the voltage component across the entire network at t = 7000 ms ...
Fig. 3
Periodic travelling wave solutions in the spiking neural field model, solved on a spatial grid of 1000 × 1000 points. Displayed is the voltage component across the entire network at t = 7000 ms (left) and t=10,000 ms (right). Using different ...

Wave Construction

To understand more fully how Ih controls the emergent scale of periodic waves seen in the simulations of Sect. 2 we now turn to a one-dimensional version of the model defined on the infinite domain. As in Sect. 2, we consider Ihd = 0, in which case the model is isotropic and (5) reduces to

ψ(x,t)=mZdyw(|xy|)η(tTm(y)),xR,t>0.
9

To reduce the model to a more mathematically convenient form we make two observations about the Ih current. The first is that Vh is typically larger than V, which suggests the approximation Vh − V ≃ Vh. The second is that the nonlinear activation function nh,∞ can be approximated by a pwl function, as illustrated in Fig. 4. Here we match the slope at VV1/2, and otherwise saturate the function to one or zero, so that

nh,(V)={1,VV12VV1/24k,V<V<V+0,VV+,,V±=V1/2±2k.
10

Fig. 4
Red line: Nonlinear activation function for nh with V1/2 = −10 mV and k = 10. Green line: Piecewise linear fit of nh,∞ given by Eq. (10)

Simulations of the model with the reduced form of the Ih are in qualitative agreement with simulations of the full nonlinear model, and indeed wherever tested the same repertoire of behaviour is always found. In both instances, travelling wave behaviour with a well-defined speed and period is easily initiated; Figs. 5 and and66 compare simulation results arising in the fully nonlinear and reduced pwl model.

Fig. 5
Simulations of a one-dimensional spiking neural field model with an Ih current solved on a spatial grid of 5000 points. Left: simulations of the model with the full nonlinear dynamics for Ih. Right: simulations of the model with the reduced pwl dynamics ...
Fig. 6
Voltage traces observed at x = 0, when the periodic travelling wave is fully developed. Top: simulations of the model with the reduced pwl dynamics for Ih given by (10). Bottom: simulations of the model with the full nonlinear dynamics for Ih. All parameters ...

If we introduce the vector X = (Vnh) ∈ ℝ2 then we may write the reduced model in a more abstract setting, namely in terms of the pwl evolution equation that governs the system behaviour between one spiking event and the next:

tX(x,t)=AX(x,t)+Ψ(x,t),Tm1(x)t<Tm(x).
11

In (11), A is a 2 × 2 matrix that is defined in a piecewise constant fashion, with dependence on the value of the voltage (in particular, which of the three domains detailed in equation (10) pertains), or whether the system is in the refractory state. In the latter case, A is defined according to

A=AR=[0001/τh],Tm1(x)t<Tm1(x)+τR,
12

while for Tm−1(x) + τR ≤ t < Tm(x)

A={A0=[1/τghgl1/τ1/(4kτh)1/τh],V<V<V+,A=[1/τghgl1/τ01/τh],VV,A+=A,VV+,
13

where we have assumed the ordering V < Vr < V+ < Vth, introduced τC/gl, and absorbed the factor Vh within gh. Similarly we define Ψ according to

Ψ=ΨR=[0(1/2(VrV1/2)/(4k))/τh],Tm1(x)t<Tm1(x)+τR,
14

and, for Tm−1(x) + τR ≤ t < Tm(x),

Ψ=gsyngl1τψ[10]+{b0,V<V<V+,b,VV,b+,VV+,
15

where

b0=[0(1/2+V1/2/(4k))/τh],b=[01/τh],b+=[00].
16

We highlight that in (12)–(16) we have introduced the subscripts {R, 0, −,  + } to indicate the state of the system, namely whether it is refractory (labelled by R) or is not refractory and has a voltage in the range (VV+) (labelled by 0), (−∞, V] (labelled by −), [V+, ∞) (labelled by +).

Travelling Wave Analysis

We now seek travelling wave solutions of (11) of the form Xˆ(ξ,t), where ξt − x/c and c is the (constant) wave speed. In this case (11) transforms to

(t+ξ)Xˆ(ξ,t)=AXˆ(ξ,t)+Ψˆ(ξ,t).
17

A stationary travelling wave Xˆ(ξ,t)=Q(ξ)=(V(ξ),nh(ξ)) satisfies the travelling wave equation

dQdξ=AQ(ξ)+Ψˆ(ξ).
18

In terms of firing events a periodic wave is described by Tm(x) = x/cmΔ, where Δ is the period of the wave such that Q(ξΔ) = Q(ξ). Substitution of this firing ansatz into (9) allows us to determine the function ψˆ(ξ)=ψ(x,t)|t=Tm(x) where Ψˆ(ξ) is obtained from (15) under the replacement of ψ by [psi]. The function [psi] is easily determined as

ψˆ(ξ)=cmZ0dsη(s)w(|c(sξ)+cmΔ|),
19

and is Δ-periodic, and can therefore be expressed in terms of a Fourier series as

ψˆ(ξ)=pZψpe2πipξ/Δ,ψp=1Δw˜(2πpcΔ)η˜(2πpΔ).
20

In (20) tildes denote the Fourier integral representation, such that for a given function a(x)

a(x)=12πdka˜(k)eikx,a˜(k)=dxa(x)eikx,
21

and we have made use of the result that cΔmeikcmΔ = 2πpδ(k − 2πp/(cΔ)).

For the bump function (7) and the α-function (8) we have

w˜(k)=w0πβsinkσsinh(πk/(2β)),η˜(k)=α2(α+ik)2.
22

Thus, given the decay properties of (22) as a function of k, the sum in (20) can be naturally truncated.

The formal solution to (18) can be constructed using variation of parameters as

Q(ξ)=G(ξ,ξ0)Q(ξ0)+ξ0ξG(ξ,ξ)Ψˆ(ξ)dξ,
23

where G is a matrix exponential given by

G(ξ,ξ)=T{exp(ξξdsA(s))}.
24

Here 𝒯 is a time-ordering operator 𝒯A(t)A(s) = H(t − s)A(t)A(s) + H(s − t)A(s)A(t). In general the issue of time-ordering makes it very difficult to evaluate G. However, in our case A is piecewise constant and so we easily may break the solution up into parts distinguished by the label μ ∈ {R, 0, −,  + }. In each case trajectories are given explicitly by (23) with G(ξξ) = G(ξ − ξ) and GGμ where Gμ(ξ) = exp(Aμξ). A global trajectory may then be obtained by patching together solutions, denoted by Qμ, from each domain. It is in this fashion that we now construct the shape of a periodic travelling wave in a self-consistent manner. Of use will be matrix exponential decomposition eAtPeΛtP−1, where Λ = diag(λ+λ) are the eigenvalues of A with associated eigenvectors q± = (1,(λ±A11)/A12,)T. Here the eigenvalues of A are given explicitly by

λ±=12(TrA±(TrA)24detA).
25

Using (20) we may then write a domain specific trajectory for μ ∈ {0,  + , −} in the form

Qμ(ξ)=Gμ(ξξ0)Qμ(ξ0)+Aμ1[Gμ(ξξ0)I2]bμ+gsyngl1τpZψpPμdiag(Zμ+(ξ,ξ0),Zμ(ξ,ξ0))Pμ1[10],
26

where

Zμ±(ξ,ξ0)=eλμ±(ξξ0)e2πipξ0/Δe2πipξ/Δλμ±+2πip/Δ.
27

Here Λμ=diag(λμ+,λμ) with λμ± representing the eigenvalues of Aμ and

Pμ=[11(λμ+[Aμ]11)/[Aμ]12(λμ[Aμ]11)/[Aμ]12].
28

When μR we adopt an alternative strategy (since AR is singular), remembering that when the system is refractory then V(ξ) is clamped at the value VVr. In this case, we only have to consider the evolution of the gating variable nh, which is obtained from (3) and (10) to give

nh(ξ) = nh(ξ0)e−(ξξ0)/τh + (1/2 − (Vr − V1/2)/(4k))[1 − e−(ξξ0)/τh].
29

Now let us consider the form of a periodic wave which elicits a single spike for every period, much like the ones seen in Fig. 5. An example of such a travelling wave orbit in the (Vnh) phase plane is shown in Fig. 7. The corresponding evolution of V(ξ) and nh(ξ) is shown in Fig. 8.

Fig. 7
Orbit of a periodic travelling wave. Parameters as in Fig. 5, with Δ = 450, c = 0.0669, nh(0) = 0.3815 and ξ1 = 225.4223. The travelling wave starts at ξ = 0 where VVr is clamped and nh evolves according to (29) until ξτR and the system ...
Fig. 8
Profile of the two components of the periodic travelling wave Q(ξ) defined by (26) and illustrated in Fig. 7. Solid lines correspond to V(ξ) (left-hand axis) and dotted ones to nh(ξ) (right-hand axis). Colour-code and parameters as in Fig.  ...

Given the translational invariance of the system we are free to choose a travelling wave origin such that ξ = 0 corresponds to the system immediately after firing. For a duration τR it will then remain clamped at Vr with nh evolving according to (29) with ξ0 = 0 (as shown in blue in Fig. 7 and Fig. 8). From here it then evolves according to (26) with μ = 0, with initial data determined by QR(ξ0) = (Vrnh(τR)), until V(ξ) reaches V±, after which we set μ = − or μ =  +  in (26) and select appropriate initial data for (26), depending on the value of V achieved first. For simplicity, and since this is reliably observed in numerical simulations for a wide range of parameters (red line in Fig. 7 and Fig. 8, though the argument is easily generalised), we assume V+ is the relevant choice. The final piece of the orbit is then obtained from (26) with μ =  +  and initial data determined by Q+(ξ0) = Q0(ξ1τR) (green line in Fig. 7 and Fig. 8) and evolving the system until V(ξ) = Vth. Denoting the time of flight for the trajectory such that Vr ≤ V < V+ by ξ1, and that, for V+ ≤ V < Vth by ξ2, the period of the orbit is given by ΔτRξ1ξ2. We note that the orbit is discontinuous because of the reset of the voltage variable after one period. Thus we have four unknowns (nh(0), cξ1Δ) related by three nonlinear algebraic equations

{V(Δ)=Vth(firing condition),nh(Δ)=nh(0)(periodicity condition),V(τR+ξ1)=V+(switching condition),
30

whose simultaneous solution determines the dispersion relationship for the wave speed as a function of the period cc(Δ). An example of a dispersion curve constructed in this way is shown in Fig. 9. Here we see that a wide range of allowed wavelengths can co-exist (with differing speeds). Note that in Fig. 9 we also include solutions that visit the region of phase-space where V < V, and these solutions typically only occur for small values of Δ. Our constructive theory does not provide a wave selection principle; however, by varying initial data in direct numerical simulations we were able to find solutions in excellent agreement with the theoretical predictions up to some maximal value of Δ. The determination of this value is the subject of the next section, where we show how to analyse wave stability.

Fig. 9
Dispersion curve cc(Δ) for a periodic travelling wave. Here c is the speed of a wave with period Δ. For small periods and on the upper branch waves are constructed that visit the domain where V < V. Parameters as in Fig. 5. Solid lines ...

Wave Stability

To determine the stability of a periodic travelling wave we must not only treat perturbations of the state variables, but also the associated effects on the times of firing. Moreover, one must remember that because the model is nonsmooth (due to the switch at VV±) and discontinuous (because of reset whenever VVth) standard approaches for analysing smooth dynamical systems cannot be immediately applied. Nonetheless, as we show below, the wave stability can in fact be explicitly determined. We do this by constructing the so-called Evans function. This has a long history of use in the analysis of wave solutions to partial differential equations, dating back to the work of Evans on the stability of action potentials in the Hodgkin–Huxley model of a nerve fibre [39], has been extended to certain classes of firing rate neural field model [40], and is developed here for spiking neural fields.

We begin our analysis by exposing the spike-train that determines the synaptic drive in (9) by writing it in the equivalent form

ψ(x,t)=mZdyw(|xy|)tdsη(ts)δ(sTm(y)),
31

where the firing times are defined according to the threshold condition V(xTm(x)) = Vth. We may relate spike times to voltage threshold conditions using the result that for fixed x

δ(t − Tm(x)) = |Vt(xTm(x))|δ(V(xt)−Vth), 
32

and Vt denotes partial differentiation of V with respect to t. Hence

ψ(x,t)=mZdyw(|xy|)×tdsη(ts)|Vt(y,Tm(y))|δ(V(y,s)Vth).
33

Consider again travelling wave solutions of the form V(x,t)=Vˆ(ξ,t), where ξt − x/c. In this co-moving frame we can define a set of firing event functions ξm(t) according to the threshold condition Vˆ(ξm(t),t)=Vth. These event times can be related to the co-moving voltage threshold condition using the result that, for fixed t,

δ(ξξm(t))=|Vˆξ(ξm(t),t)|δ(Vˆ(ξ,t)Vth).
34

Substitution into (33) and using VˆξVt close to a periodic orbit we find

ψ(x,t)=mZdyw(|xy|)tdsη(ts)δ(sy/cξm(s))=cmZ0dsη(s)w(|c(sξ)+cξm(ts)|)ψˆ(ξ,t).
35

Noting that, for a periodic wave ξm(t) = mΔ, [psi] is independent of t and Eq. (35) recovers Eq. (19) as expected.

We now analyse the stability of such a periodic wave by perturbations such that ξm(t) = mΔδξm(t), with |δξm(t)| ≪ 1. Writing the corresponding perturbation of ψˆ(ξ,t) as ψˆ(ξ,t)=ψˆ(ξ)+δψˆ(ξ,t) we find

δψˆ(ξ,t)=c2mZ0dsη(s)w(|c(sξ)+cmΔ|)δξm(ts).
36

It remains to determine the relationship between δξm(t) and the perturbations of the shape of the travelling wave. In Appendix B we show that we can relate δξm(t) to the deviation in the voltage, denoted by δVˆ(mΔ,t), via the simple relationship

δξm(t)=δVˆ(mΔ,t)Vξ(mΔ).
37

Thus for solutions of the form δVˆ(ξ,t)=δVˆ(ξ)eλt, δVˆ(ξ)=δVˆ(ξ+Δ) we find δψˆ(ξ,t)=δψˆ(ξ;λ)eλt, with δψˆ(ξ;λ)=δVˆ(0)f(ξ;λ), and

f(ξ;λ)=c2Vξ(0)mZ0dsη(s)w(|c(sξ)+cmΔ|)eλs.
38

Returning to the more abstract setting given by Eq. (11) we linearise around the travelling wave by setting Xˆ(ξ,t)=Q(ξ)+δXˆ(ξ)eλt, with δXˆ(ξ)=δXˆ(ξ+Δ). This yields the variational equation

ddξδXˆ(ξ)=A(ξ;λ)δXˆ(ξ)+δΨˆ(ξ;λ),δΨˆ(ξ;λ)=gsyngl1τδψˆ(ξ;λ)[10],
39

where A(ξλ) = A(Q(ξ))−λI2 (and we use the argument of A to emphasise that it depends on position along the periodic orbit). We may write the solution to (39) in much the same way as for the periodic orbit problem given by (18), namely with the use of a variation of parameters formula, matrix exponentials and (38):

δXˆ(ξ)={GR(ξ;λ)δXˆ(0)0ξ<τR,G0(ξτR;λ)δXˆ(τR)+τRξdξG0(ξξ;λ)Jf(ξ;λ)δXˆ(0)τRξ<τR+ξ1,G+(ξ(τR+ξ1);λ)δXˆ(τR+ξ1)+τR+ξ1ξdξG+(ξξ;λ)Jf(ξ;λ)δXˆ(0)τR+ξ1ξ<Δ.
40

Here Gμ(ξλ) = exp([Aμ − λI2]ξ) and

J=gsyngl1τ[1000].
41

However, the evolution of perturbations through the switching manifolds VV±, the firing threshold VVth and the release from the refractory state requires care, since in all these cases there is a jump in the Jacobian. The theory of nonsmooth dynamical systems gives a prescription for handling this using so-called saltation matrices dating back to the work of Aizerman and Gantmakher in the 1950s [41]. For a more recent perspective we recommend the paper by Leine et al. [42] and the book by di Bernardo et al. [43], particularly in engineering applications, and the paper by Coombes et al. [44] for applications in neuroscience. The 2 × 2 saltation matrices for handling switching, firing, and refractoriness are constructed in Appendix C and denoted Kswitch, Kfire, and Kref, respectively. In essence they map perturbations through the region of nonsmooth behaviour to give δXˆ(0+)=KfireδXˆ(0), δXˆ(τR+)=δXˆ(τR)+KrefδXˆ(0), and δXˆ((τR+ξ1)+)=KswitchδXˆ(τR+ξ1). The saltation matrices are given explicitly by KswitchI2 and

Kfire=[00(nh,ξ(0+)nh,ξ(0))/Vξ(0)1],Kref=[Vξ(τR+)/Vξ(0)000].
42

If we now introduce the function μ(ξξ0λ):

Fμ(ξ,ξ0;λ)=ξ0ξdξGμ(ξξ;λ)Jf(ξ;λ),μ{0,+,},
43

then Eq. (40) may be used to generate the perturbation after one period as δXˆ(Δ)=Γ(λ,Δ)δXˆ(0), where

Γ(λ,Δ)=F+(Δ,τR+ξ1;λ)+G+(Δ(τR+ξ1);λ)Kswitch[F0(τR+ξ1,τR;λ)+G0(ξ1;λ)[GR(τR;λ)Kfire+Kref]].
44

Enforcing that perturbations be Δ-periodic (i.e. δXˆ(Δ)=δXˆ(0)) we obtain the spectral condition ℰ(λΔ) = 0 where

ℰ(λΔ) = |Γ(λΔ)−I2|.
45

We identify (45) as the Evans function for the periodic wave. To determine (43) in a computationally useful form we use a Fourier representation to represent (38) (cf. (20) from (19)) as f(ξλ) = ∑p∈ℤfp(λ)exp(−2πipξ/Δ) where

fp(λ)=1Vξ(0)2πΔ2ipη˜(iλ2πp/Δ)w˜(2πp/(cΔ)),
46

for Re(λα) > 0. Then in a similar way to the construction of (26) we find the useful representation for (43) as

Fμ(ξ,ξ0;λ)=pfp(λ)Pμdiag(Sμ+(ξ,ξ0;λ),Sμ(ξ,ξ0;λ))Pμ1J,
47

where

Sμ±(ξ,ξ0;λ)=e(λμ±λ)(ξξ0)e2πipξ0/Δe2πipξ/Δλμ±λ+2πip/Δ.
48

The eigenvalues of the spectral problem can be practically constructed by considering the decomposition λνiω and simultaneously solving the pair of equations 𝒢(νω) = 0 and ℋ(νω) = 0, where 𝒢(νω) = Reℰ(νiωΔ) and ℋ(νω) = Imℰ(νiωΔ), subject to the constraint να > 0. Figures 10 and and1111 show the zero level sets of 𝒢 and in the (νω) plane for two different points on the dispersion curve of Fig. 9. The intercepts when να > 0 provide the zeros of the Evans function and here highlight clearly that, for Δ sufficiently large, the zeros of the Evans function can cross to the right-hand complex plane signalling a wave instability.

Fig. 10
Zeros of the Evans function (45) with Δ = 460. These occur at the intersection (green dots) of 𝒢(νω) = 0 (red curve) and ℋ(νω) = 0 (blue curve) where 𝒢 is the real part of whereas is the imaginary part. Here we can see that all the eigenvalues, ...
Fig. 11
Zeros of the Evans function (45) with Δ = 470. These occur at the intersection (green dots) of 𝒢(νω) = 0 (red curve) and ℋ(νω) = 0 (blue curve) where 𝒢 is the real part of whereas is the imaginary part. Here we see a complex conjugate pair of ...

Figure 10 also shows that there is always a zero eigenvalue. It is simple to establish the persistence of this zero under parameter variation. Differentiating (18) with respect to ξ gives

ddξdQdξ=AdQdξ+ddξΨˆ(ξ).
49

From (38) and (19) we may establish that

δψˆ(ξ;0)=δVˆ(0)Vξ(0)ddξψˆ(ξ).
50

Hence for λ = 0 we see that a solution to (39) is the eigenfunction

δXˆ(ξ)=dQdξ,
51

as expected from translation invariance of the system (so that a perturbation tangential to the travelling wave orbit is neutrally stable).

Discussion

This paper is motivated by recent work in computational neuroscience that has highlighted rebound firing as a mechanism for wave generation in spiking neural networks that can underlie the formation of grid cell firing fields in MEC [1]. We have presented a simple spiking model with inhibitory synaptic connections and an Ih current that can generate smoothly propagating activity waves via post-inhibitory rebound. These are qualitatively of the type observed in previous computational studies [36, 37], yet are amenable to an explicit mathematical analysis. This is possible because we have chosen to work with piecewise linear discontinuous models, and exploited methodologies from the theory of nonsmooth systems. In particular we have exploited the linearity of our model between events (for firing, switching, and release from a refractory state) to construct periodic solutions in a travelling wave frame. To assess the stability of these orbits we have treated the propagation of perturbations through event manifolds using saltation operators. Using this we have constructed dispersion curves showing a wide range of stable periods, with a maximum period controlled by the time-scale of the rebound current. This gives further credence to the idea that the change in grid field scale along the dorsal-ventral axis of the MEC is under local control by HCN channels.

A number of natural extensions of the work in this paper suggest themselves. We briefly outline them here, and in no particular order. For simplicity we have focussed on the analysis of waves in a homogeneous model with only one spatial dimension. The analysis of the corresponding travelling waves, with hexagonal structure, in two spatial dimensions is more challenging, though is an important requirement for a complete model of grid cell firing. Moreover, the assumption of homogeneity should be relaxed. In this regard it would be of interest to understand the effects of a heterogeneity in the time constants (for voltage response, synaptic time-scale, and the time-scale of the Ih current) on the properties of spatio-temporal patterns. Furthermore, it would be biologically more realistic to consider two sets of inhibitory interneurons, as in [1, 28, 30]. As well as depending on the Ih current, grid field spacing changes as a function of behavioural context. This is believed to occur through the activation of neuromodulators [32], and simple regulation of our Ih current model would allow a systematic study of this, even before considering the more important issue of structured input. The work in this paper has focussed on spontaneous pattern formation that occurs in the absence of such input. Given the importance of the head-direction input for driving grid cell firing fields it would be natural to consider a mathematical analysis for the case with Ihd ≠ 0. For the standard Fuhs–Touretzky mechanism of inducing firing patterns this would further require the treatment of an anisotropic interaction kernel with a dependence on a head-direction preference map. One way to address this would be via a perturbation theory around the limiting case treated in this paper, and use this to calculate a tissue response parametrised by the animal’s speed and direction of motion. The same methodology would also allow an investigation of phase precession during a grid field traversal. Our simulations have also shown the possibility of ‘lurching waves’ and it would be interesting, at least from a mathematical perspective, to analyse their properties (speed and stability) and compare them to the co-existing smoothly propagating waves. It would further be pertinent in this case to pay closer attention to any possible nonsmooth bifurcations that could give rise to wave instabilities, such as grazing bifurcations. Finally we note that grid cells are grouped in discrete modules with common grid spacing and orientation [20]. It has recently been suggested that coupling between modules or via feedback loops to the hippocampus may help to suppress noise and underpin a robust code (with a large capacity) for the representation of position [45]. Another extension of the work in this paper would thus be to consider the dynamics of networks of interacting modules, paying closer attention to the details of MEC microcircuitry [46].

Electronic Supplementary Material

Acknowledgements

SC was supported by the European Commission through the FP7 Marie Curie Initial Training Network 289146, NETT: Neural Engineering Transformative Technologies. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Tesla K40 GPU used for this research. KCAW was generously supported by the Wellcome Trust Institutional Strategic Support Award (WT105618MA). MBQ gratefully acknowledge the support of CONACyT and The University of Nottingham during her PhD studies. We thank the anonymous referees for their helpful comments on our manuscript.

Appendix A: Numerical Scheme

Here, we describe the numerical scheme that we have developed to evolve the spiking neural field model. Given the large computational overheads for simulating the latter we have focussed on implementing the algorithm on a GPU system. In essence, we exploit the piecewise linear nature of the dynamics to obtain trajectories in closed form and then use a root-finding scheme to find the timings of events (switching, firing or refractory). The dynamics is updated at an event (generating synaptic currents and resetting at firing events, changing the gating dynamics at switching events, and releasing from reset at refractory events), and the process repeated. For clarity we describe our approach for two spatial dimensions as it is easily adapted to treat just one spatial dimension.

A.1 Algorithm

We solve the system over a 2D discretised grid of size L × L with N (numerical) cells in each spatial dimension and enforce periodic boundary conditions. Cells are equi-spaced with a spatial separation of Δx = 2L/N and the state of the cells in the network are evaluated over T + 1 time steps, discretising time as tiiΔt, i = 0, …, T. To evaluate these states, we analytically integrate (11) to provide a closed-form expression for the trajectory of the system. Our closed-form expression allows us, for cell j, to write its state at time ti+1 as a function of its state at time ti:

Xj(ti+1) = φΔt(tiXj(ti)), 

where φh is the evolution operator acting over a time h. Note that this expression assumes that no firing events have occurred between ti and ti+1, but can account for switching events, as detailed below. We shall describe later how the algorithm handles firing events.

For convenience, we shall divide the state space of an individual cell into three regions, based on the instantaneous value of V: Region I, with V ≤ V; Region II, with V < V ≤ V+; and Region III, with V > V+. In addition, we define Region IV to be that where the cell is in the refractory state (recall that a cell is in the refractory state at time t if Tm−1 ≥ tτR, where Tm−1 < t is the last spiking time of that cell). In each of these regions, the equation governing the dynamics for the gating variable nh is different, and this needs to be reflected in our algorithm.

At a given time t, each cell in the network is in precisely one of the regions outlined in the preceding paragraph. We store a global array that tracks which region each of the cells is in and use this to select the appropriate equation to integrate, based on those presented in (13). Over an interval Δt, the region of a subset of the cells may change as the local voltage variable, V, crosses switching manifolds. Note that the continuity of trajectories restricts which transitions between regions are possible (for example, a transition between Region I and Region III is not permitted). Our evolution operator φh must account for this.

We identify where switching events occur by comparing the region of all of the cells before and after the operator φΔt is applied. Where switching events have occurred, we locate the switching time for cell j by searching for roots of the transcendental equation Xj1(s)Vx=0, ti ≤ s < ti+1 using Newton’s method, where the superscript denotes that we only consider the V component of the state vector Xj and Vx is the switching manifold that has been crossed. We remark that the Fréchet derivatives for use in Newton’s method are obtainable in closed form, and so the application of Newton’s method does not rely on numerical finite difference approximations for these derivatives.

After the switching time, s, has been found, the state of the cell is updated via Xj(s) = φΔs(tiXj(ti)), where Δss − ti, according to the governing equation for that region. Following this, the regional variable is updated to reflect the fact that the cell has transitioned to a new region. Finally, the remainder of the update step is taken, according to Xj(ti+1) = φΔs(sXj(s)), where Δsti+1 − s, by integrating the governing equation for new region. It is possible that between times s and ti+1, the cell passes through other switching manifolds, in which case, the procedure for identifying switching times and updating the region variable is repeated as many times as necessary.

Cells in Region IV are in the refractory state, governed by (12) and (14), in which the voltage variable is held fixed at VVr for a total duration of length τR. To account for this, we introduce for each cell, a counter that tracks how much time a cell has spent in the refractory period. Following a firing event, the spiking cell enters the refractory state, whereupon this counter is reset to zero (and the cell’s region is changed to Region IV). The cell will remain in the refractory state until the counter reaches τR. At this time, the region variable for that cell is updated to the appropriate value, based on the location of Vr relative to V and V+. For computational efficiency, we can take advantage of the fact that V is unchanging in the refractory period by replacing the voltage dynamics with V˙=1, thus letting V implement our refractory time counter in Region IV. Thus, if the cell spikes at time Tm, we have V(Tm) = 0 upon entering Region IV. To ensure continuity of solutions upon exiting the refractory period, we must also ensure to set V(TmτR) = Vr.

Since cells are independent of one another outside of firing times, this above computation to update the state of the system can be performed in parallel across the GPU architecture in the absence of firing events. Where firing events occur, we must ensure that we correctly update the state of the network at the firing time to reflect the coupling in the network. As for the switching times, we determine firing times by searching for roots of the transcendental equation Xj1(s)Vth=0, ti ≤ s < ti+1 using Newton’s method. We then amalgamate all firing events between times ti and ti+1 and find the minimum firing time, t, of this list. The states of all cells are updated to this time from ti:

X(t) = φΔt(tiX(ti)), 

where Δtt − ti and the superscript denotes that this corresponds to the state of the network in the limit as t approaches t from the left. Following this, the reset conditions are applied across the entire domain. Our approach allows for simultaneous firing events to occur by also applying reset conditions for all events that occurred with the time interval [ttε), where ε is a small positive real number. After the reset conditions have been applied, the state of the network is updated again as

X(ti+1)=φΔtc(t,X+(t)),

where Δtc=ti+1t and the superscript denotes that the state is evaluated after the reset conditions have been applied. Note that during the completion of the timestep, we must again check to see if other cells have reached threshold, and apply reset conditions as necessary. At the completion of a timestep, the state of the network is saved and the state at the next time step is computed.

We note that a further improvement in computational efficiency can be achieved in the limit of high gain for the kernel function (7). Namely in the limit β → ∞ the bump function becomes a Top-hat function:

w(r)={w0,rσ,0,r>σ.
52

Thus, cells are only coupled with a single strength w0, over a finite distance, σ, and we need only consider the behaviour of other cells within this distance when applying reset conditions.

A.2 Implementation

To take full advantage of the parallel capabilities of GPUs, we must sensibly organise the processes associated with simulating the network through the construction of appropriate kernels to best to perform tasks in parallel. In addition to this, the memory associated with state variables must be managed appropriately given that memory access can present a significant bottleneck during GPU-based computations.

As discussed in the preceding section, the evolution operator φh has to account for the differing functional forms of the governing equations as V crosses thresholds. This is done by using Newton’s method to find the times of any crossings and then by ‘patching’ together solution orbits on either side of the threshold crossing time. Associated with Newton’s method is a tolerance for finding the crossing times, which must be specified along with other parameters. In our implementation, we take advantage of the double precision capability of our GPU. Since we are using closed-form solutions for the orbits of our system, the only source of numerical error other than machine error arises due to the tolerance selected when computing the threshold crossing times for V. In all of our simulations, we prescribe this tolerance to be 10−10 so that we have precise specification of these times.

During an update step, there may be multiple cells that reach the firing threshold at differing times. To preserve the correct dynamics in spite of this, we need to update the state of the network to the minimum of this set of firing times. To identify this, we define a variable that stores the minimum firing time across the network, which is updated by the subset of spiking cells as firing events are found. Without intervention, this procedure can lead to race conditions, in which one spiking cell will query the stored value of the minimum firing time (to compare with its own), whilst another spiking cell is replacing that same value. Race conditions mean that we cannot guarantee that the value stored in our firing time variable truly represents the minimum over all firing events. To address this problem, we use atomic events, which allow only one cell to read from and write to the stored firing time variable at once. In this way, we are guaranteed to find the true minimum across the set of firing times.

In Algorithms 1 and 2, we detail how to evolve the dynamics of the 2D network. Note that all of the for loops, and the use of φh to update the state of the network are performed in parallel using the GPU. The full code for our simulations, written in the C++/CUDA programming language, is available in the supplementary material.

Algorithm 1

An external file that holds a picture, illustration, etc.
Object name is 13408_2017_51_Figa_HTML.jpg

Evolution of 2D network

Algorithm 2

An external file that holds a picture, illustration, etc.
Object name is 13408_2017_51_Figb_HTML.jpg

EvolutionOperator

Appendix B: Perturbations at Switching Conditions

The model undergoes switching in its dynamics as the voltage variable passes through the values V+, V, and Vth. With the introduction of a set of indicator functions h(X(ξt); μ) = V(ξt)−Vμ, where μ ∈ { + , −, th} we can define the travelling wave coordinate values at which these switching events occur according to h(X(ξt); μ) = 0. Now suppose that we have two trajectories: an unperturbed trajectory X(ξt) = (V(ξt), nh(ξt)) and a perturbed trajectory X˜(ξ,t) such that δX(ξ,t)=X˜(ξ,t)X(ξ,t), with δX small. Moreover, let us consider the unperturbed trajectory to pass through the switching manifold when ξξm(t), m ∈ ℤ. Similarly we shall consider the perturbed trajectory to switch when ξ=ξ˜m(t)=ξm(t)+δξm(t). The indicator function for the perturbed trajectory may be Taylor expanded as (suppressing the dependence on μ for clarity):

h(X˜(ξ˜m,t))=h(X˜(ξm+δξm),t)=h(X(ξm+δξm,t)+δX(ξm+δξm,t))h(X(ξm,t)+Xξ(ξm,t)δξm)+Xh(X(ξm+δξm,t))δX(ξm+δξm,t)h(X(ξm,t))+Xh(X(ξm,t))Xξ(ξm,t)δξm+Xh(X(ξm,t))δX(ξm,t).
53

Here we have introduced the notation X(ξ±m,t)=limϵ0X(ξm±ϵ,t), to make sure that the partial derivative in ξ is well defined. Using the fact that h(X˜(ξ˜m,t))=0=h(X(ξm,t)) we obtain

Xh(X(ξm,t))[δX(ξm,t)+Xξ(ξm,t)δξm]=0.
54

Using the result that Xh(Xμ) = (∂V, ∂nh)(V − Vμ) = (1, 0) the above can be re-arranged to give the perturbation in the switching coordinate in terms of the difference between the perturbed and unperturbed trajectories as

δξm(t)=δV(ξm,t)Vξ(ξm,t).
55

Appendix C: Saltation Matrices

Saltation matrices allow us to handle any jumps in the system (or its linearisation) when it changes from one dynamical regime to another. As well as occurring at the switching manifolds this also happens when the voltage is unclamped and released from the refractory state. Using the notation of Appendix B let us first consider the deviation between the two trajectories at a switching event defined by ξξs, with a set of perturbed switching events at ξξsδξ, as

δX(ξs+δξ,t)=X˜(ξs+δξ,t)X(ξs+δξ,t)X˜(ξs,t)+X˜ξ(ξs,t)δξ[X(ξs,t)+Xξ(ξs,t)δξ]=δX(ξs,t)+[X˜ξ(ξs,t)Xξ(ξs,t)]δξ.
56

If δξ > 0 then the unperturbed trajectory will already have transitioned through the switch (from below), in which case the two trajectories are governed by different dynamics. A similar argument holds for δξ < 0. Thus we may write

δX(ξsδξt) ≃ δX(ξst) + [Xξ(ξst)−Xξ(ξs+t)]δξ.
57

Combining (55) and (57) gives

[δV(ξs+δξ,t)δnh(ξs+δξ,t)]=[δV(ξs,t)δnh(ξs,t)]δV(ξs,t)Vξ(ξs,t)[Vξ(ξs,t)Vξ(ξs+,t)nh,ξ(ξs,t)nh,ξ(ξs+,t)].
58

We may write the above in matrix form as

[δV(ξs+δξ,t)δnh(ξs+δξ,t)]=[1Vξ(ξs,t)Vξ(ξs+,t)Vξ(ξs,t)0nh,ξ(ξs,t)nh,ξ(ξs+,t)Vξ(ξs,t)1][δV(ξs,t)δnh(ξs,t)],
59

or equivalently as δX(ξsδξt) = K(ξs)δX(ξst) with

K(ξs)=[Vξ(ξs+,t)/Vξ(ξs,t)0(nh,ξ(ξs+,t)nh,ξ(ξs,t))/Vξ(ξs,t)1].
60

Now since the voltage is clamped immediately after a firing event (so that Vξ(ξs+t) = 0) and the dynamics for nh jumps (since it depends on V which is discontinuously reset from Vth to Vr) then the saltation matrix for firing is given by

Kfire(ξs)=[00(nh,ξ(ξs+,t)nh,ξ(ξs,t))/Vξ(ξs,t)1].
61

At a switching event whenever VV± we note that the voltage and gating dynamics are both continuous. Thus the saltation matrix for switching is given simply by Kswitch(ξs) = I2, namely there is no effect.

The use of saltation matrices to propagate perturbations through the refractory state is a little more subtle than through switching and firing events, since the former occur over a finite time-scale τR whilst the latter are instantaneous. In this case the perturbation δξ at a firing event ξξf is propagated for a time τR before a new dynamical regime is encountered. From (55) δξ = −δV(ξft)/Vξ(ξft). Setting ξsξfτR and combining the above with (57) gives

δX(ξs+δξ,t)δX(ξs,t)δV(ξf,t)Vξ(ξf,t)[Vξ(τR,t)Vξ(τR+,t)nh,ξ(τR,t)nh,ξ(τR+,t)].
62

Hence we see that δX(ξs+t) = δX(ξst) + Kref(ξs)δX(ξft) where

Kref(ξs)=[Vξ(ξs+)/Vξ(ξf)000].
63

Here we have used the fact that Vξ(ξst) = 0 (since the system is in its refractory state), and that the dynamics for nh is continuous at ξξs.

Footnotes

Electronic Supplementary Material

The online version of this article (doi:10.1186/s13408-017-0051-7) contains supplementary material.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

MBQ, KCAW, RDO and SC contributed equally. All authors read and approved the final manuscript.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Contributor Information

Mayte Bonilla-Quintana, ku.ca.mahgnitton@anatniuQallinoB.etyaM.

Kyle C. A. Wedgwood, ku.ca.retexe@doowgdeW.A.C.K.

Reuben D. O’Dea, ku.ca.mahgnitton@aed'o.nebuer.

Stephen Coombes, ku.ca.mahgnitton@sebmooc.nehpets.

References

1. Hasselmo ME. Neuronal rebound spiking, resonance frequency and theta cycle skipping may contribute to grid cell firing in medial entorhinal cortex. Philos Trans R Soc Lond B, Biol Sci. 2013;369 doi: 10.1098/rstb.2012.0523. [PMC free article] [PubMed] [Cross Ref]
2. Fyhn M, Molden S, Witter MP, Moser EI, Moser M-B. Spatial representation in the entorhinal cortex. Science. 2004;305:1258–1264. doi: 10.1126/science.1099901. [PubMed] [Cross Ref]
3. Tolman EC. Cognitive maps in rats and men. Psychol Rev. 1948;55:189–208. doi: 10.1037/h0061626. [PubMed] [Cross Ref]
4. O’Keefe J, Nadal L. The hippocampus as a cognitive map. Oxford: Clarendon; 1978.
5. Hafting T, Fyhn M, Molden S, Moser M-B, Moser EI. Microstructure of a spatial map in the entorhinal cortex. Nature. 2005;436:801–806. doi: 10.1038/nature03721. [PubMed] [Cross Ref]
6. Brun VH, Solstad T, Kjelstrup KB, Fyhn M, Witter MP, Moser EI, Moser M-B. Progressive increase in grid scale from dorsal to ventral medial entorhinal cortex. Hippocampus. 2008;18:1200–1212. doi: 10.1002/hipo.20504. [PubMed] [Cross Ref]
7. Rowland DC, Roudi Y, Moser M-B, Moser EI. Ten years of grid cells. Annu Rev Neurosci. 2016;39:19–40. doi: 10.1146/annurev-neuro-070815-013824. [PubMed] [Cross Ref]
8. Giocomo LM, Moser M-B, Moser EI. Computational models of grid cells. Neuron. 2011;71(4):589–603. doi: 10.1016/j.neuron.2011.07.023. [PubMed] [Cross Ref]
9. Schmidt-Hieber C, Häusser M. How to build a grid cell. Philos Trans R Soc Lond B, Biol Sci. 2014;369 doi: 10.1098/rstb.2012.0520. [PMC free article] [PubMed] [Cross Ref]
10. Burgess N, Barry C, O’Keefe J. An oscillatory interference model of grid cell firing. Hippocampus. 2007;17:801–812. doi: 10.1002/hipo.20327. [PMC free article] [PubMed] [Cross Ref]
11. Reifenstein ET, Kempter R, Schreiber S, Stemmler MB, Herz AVM. Grid cells in rat entorhinal cortex encode physical space with independent firing fields and phase precession at the single-trial level. Proc Natl Acad Sci USA. 2012;109:6301–6306. doi: 10.1073/pnas.1109599109. [PubMed] [Cross Ref]
12. Fuhs MC, Touretzky DS. A spin glass model of path integration in rat medial entorhinal cortex. J Neurosci. 2006;26:4266–4276. doi: 10.1523/JNEUROSCI.4353-05.2006. [PubMed] [Cross Ref]
13. Burak Y, Fiete IR. Accurate path integration in continuous attractor network models of grid cells. PLoS Comput Biol. 2009;5 doi: 10.1371/journal.pcbi.1000291. [PMC free article] [PubMed] [Cross Ref]
14. McNaughton BL, Battaglia FP, Jensen O, Moser EI, Moser M-B. Path integration and the neural basis of the ‘cognitive map’ Nat Rev Neurosci. 2006;7:663–678. doi: 10.1038/nrn1932. [PubMed] [Cross Ref]
15. Guanella A, Kiper D, Verschure P. A model of grid cells based on a twisted torus topology. Int J Neural Syst. 2007;17:231–240. doi: 10.1142/S0129065707001093. [PubMed] [Cross Ref]
16. Mhatre H, Gorchetchnikov A, Grossberg S. Grid cell hexagonal patterns formed by fast self-organized learning within entorhinal cortex. Hippocampus. 2012;22:320–334. doi: 10.1002/hipo.20901. [PubMed] [Cross Ref]
17. Pilly PK, Grossberg S. Spiking neurons in a hierarchical self-organizing map model can learn to develop spatial and temporal properties of entorhinal grid cells and hippocampal place cells. PLoS ONE. 2013;8(4) doi: 10.1371/journal.pone.0060599. [PMC free article] [PubMed] [Cross Ref]
18. Pilly PK, Grossberg S. How does the modular organization of entorhinal grid cells develop? Front Human Neurosci. 2014;8 doi: 10.3389/fnhum.2014.00337. [PMC free article] [PubMed] [Cross Ref]
19. Barry C, Hayman R, Burgess N, Jeffery KJ. Experience-dependent rescaling of entorhinal grids. Nat Neurosci. 2007;10:682–684. doi: 10.1038/nn1905. [PubMed] [Cross Ref]
20. Stensola H, Stensola T, Solstad T, Frøland K, Moser M-B, Moser EI. The entorhinal grid map is discretized. Nature. 2012;492:72–78. doi: 10.1038/nature11649. [PubMed] [Cross Ref]
21. Couey JJ, Witoelar A, Zhang S-J, Zheng K, Ye J, Dunn B, Czajkowski R, Moser M-B, Moser EI, Roudi Y, Witter MP. Recurrent inhibitory circuitry as a mechanism for grid formation. Nat Neurosci. 2013;16:318–324. doi: 10.1038/nn.3310. [PubMed] [Cross Ref]
22. Buetfering C, Allen K, Monyer H. Parvalbumin interneurons provide grid cell-driven recurrent inhibition in the medial entorhinal cortex. Nat Neurosci. 2014;17:710–718. doi: 10.1038/nn.3696. [PubMed] [Cross Ref]
23. Alonso A, Llinás RR. Subthreshold Na+-dependent theta-like rhythmicity in stellate cells of entorhinal cortex layer II. Nature. 1989;342:175–177. doi: 10.1038/342175a0. [PubMed] [Cross Ref]
24. Pastoll H, Ramsden HL, Nolan MF. Intrinsic electrophysiological properties of entorhinal cortex stellate cells and their contribution to grid cell firing fields. Front Neural Circuits. 2012;6 doi: 10.3389/fncir.2012.00017. [PMC free article] [PubMed] [Cross Ref]
25. Giocomo LM, Zilli EA, Fransén E, Hasselmo ME. Temporal frequency of subthreshold oscillations scales with entorhinal grid cell field spacing. Science. 2007;315:1719–1722. doi: 10.1126/science.1139207. [PMC free article] [PubMed] [Cross Ref]
26. Nolan MF, Dudman JT, Dodson PD, Santoro B. HCN1 channels control resting and active integrative properties of stellate cells from layer II of the entorhinal cortex. J Neurosci. 2007;27:12440–12451. doi: 10.1523/JNEUROSCI.2358-07.2007. [PubMed] [Cross Ref]
27. Navratilova Z, Giocomo LM, Fellous J-M, Hasselmo ME, McNaughton BL. Phase precession and variable spatial scaling in a periodic attractor map model of medial entorhinal grid cells with realistic after-spike dynamics. Hippocampus. 2012;22:772–789. doi: 10.1002/hipo.20939. [PubMed] [Cross Ref]
28. Hasselmo ME, Shay CF. Grid cell firing patterns may arise from feedback interaction between intrinsic rebound spiking and transverse traveling waves with multiple heading angles. Front Syst Neurosci. 2014;8 doi: 10.3389/fnsys.2014.00201. [PMC free article] [PubMed] [Cross Ref]
29. Tsuno Y, Chapman GW, Hasselmo ME. Rebound spiking properties of mouse medial entorhinal cortex neurons in vivo. Eur J Neurosci. 2015;42:2974–2984. doi: 10.1111/ejn.13097. [PMC free article] [PubMed] [Cross Ref]
30. Shay CF, Ferrante M, Chapman GW, Hasselmo ME. Rebound spiking in layer II medial entorhinal cortex stellate cells: possible mechanism of grid cell function. Neurobiol Learn Mem. 2015;129:83–98. doi: 10.1016/j.nlm.2015.09.004. [PMC free article] [PubMed] [Cross Ref]
31. Ferrante M, Shay CF, Tsuno Y, Chapman GW, Hasselmo ME. Post-inhibitory rebound spikes in rat medial entorhinal layer II/III principal cells: in vivo, in vitro, and computational modeling characterization. Cereb Cortex. 2017;27(3):2111–2125. [PubMed]
32. Heys JG, Hasselmo ME. Neuromodulation of IhIh in layer II medial entorhinal cortex stellate cells: a voltage clamp study. J Neurosci. 2012;32:9066–9072. doi: 10.1523/JNEUROSCI.0868-12.2012. [PMC free article] [PubMed] [Cross Ref]
33. Shay CF, Boardman IS, James NM, Hasselmo ME. Voltage dependence of subthreshold resonance frequency in layer II of medial entorhinal cortex. Hippocampus. 2012;22:1733–1749. doi: 10.1002/hipo.22008. [PMC free article] [PubMed] [Cross Ref]
34. Giocomo LM, Hasselmo ME. Time constants of h current in layer II stellate cells differ along the dorsal to ventral axis of medial entorhinal cortex. J Neurosci. 2008;28:9414–9425. doi: 10.1523/JNEUROSCI.3196-08.2008. [PMC free article] [PubMed] [Cross Ref]
35. Raudies F, Brandon MP, Chapman GW, Hasselmo ME. Head direction is coded more strongly than movement direction in a population of entorhinal neurons. Brain Res. 2015;1621:355–367. doi: 10.1016/j.brainres.2014.10.053. [PMC free article] [PubMed] [Cross Ref]
36. Rinzel J, Terman D, Wang X-J, Ermentrout B. Propagating activity patterns in large-scale inhibitory neuronal networks. Science. 1998;279:1351–1355. doi: 10.1126/science.279.5355.1351. [PubMed] [Cross Ref]
37. Coombes S. Dynamics of synaptically coupled integrate-and-fire-or-burst neurons. Phys Rev E. 2003;67 doi: 10.1103/PhysRevE.67.041910. [PubMed] [Cross Ref]
38. Wasylenko TM, Cisternas JE, Laing CR, Kevrekidis IG. Bifurcations of lurching waves in a thalamic neuronal network. Biol Cybern. 2010;103:447–462. doi: 10.1007/s00422-010-0409-3. [PubMed] [Cross Ref]
39. Evans J. Nerve axon equations: IV the stable and unstable impulse. Indiana Univ Math J. 1975;24:1169–1190. doi: 10.1512/iumj.1975.24.24096. [Cross Ref]
40. Coombes S, Owen MR. Evans functions for integral neural field equations with Heaviside firing rate function. SIAM J Appl Dyn Syst. 2004;34:574–600. doi: 10.1137/040605953. [Cross Ref]
41. Aizerman MA, Gantmakher FR. On the stability of periodic motions. J Appl Math Mech. 1958;22:1065–1078. doi: 10.1016/0021-8928(58)90033-9. [Cross Ref]
42. Leine RI, Van Campen DH, Van de Vrande BL. Bifurcations in nonlinear discontinuous systems. Nonlinear Dyn. 2000;23:105–164. doi: 10.1023/A:1008384928636. [Cross Ref]
43. di Bernardo M, Budd C, Champneys AR, Kowalczyk P. Piecewise-smooth dynamical systems: theory and applications. London: Springer; 2008.
44. Coombes S, Thul R, Wedgwood KCA. Nonsmooth dynamics in spiking neuron models. Physica D. 2012;241:2042–2057. doi: 10.1016/j.physd.2011.05.012. [Cross Ref]
45. Sreenivasan S, Fiete I. Grid cells generate an analog error-correcting code for singularly precise neural computation. Nat Neurosci. 2011;14:1330–1337. doi: 10.1038/nn.2901. [PubMed] [Cross Ref]
46. Burgalossi A, Herfst L, von Heimendahl M, Förste H, Haskic K, Schmidt M, Brecht M. Microcircuits of functionally identified neurons in the rat medial entorhinal cortex. Neuron. 2011;70:773–786. doi: 10.1016/j.neuron.2011.04.003. [PubMed] [Cross Ref]

Articles from Journal of Mathematical Neuroscience are provided here courtesy of Springer-Verlag