A neural mass model of place cell activity: theta phase precession, replay and imagination of never experienced paths

Abstract

Recent results on hippocampal place cells show that the replay of behavioral sequences does not simply reflect previously experienced trajectories, but may also occur in the reverse direction, or may even include never experienced paths. In order to elucidate the possible mechanisms at the basis of this phenomenon, we have developed a model of sequence learning. The present model consists of two layers of place cell units. Long-range connections among units implement heteroassociation between the two layers, trained with a temporal Hebb rule. The network was trained assuming that a virtual rat moves within a virtual maze. This training leads to the formation of bidirectional synapses between the two layers, i.e. synapses connecting a neuron both with its previous and subsequent element in the path. Subsequently, two distinct conditions were simulated with the trained network. During an exploratory phase, characterized by a similar consideration to the external environment and to the internal representation, the model simulates the occurrence of theta precession in the forward path and the temporal compression. During an imagination phase, when there is no consideration to the external location, the model produces trains of gamma oscillations, without the presence of a theta rhythm, and simulates the occurrence of both direct and reverse replay, and the imagination of never experienced paths. The new paths are built by combining bunches of previous trajectories. The main mechanisms at the basis of this behavior are explained in detail, and lines for future improvements (e.g., to simulate preplay) are discussed.

Appendices

Appendix A: Equations of a single unit

The model of a single cortical column consists of four neural populations (Fig. 1), which represent pyramidal neurons (subscript p), excitatory interneurons (subscript e), and inhibitory interneurons with slow and fast synaptic kinetics (GABA_A, slow and GABA_A, fast, subscripts s and f, respectively). All populations are described with a similar mathematical formalism. Briefly, each population computes an average post-synaptic membrane potential (v, in (A4), (A8), (A12) and (A16)) from other neural populations and from external inputs, and converts this membrane potential into an average density of spikes fired by the neurons (z, in (A3), (A7), (A11) and (A15)). In order to account for the presence of inhibition (when potential is below a given threshold) and saturation (when potential is high) this conversion is simulated with a static sigmoidal relationship.

To model a whole cortical column, the four populations are connected via excitatory and inhibitory synapses, with impulse response h _e (t), h _s (t) or h _f (t), assuming that synapses from pyramidal neurons and excitatory interneurons have similar dynamics. The kinetics of each synapse are described with a second-order system, but with different parameter values. By denoting with z the input to a synapse, and with y the synapse output, its dynamics can be described through a second order differential equation (equivalent to two first-order differential equations, i.e., Eqs (A1)-(A2), (A5)-(A6), (A9)-(A10) and (A13)-(A14) for the four populations.

This corresponds to the following equations

Pyramidal neurons

$$ \frac{d{y}_p(t)}{dt}={x}_p(t) $$

(A1)

$$ \frac{d{x}_p(t)}{dt}=\frac{G_e}{\tau_e}{z}_p(t)-\frac{2}{\tau_e}{x}_p(t)-\frac{y_p(t)}{\tau_e^2} $$

(A2)

$$ {z}_p(t)=\frac{2{e}_o}{1+{e}^{r\left( so-{v}_p\right)}} $$

(A3)

$$ {v}_p(t)={C}_{pe}{y}_e(t)-{C}_{ps}{y}_s(t)-{C}_{pf}{y}_f(t)+E(t) $$

(A4)

Excitatory interneurons

$$ \frac{d{y}_e(t)}{dt}={x}_e(t) $$

(A5)

$$ \frac{d{x}_e(t)}{dt}=\frac{G_e}{\tau_e}\left({z}_e(t)+\frac{u_p(t)}{C_{pe}}\right)-\frac{2}{\tau_e}{x}_e(t)-\frac{y_e(t)}{\tau_e^2} $$

(A6)

$$ {z}_e(t)=\frac{2{e}_o}{1+{e}^{r\left({s}_o-{v}_e\right)}} $$

(A7)

$$ {v}_e(t)={C}_{ep}{y}_p(t) $$

(A8)

Slow inhibitory interneurons

$$ \frac{d{y}_s(t)}{dt}={x}_s(t) $$

(A9)

$$ \frac{d{x}_s(t)}{dt}=\frac{G_s}{\tau_s}{z}_s(t)-\frac{2}{\tau_s}{x}_s(t)-\frac{y_s(t)}{\tau_s^2} $$

(A10)

$$ {z}_s(t)=\frac{2{e}_o}{1+{e}^{r\left({s}_o-{v}_s\right)}} $$

(A11)

$$ {v}_s(t)={C}_{sp}{y}_p(t) $$

(A12)

Fast inhibitory interneurons

$$ \frac{d{y}_f(t)}{dt}={x}_f(t) $$

(A13)

$$ \frac{d{x}_f(t)}{dt}=\frac{G_f}{\tau_f}{z}_f(t)-\frac{2}{\tau_f}{x}_f(t)-\frac{y_f(t)}{\tau_f^2} $$

(A14)

$$ {z}_f(t)=\frac{2{e}_o}{1+{e}^{r\left({s}_o-{v}_f\right)}} $$

(A15)

$$ {v}_f(t)={C}_{fp}{y}_p(t)-{C}_{fs}{y}_s(t)-{C}_{ff}{y}_f(t)-{y}_l(t)+I(t) $$

(A16)

$$ \frac{d{y}_l(t)}{dt}={x}_l(t) $$

(A17)

$$ \frac{d{x}_l(t)}{dt}=\frac{G_e}{\tau_e}{u}_f(t)-\frac{2}{\tau_e}{x}_l(t)-\frac{y_l(t)}{\tau_e^2} $$

(A18)

In the previous equations, the term E(t) (Eq. (A4)) represents the overall input entering into the pyramidal neurons of that unit, from units in the other layer (inter-layer connections). The term I(t) (Eq. (A16)) represents the overall input entering into GABA_A,fast interneurons of the unit from pyramidal neurons of other units in the same layer (intra-layer lateral connections). The terms u _p (t) (Eq. (A6)) and u _f (t) (Eq. (A18)) represent exogenous inputs. They include Gaussian noise and (when available) information on the present position of the rat (see Appendix C).

To improve model reliability, we assumed that the two inputs (u _p (t) and u _f (t)) affect the target neurons through the typical dynamics of excitatory synapses. This is justified by the idea that all long range connections in the brain (including exogenous inputs) come from pyramidal neurons. To this end, the quantity u _p (t) is included in Eq. (A6), i.e., we exploit the same dynamics of an excitatory interneuron to affect the target pyramidal neurons. This truck avoids the introduction of two further state variables. Conversely, the same truck was not possible for what concerns the input quantity u _f (t) affecting the fast inhibitory interneurons (in fact, these interneurons do not receive any specific excitation). Hence, we introduced two additional state variables (Eqs. (A17) and (A18)) to describe the synaptic kinetics of the u _f (t) input.

For the meaning of other quantities, see Tables 1 and 2.

Table 1 Parameters of units

Table 2 Variables of units

In the subsequent paragraph, quantities belonging to the first and second layers will be denoted through the superscripts L₁ and L₂, respectively. The position within the layer will be denoted with a subscript (i or j).

Appendix B: The connections between units

Two types of long-range connections are used, both coming from pyramidal neurons of the pre-synaptic unit. Connections of type W affect the membrane potential v ^Lx _p,i of the pyramidal neurons in the target unit by an additive term E ^Lx _i (t). Connections of type A, instead, affect the membrane potential v ^Lx _f,i of the GABA_A,fast interneurons in the target unit by an additive term I ^Lx _i (t) (see Eqs. (A4) and (A16) and also Cona et al. (2012) for more mathematical details). In the following, two superscripts and two subscripts are used for each connection. The superscripts represent the layers to which the postsynaptic and the presynaptic units respectively belong. The two subscripts represent the positions of the postsynaptic and of the presynaptic units within these layers.

The expressions of E ^Lx _i (t) and I ^Lx _i (t) for the two layers are:

$$ {E}_i^{L1}(t)={\sum}_j{W}_{ij}^{L1L2}\cdot {y}_{p,j}^{L2}(t) $$

(B1)

$$ {I}_i^{L1}(t)={\sum}_{j\ne i}{A}_{ij}^{L1L1}\cdot {z}_{p,j}^{L1}(t) $$

(B2)

$$ {E}_i^{L2}(t)={W}_{ii}^{L2L1}\cdot {y}_{p,i}^{L1}(t) $$

(B3)

$$ {I}_i^{L2}(t)={\sum}_{j\ne i}{A}_{ij}^{L2L2}\cdot {z}_{p,j}^{L2}(t) $$

(B4)

The previous equations can be explained as follows. Eq. (B1) signifies that pyramidal neurons in L₁ receive excitatory inputs from pyramidal neurons in the L₂. These connections are learnt during the maze exploration (see below) and implement a heteroassociative network that enables units in L₂ to activate units in L₁ encoding for near areas of the maze. According to Eq. (B2), GABA_A,fast interneurons in L₁ receive inputs through synapses with fast (AMPA) kinetics from the pyramidal neurons of other units in the same layer. Eq. (B3) means that pyramidal neurons in L₂ receive an input from pyramidal neurons in L₁ located at the same position. Eq. (B4) has the same meaning as Eq. (B2), but with reference to L₂.

These expressions describe how the firing rates in source populations are mixed with different connection strengths to provide depolarizing or hyperpolarizing (excitatory or inhibitory) inputs to target populations. For what concerns the connections W, they are mediated by glutamatergic synapses, hence the postsynaptic potential y _p (t) was used as driving quantity (Eqs. (B1) and (B3)). Conversely, in order to characterize the faster AMPA synapses that mediate connections of type A in Eqs. (B2) and (B4), we used directly the firing rate z _p (t) of the presynaptic neuron, thus simulating a synapse with a negligible time constant (Cona and Ursino 2013; Cona et al. 2012).

Appendix C: The external input

As shown in Eqs (A6) and (A18) above, both pyramidal neurons and GABA_A,fast interneurons in all the units in both layers receive an additional external input. In particular, all units receive a white Gaussian noise input with zero mean value, as in previous works (Cona and Ursino 2013; Cona et al. 2012). Moreover, pyramidal neurons in layer L₁ also receive a sensory input, (say S _i(t)) which indicates where the rat is located at a given time.

Hence, the quantities u _p (t) and u _f (t) have the following expressions:

$$ {u}_{p,i}^{L1}(t)={S}_i(t)+{n}_{p,i}^{L1}(t) $$

(C1)

$$ {u}_{p,i}^{L2}(t)={n}_{p,i}^{L2}(t) $$

(C2)

$$ {u}_{f,i}^{L1}(t)={n}_{f,i}^{L1}(t) $$

(C3)

$$ {u}_{f,i}^{L2}(t)={n}_{f,i}^{L2}(t) $$

(C4)

where the quantities n(t) represent independent Gaussian noises with zero mean and assigned variance (see Table 3 for the values used during the present simulations).

Table 3 Parameters involving maze exploration and inputs

The input S _i (t) depends on the simulation phase (training, exploration or imagination) as described below.

Training phase - While the rat moves along the maze during the training phase, the layer L₁ receives a sensory stimulus ST(x, y; t) that indicates where the rat is located at a given time, and has therefore the shape of a bi-dimensional impulse:

$$ ST\left(x,y;t\right)=\delta \left(x-{r}_x(t),y-{r}_y(t)\right) $$

(C4)

where r _x (t) and r _y (t) are the rat coordinates at time t, and δ(x, y) represents a bidimensional Dirac delta function.

Each unit in L₁ responds to this stimulus according to a Gaussian shaped receptive field:

$$ R{F}_i\left(x,y\right)={G}_{RF}\cdot exp\left[-\frac{{\left(x-R{F}_{x,i}\right)}^2+{\left(y-R{F}_{y,i}\right)}^2}{2{\sigma}_{RF}^2}\right] $$

(C5)

where RF _x,i and RF _y,i are the central coordinates of the receptive field of unit i, while G _RF and σ _RF are the receptive field’s gain and standard deviation. RF _x,i and RF _y,i for the units in L₁ have been set to fully cover the maze with a rectangular grid (Fig. 2).

The input S _i (t) during training is then computed by calculating the convolution between the sensory input and the receptive field:

$$ {S}_i(t)=\int \int ST\left(x,y;t\right)\cdot R{F}_i\left(x,y\right) dxdy={G}_{RF}\cdot exp\left[-\frac{{\left({r}_x(t)-R{F}_{x,i}\right)}^2+{\left({r}_y(t)-R{F}_{y,i}\right)}^2}{2{\sigma}_{RF}^2}\right] $$

(C6)

which is sufficient to activate a few neurons that encode for areas around the point in the maze where the rat is at time t, while not affecting neurons that encode for distant areas.

Exploration phase - As described in the text, to mimic the exploration phase, when the rat is familiar with the environment, we assumed a decrease in the consideration toward the external stimulus, and an increased impact of the internal representation (i.e., on the activity of pyramidal neurons in L₁). This was simulated by reducing the gain G _RF to half the value used during the training phase, and adding an offset term, $\frac{1}{2}{S}_{offset}$, to all units. It is worth noting that the presence of this offset is the same as to assume a reduced threshold for all neurons. Moreover, when exploring a known environment, the rat must choose a given direction of movement to imagine a future path. For instance, this may correspond to head direction. To this end, we included a further receptive field obtained as the directional derivative of RF _i (x, y) along the direction in which the rat is moving $\widehat{v}(t)$. Thus, the input in this case becomes

$$ {S}_i(t)=\frac{1}{2}\int \int ST\left(x,y;t\right)\cdot \left[R{F}_i\left(x,y\right)-{D_{\hat{v}}}_{(t)}R{F}_i\left(x,y\right)\right] dxdy+\frac{1}{2}{S}_{offset} $$

(C7)

$$ \hat{v}(t)=\frac{\left[\frac{\partial {r}_x(t)}{\partial t},\frac{\partial {r}_y(t)}{\partial t}\right]}{\sqrt{\left(\frac{\partial {r}_x{(t)}^2}{\partial t}+\frac{\partial {r}_y{(t)}^2}{\partial t}\right)}} $$

(C8)

where −${D}_{\widehat{v}(t)}R{F}_i\left(x,y\right)$ gives a negative contribution to the units that encode for the maze locations visited in the near past and thus warrants that the sequence is evoked in the same direction in which the rat is moving.

Imagination phase - During the imagination phase, we assumed that the rat completely neglects the external stimulus, but it pays stronger consideration to its internal representation of the maze, to spontaneously evoke neuronal activity. This has been simulated by doubling the offset term, compared with that used in the exploration phase. We have

$$ {S}_i(t)={S}_{offset} $$

(C9)

The increase of this offset term is the same as to assume a moderate leftward shift in the sigmoidal relationship of the pyramidal populations, i.e., all the neurons are more excitable.

The values of all of the parameters are found in Table 3.

Appendix D: The Hebbian learning

While the virtual rat moves through the virtual maze its units in L₁ get activated by the inputs given in Eq. (C1). This activation is transferred to L₂ through the W ^L2L1 _ii connections, according to Eq. (B3). When units in the two layers are activated in close time synchronization, Hebbian learning between the units is possible. In this work Hebbian learning is implemented to train excitatory connections from L₂ to L₁ according to the following equations:

$$ {\tau}_{train}\frac{d{W}_{ij}^{L1L2}}{dt}={\gamma}_{ij}{\left[{z}_{p,i}^{L1}-\theta \right]}^{+}{\left[{m}_{p,j}^{L2}-\theta \right]}^{+} $$

(D1)

where γ _ij is a learning rate (depending on the particular synapse value), z ^L1 _p,i (t) is the spike density of the pyramidal population in the postsynaptic unit (in L₁), m ^L2 _p,j (t) is the average activity of pyramidal populations in the presynaptic unit (in L₂) computed in a previous temporal window consisting of N _s samples. Hence

$$ {m}_{p,j}^{L2}(t)=\frac{{\displaystyle {\sum}_{i=0}^{Ns-1}{z}_{p,j}^{L2}}\left(t-l\cdot Ts\right)}{N_s} $$

(D2)

It is worth noting that the right hand member of Eq. (D1) can be rewritten introducing the synchronization between the presynaptic and postsynaptic units, defined as follows

$$ sy{n}_{ij}^{L1L2}(t)=\frac{{\left[\frac{z_{p,i}^{L1}}{2{e}_o}-{\theta}_L\right]}^{+}{\left[\frac{m_{p,j}^{L2}}{2{e}_o}-{\theta}_L\right]}^{+}}{{\left(1-{\theta}_L\right)}^2} $$

(D3)

where 2e _o is the maximal activity of a population (see Eqs. (A3), (A7), (A11) and (A15)). Eq. (D3) means that if the activity of the postsynaptic unit at time t and the average activity of the presynaptic unit in the near past are both above a given threshold, then the synchronization term gets a high value. The denominator is a normalizing factor to have values in the [0;1] range. Comparing Eqs. (D1) and (D3), we have

$$ {\tau}_{train}\frac{d{W}_{ij}^{L1L2}}{dt}=\left({W}_{\max }-{W}_{ij}^{L1L2}\right)sy{n}_{ij}^{L1L2}(t) $$

(D4)

with θ = 2e _o θ _L and ${\gamma}_{ij}=\frac{W_{\max }-{W}_{ij}^{L1L2}}{{\left(2{e}_o-\theta \right)}^2}$.

In Eq. (D4) we assumed that the learning rate γ _ij decreases progressively as the connection reaches its maximum value. Finally, at each step, once the values of W ^L1L2 _ij for all the i,j pairs are updated, these values are normalized so that the sum of all connections entering each unit i has a fixed value W _norm :

$$ {W}_{ij}^{L1L2}=\frac{W_{ij}^{L1L2}}{{\displaystyle {\sum}_j{W}_{ij}^{L1L2}}}{W}_{norm} $$

(D5)

Appendix E: Parameters assignment

The values of most parameters of the units (Table 1) are identical to those used in our previous works (Cona and Ursino 2013; Cona et al. 2012). In more detail, the only parameters of the units that changed in this work are the connection strengths from GABA_A,slow interneurons to pyramidal neurons C _ps , and from pyramidal neurons to GABA_A,fast interneurons C _fp . Both these changes were necessary to account for differences in the network architecture as commented in section Discussion.

The parameters concerning the maze, its exploration and the corresponding inputs (Table 3), such as the dimension of the receptive fields or the rat velocity, were mostly assigned following indications from real values measured in experimental settings (Dragoi and Buzsáki 2006). The value of the sensory input (G _RF ) during training was assigned to have a significant activation of units in L₁ as soon as the input enters the periphery of the receptive field. The value of neural activation during the imagination phase (S _offset ) was chosen so that a few neurons in each simulation can be activated randomly by the internal noise, thus starting an imagined sequence; this value must not be excessively high, to allow inhibition from slow GABAergic synapses to induce a refractory period which avoids the presence of reflected activation. The values of the inputs used for the exploration of the familiar maze were arbitrarily assigned as G _RF /2 and S _offset /2 respectively, so that exploration of the familiar maze can be considered an average behavior between the exploration of the unknown maze (when the rat is just paying consideration to the external world) and the imagination of the known maze (when neuron activity completely depends on the internal dynamics). This simple choice provided satisfactory results that agree with what was expected from this modality.

The values of the long range connections (Table 4) were mostly inherited from our previous work (Cona and Ursino 2013). It should be noted that these connection strengths have been multiplied by a factor of 36, since in the previous work each unit received connections from neural assemblies consisting of 36 units, while in the present work each neural assembly consists of a single unit.

Table 4 Parameters of long range connections

Further slight adjustments were then required to compensate for the differences in the connectivity pattern of the present network, such as the absence of the inter-layer connections between units encoding for different features of the same object and the connection matrix W ^L1L2 _ij that is now symmetric.