A robotic model of hippocampal reverse replay for reinforcement learning

The network model of place cells represents a simplified hippocampal CA3 network, where CA3 stands for the hippocampal Cornu Ammonis 3 region, capable of generating reverse replays of recent place cell sequence trajectories. We presented this model of reverse replays in [64], but with one minor modification. Whereas the reverse replay model in [64] has a global inhibitory term acting on all place cells, in the present version, the place cells have those inhibitory inputs removed from their dynamics. Instead, our model uses a binary parameter to control synaptic activity, which functions similar to inhibition, see equation (5) below; therefore, the inhibitory inputs are not necessary. This modification does not affect the ability of the network to produce reverse replays (see supplementary material), where we compare reverse replays both with and without global inhibition.

In more detail, the place cells consist of a network of 100 neurons, each of which is bidirectionally connected to its eight nearest neighbours as determined by the positioning of their place fields. Hence, place cells with neighbouring place fields are bidirectionally connected (figure 1(B)), whereas place cells whose place fields are further than one place field apart are not. In this manner, the network's connectivity represents a map of the environment. This network approach is similar to the one taken by Haga and Fukai [22] in their model of reverse replay, except their weights are plastic whilst we keep ours static. Here, we keep only essential mechanisms to study the interplay between reverse replay and RL. The static weights for each cell, represented by $w_{jk}^{\,place}$ indicating the weight projecting from neuron k onto neuron j, are all set to 1, with no cells self-projecting to themselves. Figure 1(B) displays the full connectivity schema for the bidirectionally connected place cell network.

The rate for each place cell neuron, represented by x_j , is given as a linearly rectified rate with upper and lower bounds,

$$\begin{align} x_j = \begin{cases} 0 &\mathrm{if}\ x_j^{^{\prime}} \lt 0\\ 100 &\mathrm{if}\ x_j^{^{\prime}} \gt 100\\ x_j^{^{\prime}} &\mathrm{otherwise} \end{cases}. \end{align} \tag{ 3 }$$

The variable $x_j^{^{\prime}}$ is defined as,

$$\begin{align*} x_j^{^{\prime}} = \alpha \left(I_j - \epsilon \right), \end{align*}$$

where α and ε are constants determining the scaling factor and threshold of the linear rectifier, respectively. I_j is the cell's activity, which evolves according to time-decaying first-order dynamics,

$$\begin{align} \unicode{x03C4}_I\frac{d}{dt}I_j = -I_j + \psi_jI_j^{syn} + I_j^{\,place}, \end{align} \tag{ 4 }$$

where $\unicode{x03C4}_I$ is the time constant, $I_j^{syn}$ is the synaptic input from the cell's neighbouring neurons, and $I_j^{\,place}$ is the place-specific input calculated as per a normal distribution of MiRo's position from the place field's centre point. ψ_j represents the place cell's intrinsic plasticity, detailed further below.

Each place cell has been associated with a field in the environment defined by its centre point and width, with place fields distributed evenly across the environment (100 in total). As stated, the place-specific input, $I_j^{\,place}$, is computed from a two-dimensional normal distribution determined by MiRo's distance from the place field's centre point,

$$\begin{align} I^{\,place}_j = I^{\,p}_{max} \mathrm{exp}\left[- \frac{(x_{MiRo}^c - x_j^c)^2 + (y_{MiRo}^c - y_j^c)^2}{2d^2} \right], \end{align} \tag{ 5 }$$

where $I^{\,p}_{max}$ determines the max value for the place cell input. The coordinates $(x_{MiRo}^c, y_{MiRo}^c)$ represent MiRo's (x, y) position in the environment, whilst $(x_j^c, y_j^c)$ is the location of the place field's centre point. The term d in the denominator is a constant determining the width of the place field. For simplicity, we do not model the formation of the place cells from the visual input. We assume the robot coordinates are known, hence the place-cell activity defined by equation (5). From a machine learning point of view, this equation converts a low dimensional representation (coordinates) to a high dimensional representation (place cells activity).

The synaptic inputs, $I_j^{syn}$, are computed as a sum over neighbouring synaptic inputs modulated by the effects of short-term depression and facilitation, D_k and F_k , respectively,

$$\begin{align} I^{syn}_j = \lambda\sum_{k = 1}^8 w_{jk}^{\,place} x_k D_k F_k, \end{align} \tag{ 6 }$$

where $w_{jk}^{\,place}$ is the weight projecting from place cell k onto place cell j. In this model, all these weights are fixed at a value of 1. Parameter λ takes on a value of 0 or 1 depending on whether MiRo is exploring (λ = 0) or is at the reward (λ = 1). It prevents synaptic transmissions during exploration but not whilst MiRo is at the reward (the point at which reverse replays occur). Therefore, while the robot moves in the environment, λ = 0 and thus $I^{syn}_j = 0$. When it receives a reward, and during reverse replays, $I^{syn}_j$ has a non-zero value. A similar two-stage approach can be found in other models as a means to separate an encoding stage during exploration from a retrieval stage [52], and was a key feature of some of the early associative memory models [25]. Experimental evidence also supports this two-stage process due to the effects of acetylcholine. Acetylcholine levels are high during exploration but drop during rest [29]. Acetylcholine suppresses the recurrent synaptic transmissions in the hippocampal CA3 region [24]. We do not explicitly model this process. Instead, we consider the λ parameter conceptually corresponding to high acetylcholine levels (λ = 0) and low acetylcholine levels (λ = 1). We want to underline that the global inhibitory inputs found in our earlier work [64] were unnecessary. The λ term effectively plays the role of inhibitory inputs (inhibition is decreased during reverse replays, thus increasing synaptic transmission), yet is simpler to implement.

D_k and F_k in equation (6) are respectively the short-term depression and short-term facilitation terms, and for each place cell these are computed as (as in [22], but see [11, 58, 60, 61]),

$$\begin{align} \frac{d}{dt}D_k = \frac{1-D_k}{\unicode{x03C4}_{STD}} - x_k D_k F_k, \end{align} \tag{ 7 }$$

$$\begin{align} \frac{d}{dt}F_k = \frac{U-F_k}{\unicode{x03C4}_{STF}} + U \left( 1-F_k \right) x_k, \end{align} \tag{ 8 }$$

where $\unicode{x03C4}_{STD}$ and $\unicode{x03C4}_{STF}$ are the time constants, and U is a constant representing the steady-state value for short-term facilitation when there is no neuron activity $(x_k = 0)$. D_k and F_k each take on values in the range $[0,1]$. Notice that when $x_k \gt 0$, short-term depression is driven steadily towards 0, whereas short-term facilitation is driven steadily upwards towards 1. Modifying the time constants allows short-term depression or short-term facilitation effects to dominate. In this model, the time constants are chosen so that depression is the primary short-term effect. Our choice ensures that activity propagating from one neuron to the next during reverse replay events dissipates quickly, allowing for stable replays without activity exploding in the network. We note that while the reverse replay model has been adopted from [22], the parameters used here (table 1) differ. They are higher than typical values, albeit [60, 61] uses short-term plasticity parameters fitted to biological data that led to time constants as high as 900 ms.

Table 1. Summarising the model parameter values for the hippocampal network used in the experiments. All these parameters are kept constant across all experiments.

Parameter	Value
α	1 C−1
ε	2 A
$\unicode{x03C4}_I$	0.05 s
$I^p_{max}$	50 A
d	0.1 m
λ	0 or 1, see text
$\unicode{x03C4}_{STD}$	1.5 s
$\unicode{x03C4}_{STF}$	1 s
U	0.6
ψss	0.1
ψmax	4
$\unicode{x03C4}_{\psi}$	10 s
β	1
xψ	10 Hz

We turn finally to the intrinsic plasticity term in equation (4), represented by ψ_j . As observed in equation (4), its behaviour is to scale all incoming synaptic inputs. In [42], Pang and Fairhall used a heuristically developed sigmoid whose output was a function of the neuron's rate. Intrinsic plasticity in their model did not decay after its activation. Since our robot often travels across most of the environment, we needed a time-decaying form of intrinsic plasticity to avoid potentiating all cells in the network. The simplest form of such time-decaying intrinsic plasticity is, therefore,

$$\begin{align} \frac{d}{dt}\psi_j = \frac{\psi_{ss} - \psi_j}{\unicode{x03C4}_{\psi}} + \frac{\psi_{max} - 1}{1 + \mathrm{exp} \left[ -\beta (x_j - x_\psi) \right]}, \end{align} \tag{ 9 }$$

with again, $\unicode{x03C4}_{\psi}$ being its time constant, and ψ_ss being a constant that determines the steady state value for when the sigmoidal term on the right is 0. All of ψ_max , β and x_ψ are constants that determine the shape of the sigmoid. Since ψ_j could potentially grow beyond the value of ψ_max , we restrict ψ_j so that if $\psi_j \gt \psi_{max}$, then ψ_j is set to ψ_max .

To initiate a replay event, place cell inputs, computed using equation (5) with MiRo's current location at the reward, are input into the place cell dynamics (see equation (4)) one second after MiRo reaches the reward, for a duration of 100 ms. Intrinsic plasticity for those most recently active cells during the trajectory is increased, whilst synaptic conductance in the place cell network is turned on by setting λ = 1. Therefore, the place cell input activates only its adjacent cells that were recently active. This effect continues throughout all recently active cells, thus resulting in a reverse replay. Short-term depression ensures that the activity dissipates fast as it propagates from one neuron to the next.

The action cell values determine how MiRo moves in the environment. All place cells project feedforward through a set of plastic synapses to all action cells, as shown in figure 1(B). There are 72 action cells, the value of each drawn from a Gaussian distribution with mean $\tilde{y}_i$ and variance σ²,

$$\begin{align} y_i \sim \mathcal{N}\left( \tilde{y}_i, \sigma^2 \right). \end{align} \tag{ 10 }$$

The mean value $\tilde{y}_i$ is calculated as follows,

$$\begin{align} \tilde{y}_i = \frac{1}{1 + \mathrm{exp} \left[ -c_1 \sum_{j = 1}^{100} w_{ij}^{\,PC\textrm{-}AC} x_j - c_2 \right]}, \end{align} \tag{ 11 }$$

with c₁ and c₂ determining the shape of the sigmoid. $w_{ij}^{\,PC\textrm{-}AC}$ represents the weight projecting from place cell j onto action cell i. The sigmoidal function is one possible choice which results in saturating terms in the RL learning rule (appendix section 'Mathematical derivation of the place-action cell synaptic learning rule'); an alternative option, for instance, could have been a linear function. The action cells are restricted to take values between 0 and 1, i.e. $y_i \rightarrow [0,1]$, and be interpretable as normalised firing rates.

MiRo moves at a constant forward velocity, whereas the output of the action cells sets a target heading for MiRo to move in. This target heading is allocentric in that the heading is relative to the arena. The activity for each active cell is y_i and the target heading θ_target . We compute the population vector of the action cell values to find the heading from the cells:

$$\begin{align} \theta_{target} = \arctan \left( \frac{\sum_i y_i \sin{\theta_i}}{\sum_i y_i \cos{\theta_i}} \right), \end{align} \tag{ 12 }$$

where θ_i is the angle coded for by action cell i. It is also possible to compute the magnitude of the population vector, which denotes how strongly the action cell activities are promoting a particular heading,

$$\begin{align} m_{target} = \sqrt{\left( \sum_i y_i \sin{\theta_i} \right)^2 + \left( \sum_i y_i \cos{\theta_i} \right)^2}. \end{align} \tag{ 13 }$$

For practical reasons, the action cells are computed not only from place cell inputs but also by a separate module, termed a semi-random walk module. The reason for using such a module is that the network, particularly in the early stages of exploration when the weights are random, is often unable to make useable directional decisions. Therefore, implementing a semi-random walk module allows MiRo to explore the environment sensibly instead of erratically when we use the randomised network weights. Below, we provide the details of the semi random walk implementation.

2.2.2.1. Semi-random walk module

In cases where the signal provided by the action cells, as computed by equation (13) is not strong enough (i.e. less than 1), then MiRo takes a random walk than following the direction selected by the action cells. To compute the heading, a small but random value, θ_noise , is added to MiRo's current heading,

$$\begin{align} \theta_{random\_walk} = \theta_{current} + \theta_{noise}, \end{align} \tag{ 14 }$$

where θ_noise is a random variable taken from the uniform distribution $\theta_{noise} \sim \mathrm{unif}(-50^{\circ}, 50^{\circ})$. It ensures that MiRo generally keeps moving in its current direction but can change slightly to the left or right by no more than 50^∘. Generally, it is possible to achieve a random walk without implementing an explicit exception for low activity values (e.g. [59]). However, this requires careful tuning of the noise levels for its achievements and would lead to more erratic movement at the beginning of the training.

To convert this direction to action cell values, we compute each action cell as a function of its angular distance from $\theta_{random\_walk}$. We do this similarly to how we compute the place cell activities, i.e. as the Cartesian distance of MiRo from the place cell centres,

$$\begin{align} y_i^{random\_walk} = y_i^{max} \exp \left[ -\frac{\left( \theta_{random\_walk} - \theta_i \right)^2}{2\theta_d^2} \right], \end{align} \tag{ 15 }$$

where $y_i^{max}$ determines the maximum value for y_i , in this case 1, and θ_d determines the distribution width, and θ_i is the angle corresponding to action cell i.

To state this more formally, let the magnitude of the place cell network proposal be (see equation (13)),

$$\begin{align} m_{PC\_proposal} = \sqrt{\left( \sum_i \tilde{y}_i \sin{\theta_i} \right)^2 + \left( \sum_i \tilde{y}_i \cos{\theta_i} \right)^2}, \end{align} \tag{ 16 }$$

then the final action cell values are only changed to $y_i = y_i^{random\_walk}$ if $m_{PC\_proposal} \lt 1$. Else they stay as they are from equation (10).

2.2.2.2. Computing action cells during reverse replays

The computation for y_i in equation (10) is suitable for the exploration stage. Still, it requires a minor modification for the action cells to replay properly during reverse replay events. Thus far, y_i is computed either by taking the network's output as determined by the place cell inputs or, if this output is weak, by using a semi-random walk. For the y_i term to compute properly in the reverse replay case then, we perform the following,

$$\begin{align} & y_i^{replay}\nonumber\\ & \;\;\; = \frac{1}{1 + \mathrm{exp} \left[ -c_1 \sum_{j = 1}^{100} \left( w_{ij}^{\,PC\textrm{-}AC} + \omega~sgn( {e_{ij}^{r}}) \right) x_j - c_2 \right]}, \end{align} \tag{ 17 }$$

which is the same computation as equation (11), with the only difference being that we have added to the place cell to action cell weights the value ω = 0.1 multiplied by the sign of the eligibility trace for that synapse. The term $e_{ij}^{r}$ represents the value of e_ij , i.e. a trace of the potential synaptic changes at the moment of reward retrieval. This term effectively stores the history of synaptic activity and adds a transient weight increase to recently active synapses. We describe the computation of this eligibility trace in appendix section 'Mathematical derivation of the place-action cell synaptic learning rule'.

The scaling factor of 0.1 is heuristically selected so that the sign of the eligibility trace (i.e. $-1/+1$) will not over-dominate the weight term. A positive eligibility trace implies that weights should increase, and a negative eligibility trace indicates that weights should decrease. By adding 0.1 of the eligibility, we modulate the output of the striatal cells in the replay in producing more (or less) activity according to the desired direction. In other words, similar to the eligibility trace, the reply conveys information about the desirable synaptic change, which leads to improved performance.

Modifying the action cells during replays is necessary so that a reverse replay of the place cells can appropriately reinstate the activity in the action cells [18]. Without this change, the reverse replays would offer no additional benefits. This modification acts like a synaptic tag that activates at reward retrieval only and provides temporary synaptic modifications, according to the sign of the eligibility trace, during the reverse replay stage. Despite this assumption, this temporary change in synaptic strengths is similar to that of acetylcholine levels modifying synaptic conductances during replay events in the hippocampus [24]. In other words, synaptic weights (and their modifications) are suppressed during exploration but are manifest during the replay stage.

We also tested the rule using a weaker assumption, adding only a value of ω = 0.1 for any synapse in which $e_{ij} \gt 0$, whilst adding nothing for synapses where $e_{ij} \lt 0$. However, this performs worse than even the non-replay case. Since replays activate multiple cells simultaneously, neighbouring place and action cell activities will influence synapses that may have had $e_{ij} \lt 0$. This influence caused them to increase their weights instead of decreasing, which would be the proper direction given a negative value for e_ij .

The learning rule we derived is a policy-gradient RL method [57]. Its form is that of a three-factor learning rule with an eligibility trace [15]. The complete derivation for the learning rule is in the appendix.

When MiRo is exploring, a learning rule of the following form is active:

$$\begin{align} \frac{dw_{ij}^{\,PC\textrm{-}AC}}{dt} = \frac{\eta}{\sigma^2} R e_{ij}, \end{align} \tag{ 18 }$$

where R is a reward value, whilst the term e_ij represents the eligibility trace and is a time-decaying function of the potential weight changes, determined by,

$$\begin{align} \frac{de_{ij}}{dt} = -\frac{e_{ij}}{\unicode{x03C4}_e} + \left(y_i - \tilde{y}_i\right) \left( 1-\tilde{y}_i \right) \tilde{y}_i x_j. \end{align} \tag{ 19 }$$

During reverse replays, however, the action cells' target activity is given by $y_i^{replay}$, making this a supervised learning scenario. We, therefore, derived a learning rule structurally similar to the RL rule from the supervised learning framework (minimisation of an error function, see appendix), that is:

$$\begin{align} \frac{dw_{ij}^{\,PC\textrm{-}AC}}{dt} = \eta^{\prime} e_{ij}, \end{align} \tag{ 20 }$$

where the eligibility trace is determined by,

$$\begin{align} \frac{de_{ij}}{dt} = -\frac{e_{ij}}{\unicode{x03C4}_e} + \left(y_i^{replay} - \tilde{y}_i\right) \left( 1-\tilde{y}_i \right) \tilde{y}_i x_j. \end{align} \tag{ 21 }$$

We have set $\eta^{^{\prime}} = \eta / \sigma^2$ and let R = 1 at the reward location in our simulations, which renders the RL rule and the supervised learning rule equivalent.

We compute the population weight vector for a single place cell,

$$\begin{align} \left(w_j^{\,x}, w_j^{\,y}\right) = \left( \sum_{i = 1}^{72} w_{ij}^{\,PC\textrm{-}AC} \cos{\theta_i}\, , \, \sum_{i = 1}^{72} w_{ij}^{\,PC\textrm{-}AC} \sin{\theta_i} \right), \end{align} \tag{ 22 }$$

where $(w_j^{\,x}, w_j^{\,y})$ represents the x and y components for the weight population vector of the j^th place cell, $w_{ij}^{\,PC\textrm{-}AC}$ is the value of the weight from place cell j onto action cell i, and θ_i is the heading direction that action cell i codes for. The magnitude of the population weight vector is given by,

$$\begin{align} M_{w_j} = \sqrt{\left(w_j^{\,x}\right)^2 + \left(w_j^{\,y}\right)^2}. \end{align} \tag{ 23 }$$

The population weight vector depicts the preferred direction of MiRo when placed at the centre of the location of the place cell.

The entire implementation process is described here, with an overview of the algorithmic implementation presented in box 1.

2.2.5.1. Initialisation

At the start of a new experiment, the weights that connect the place cells to the action cells are initialised according to a uniform distribution and then normalised,

$$\begin{align} w_{ij}^{\,PC\textrm{-}AC} \leftarrow \frac{w_{ij}^{\,PC\textrm{-}AC}}{\sum_i w_{ij}^{\,PC\textrm{-}AC}}. \end{align} \tag{ 24 }$$

All the variables for the place cells are set to their steady-state conditions for when no place-specific inputs are present, and the action cells are all set to zero. MiRo is then placed in a random location in the arena.

2.2.5.2. Taking actions

2.2.5.3. Determining reward values

2.2.5.4. Initiating reverse replays

2.2.5.5. Updating network variables

All parameter values used in the Hippocampal network are in table 1, and those for the Striatal network in table 2. Values for η and $\unicode{x03C4}_e$ are specified appropriately in the results since they vary across different experiments. Unless otherwise stated in the text, parameters have been selected close to biological values where appropriate but tuned to achieve a good model performance.

Table 2. Summarising the model parameter values for the striatal network used in the experiments. Except for the learning rate, η, and the eligibility trace time constant, $\unicode{x03C4}_e$, all other parameters are kept constant for all experiments.

Parameter	Value
c1	0.1
c2	20
σ	0.1
θd	10
$\unicode{x03C4}_e$	See text
η	See text

The non-replay model requires the recent history to be stored in the eligibility trace. Having too short an eligibility trace time constant might negatively impact the model's performance. The time constant reflects how far back the information about the Reward will be 'transmitted'. Reverse replays, however, have the potential to compensate for this issue since the recent history is also stored and then replayed in the place cell network. Figure 3 shows the effects on performance of significantly reducing the eligibility trace time constant (to τ = 0.04 s). Both cases, with and without reverse replays, are compared. If the learning rate is too low (η = 0.01), then for neither case is there any learning. But as the learning rate increases, having reverse replays has significantly improved performance. Similar results are true for a larger, but not too large, eligibility trace time constant of $\unicode{x03C4}_e = 0.2$ s (see supplementary material). Replays offer the most significant advantage when the eligibility trace time constant, $\unicode{x03C4}_e$, is relatively small. As this time constant gets larger, replays offer little to no performance advantage over non-replays (see supplementary material) when the maximum learning time is 20 trials (but see section 3.2.4 for a higher number of trials).

Zoom In Zoom Out Reset image size

There is an interesting comparison between the magnitudes of weight changes for the replay and non-replay cases. Figure 4 shows the population vectors of the weights after reward retrieval. Population vectors for the weights are computed according to equations (22) and (23). There are two observations to be made here. First, the weight magnitudes are greater with reverse replays, which is not surprising since activity replay results in more synaptic changes. And second is that the direction of the population weight vectors themselves is slightly different, particularly in the location at the start of the trajectory. In particular, the weight vectors point more towards the goal location in the replay case, whereas the non-replay case has weight vectors pointing along the direction of the path taken by the robot. Whilst we depict only one case here, this is representative of several cases for various parameter values.

Zoom In Zoom Out Reset image size

We investigated the robustness of the performance across various values of $\unicode{x03C4}_e$ and η. Figure 5 displays the average performance over trials 11–20, comparing again with replays versus without replays. There are perhaps two noticeable observations to make here. Firstly, when the eligibility trace time constant is small, employing reverse replays shows considerable improvements in performance over the non-replay case across the various values of learning rates. Learning still exists in the non-replay case; however, it is noticeably diminished compared with the replay case. Secondly, although this marked performance improvement vanishes for larger eligibility trace time constants, reverse replays do not hinder performance at the very least. To better interpret the results, it is essential to recall that learning rate and eligibility trace are multiplicative terms defining the weight updates; see equation (19). The eligibility trace should be long enough to enable propagation to the initial locations of the trajectory. High values in synaptic changes can also lead to instabilities, so small learning rates compensate for this. In general, the most successful combinations are large eligibility trace time constants and low learning rates or vice versa. However, for the reverse replay case, the additional information from the replay helps to perform well even in the case of a small eligibility trace time constant, as long as the learning rate is sufficiently large.

Zoom In Zoom Out Reset image size

Figure 6 compares the results for the best cases with and without reverse replays. We optimised $\unicode{x03C4}_e$ and η independently for each case and ran 30 trials to achieve these results. The reason for this was a suspected instability in the non-replay case, i.e. a drop in performance as learning continues above the 20 trials. We interpret this as an advantage of the replay case in terms of stability rather than performance.

Zoom In Zoom Out Reset image size

In figure 2, where in trial 20 and above, the time to complete the task increases without replay but not in the case with replay. A Wilcoxon signed-rank test on the trials supported that data for the two conditions do not have the same median for 8 of the last 12 trials (p < 0.05, the complete table of results is in the appendix). The test failed to reject the null hypothesis that data have the same median in trials 0–18 (also evident by eye inspection).

We also note the significant difference in the optimal parameters for the best cases with/without replay. With reverse replays the parameters are $\unicode{x03C4}_e = 0.04$ s, η = 1, whereas without reverse replays they are $\unicode{x03C4}_e = 1$ s, η = 0.01. We speculate that the necessary choice in the eligibility time constant for the non-replay case (i.e. it needs to be large enough to store the trajectory history) underlines the cause of this instability. On the contrary, the reversal reply introduces additional synaptic changes, allowing for smaller eligibility time constants and helping the rule stability.

We derive a policy gradient rule [57] following [59], but here we use continuous valued neurons instead of spiking neurons. The expectation for the rewards earned in an episode of duration T is given by,

$$\begin{align} \left\langle R\right\rangle_T = \int_X \int_Y \ R \left({\mathbf{x}}, {\mathbf{y}}\right)P_w\left({\mathbf{x}},{\mathbf{y}}\right) d{\mathbf{y}}d{\mathbf{x}}, \end{align} \tag{ 25 }$$

where X is the space of the inputs of and Y the space of the output of the network, and $P_w\left({\mathbf{x}},{\mathbf{y}}\right)$ the probability that the network has input x and output y, parametrised by the weights.

We can decompose the probability, $P_w\left({\mathbf{x}},{\mathbf{y}}\right)$ (see decomposition of the probability in [59]) as,

$$\begin{align} P_w\left({\mathbf{x}},{\mathbf{y}}\right) = \prod_{j} g_j\left({\mathbf{x}},{\mathbf{y}}\right) \prod_i h_i\left({\mathbf{x}},{\mathbf{y}}\right), \end{align} \tag{ 26 }$$

where h_i is the probability the ith action cell generates output ${\mathbf{y}}_j$ contained in y, when the network receives input x. Similarly g_j is the probability for the activity produced by the jth place cell given its input. We then wish to calculate the partial derivative over a weight w_kl of the expected reward,

$$\begin{align} \frac{\partial \left\langle R_T\right\rangle}{\partial w_{kl}} = \int_X \int_Y \ R \left({\mathbf{x}}, {\mathbf{y}}\right) \frac{\partial P_w\left({\mathbf{x}},{\mathbf{y}}\right)}{\partial w_{kl}}~d{\mathbf{y}}d{\mathbf{x}} . \end{align} \tag{ 27 }$$

To do so, we take into account that $P_w\left({\mathbf{x}},{\mathbf{y}}\right) = \left[ \frac{P_w\left({\mathbf{x}},{\mathbf{y}}\right)}{h_k\left({\mathbf{x}},{\mathbf{y}} \right)} \right]h_k\left({\mathbf{x}},{\mathbf{y}}\right)$, where the term in square brackets does not depend on w_kl since we remove its contribution from $P_w\left({\mathbf{x}},{\mathbf{y}}\right)$ by dividing with $h_k\left({\mathbf{x}},{\mathbf{y}}\right)$. We can then write,

$$\begin{align} \frac{\partial P_w\left({\mathbf{x}},{\mathbf{y}}\right)}{\partial w_{kl}} = P_w\left({\mathbf{x}},{\mathbf{y}}\right)\frac{\partial~log~h_k\left({\mathbf{x}},{\mathbf{y}}\right)}{\partial w_{kl}}. \end{align} \tag{ 28 }$$

This leads to,

$$\begin{align} \frac{\partial \left\langle R_T\right\rangle}{\partial w_{kl}} = \int_X \int_Y \ R \left({\mathbf{x}}, {\mathbf{y}}\right) P_w\left({\mathbf{x}},{\mathbf{y}}\right)\frac{\partial~log~h_k\left({\mathbf{x}},{\mathbf{y}}\right)}{\partial w_{kl}}~d{\mathbf{y}}d{\mathbf{x}} . \end{align} \tag{ 29 }$$

To proceed, we need to consider the distribution of the activities of the action cells h_k . This we choose to be a Gaussian function with mean $\bar{y}_k$ and variance σ² (see also section 'Striatal Action Cells'),

$$\begin{align} h_k\left(X,Y\right) = \frac{1}{\sigma \sqrt{2\pi}}\exp{\left(-\frac{\left(y_k - \tilde{y}_k\right)^2}{2\sigma^2}\right).} \end{align} \tag{ 30 }$$

The mean of the distribution is calculated by $\tilde{y}_k = f_s\left( c_1 \sum_j w_{kj} x_j +c_2 \right)$, see also equation (11), where f_s is a sigmoidal function. We note that a different choice of function would have resulted in a variant of this rule. Therefore,

$$\begin{align} \frac{\partial~log~h_k\left({\mathbf{x}},{\mathbf{y}}\right)}{\partial w_{kl}} = c_1 \frac{y_k-\tilde{y}_k}{\sigma^2} \left( 1 - \tilde{y}_k \right) \tilde{y}_k x_l. \end{align} \tag{ 31 }$$

Replacing (31) in (29) we end up with,

$$\begin{align} \frac{\partial \left\langle R_T\right\rangle}{\partial w_{kl}}& = \int_X \int_Y c_1 \ R \left({\mathbf{x}}, {\mathbf{y}}\right) P_w\left({\mathbf{x}},{\mathbf{y}}\right) \frac{y_k-\tilde{y}_k}{\sigma^2}\nonumber\\ & \quad \times \left( 1 - \tilde{y}_k \right) \tilde{y}_k x_l ~d{\mathbf{y}}d{\mathbf{x}}. \end{align} \tag{ 32 }$$

Then the batch update rule is given by:

$$\begin{align} \frac{dw_{kl}}{dt}& = \eta \int_X \int_Y \ R \left({\mathbf{x}}, {\mathbf{y}}\right) P_w\left({\mathbf{x}},{\mathbf{y}}\right) \frac{y_k-\tilde{y}_k}{\sigma^2}\nonumber\\ & \quad \times \left( 1 - \tilde{y}_k \right) \tilde{y}_k x_l ~d{\mathbf{y}}d{\mathbf{x}}, \end{align} \tag{ 33 }$$

with the factor c₁ absorbed in the learning rate.

The batch rule indicates that we need to average the term $R \left({\mathbf{x}}, {\mathbf{y}}\right) \frac{y_k-\tilde{y}_k}{\sigma^2} \left( 1 - \tilde{y}_k \right) \tilde{y}_k x_l$ across many trials. When an on-line setting is considered, the average is naturally rising from sampling throughout the episodes. Hence the on-line version of this rule is given by:

$$\begin{align} \frac{dw_{kl}}{dt} = \eta~R \left({\mathbf{x}}, {\mathbf{y}}\right)~\frac{y_k-\tilde{y}_k}{\sigma^2} \left( 1 - \tilde{y}_k \right) \tilde{y}_k x_l. \end{align} \tag{ 34 }$$

We note however that this rule is appropriate for scenarios where reward is immediate. To deal with cases of distant rewards, such as ours where reward comes at the end of a sequence of actions, we need to resort to eligibility traces. Our rule is similar to REINFORCE with multiparameter distribution [65]; we differ by having a continuous time formulation and a different parametrisation of the neuronal probability density function. Further, in our case we do not learn the variance of the probability density function.

We introduce an eligibility trace by updating the weights connecting the place cells to the action cells, $W^{\,PC\textrm{-}AC}$ by:

$$\begin{align} \frac{dw_{ij}^{\,PC\textrm{-}AC}}{dt} = \frac{\eta}{\sigma^2} R \left({\mathbf{x}}, {\mathbf{y}}\right) e_{ij}. \end{align} \tag{ 35 }$$

The term e_ij represents the eligibility trace, see also [57], and is a time decaying function of the potential weight changes, determined by:

$$\begin{align} \frac{de_{ij}}{dt} = -\frac{e_{ij}}{\unicode{x03C4}_e} + \left(y_i - \tilde{y}_i\right) \left( 1-\tilde{y}_i \right) \tilde{y}_i x_j. \end{align} \tag{ 36 }$$

During replays, we assume that synapses between place and action cells change to minimise the function:

$$\begin{align} E = \frac{1}{2}\sum_i\left(y_i^{replay} - \tilde{y}_i\right)^2, \end{align} \tag{ 37 }$$

in other words, we assume that during the replay equation (17) provides a fixed target value for the mean of the Gaussian distribution of the action cells at time t. In what follows we consider the target constant for the sake of the derivation and consistency with the form of the RL rule, but in fact this target changes as time and consequently the weights from place to action cells change, making the rule unstable, but stabilising under a short, fixed length of reply time. Taking the gradient over the error function with respect to the weight w_kl , when considering the 'target' activity for the action cells fixed, leads to the backpropagation update rule for a single layer network:

$$\begin{align} \frac{dw_{kl}}{dt} = \eta^{\prime}~\left(y_k^{replay}-\tilde{y}_k\right) \tilde{y}_k \left( 1 - \tilde{y}_k \right) x_l, \end{align} \tag{ 38 }$$

where $\eta^{\prime}$ is the learning rule, in our simulations $\eta^{\prime} = \eta/\sigma^2$ similar to the RL rule. Also for consistency with the RL rule formulation, we introduce an eligibility trace by updating the weights connecting the place cells to the action cells, $W^{\,PC\textrm{-}AC}$ by:

$$\begin{align} \frac{dw_{ij}^{\,PC\textrm{-}AC}}{dt} = \eta^{\prime} e_{ij}, \end{align} \tag{ 39 }$$

where the eligibility trace is determined by:

$$\begin{align} \frac{de_{ij}}{dt} = -\frac{e_{ij}}{\unicode{x03C4}_e} + \left(y_i^{replay} - \tilde{y}_i\right) \left( 1-\tilde{y}_i \right) \tilde{y}_i x_j, \end{align} \tag{ 40 }$$

where again the time constant τ_e is the same as in the RL rule.

In the case of replays then, when the robot has reached its target, it first learns using the standard learning rule as in equations (35) and (36). After 1 s, a replay event is initiated, and learning is then done using the supervised learning rule here, using equations (39) and (40). By setting the reward value to R = 1, we can ensure that both the RL learning rule and the supervised learning rule become identical. Alternatively, a different value of the reward would require setting $\eta^{\prime} = \eta R /\sigma^2$.