Smooth and temporally continuous random motion can be difficult to model. To be smooth (and therefore physically plausible), a trajectory must be continuous in both position and velocity. To be temporally continuous, the statistics of the motion must be independent of the integration timestep being used. To be random, position and velocity at one time must not be reliable predictors of position and velocity at another time, provided these times are seperated by a sufficiently long interval. Implementations of random motion models typically fail to satisfy one, or sometimes two, of these principles (Raudies and Hasselmo, 2012; Benna and Fusi, 2021).

Ornstein-Uhlenbeck processes, which sit at the heart of the RatInABox random motion model, are continuous-in-time random walks with a tendency to return to a central drift value. The decorrelation timescale can be also be controlled. We use these to update the velocity vector (linear and rotational velocities are updated independently) on each update step. Position is then updated by taking a step along the velocity vector with some additional considerations to avoid walls. This method ensures both position and velocity are continuous, yet evolve ‘randomly’ (Figure 1a, d), and the statistics of the motion is independent of the size of the discretisation timestep being used.

Reanalysing rat locomotion data from Sargolini et al., 2006 (as has been done before, by Raudies and Hasselmo, 2012) we found that the histograms of linear speeds are well fit by a Rayleigh distributions whereas rotational velocities are approximately fit by normal distributions (Figure 2a). Unlike Raudies and Hasselmo, 2012, we also extract the decorrelation timescale of these variables and observe that rotational velocity in real locomotion data decorrelates nearly an order of magnitude faster than linear velocity (0.08 s vs. 0.7 s). We set the default parameters of our Ornstein-Uhlenbeck processes (including applying a transform on the linear velocity so its long-run distribution also follows a Rayleigh distribution, see Methods) to those measured from the Sargolini et al., 2006 dataset (Figure 2b).

Animals rarely charge head-first into a wall, turn around, then continue in the opposite direction. Instead, they slow down smoothly and turn to avoid a collision. Additionally, during random foraging, rodents are observed to show a bias towards following walls, a behaviour known as thigmotaxis (Satoh et al., 2011; Figure 2c). To replicate these observations, walls in the Environment lightly repel the Agent when it is close. Coupled with the finite turning speed this creates (somewhat counter-intuitively) a thigmotactic effect where the agent over-explores walls and corners, matching what is observed in the data (Figure 2e). A user-defined parameter called ‘thigmotaxis’ can be used to control the strength of this emergent effect (Figure 2d).

RatInABox supports importing trajectory data which can be used instead of the inbuilt random motion model. Imported trajectory data points which may be of low temporal-resolution are interpolated using cubic splines and smoothly upsampled to user-define temporal precision (Figure 3a). This upsampling is essential if one wishes to use low temporal resolution trajectory data to generate high temporal resolution neural data.

RatInABox supports online velocity control. At each integration step a target drift velocity can be specified, towards which the Agent accelerates. We anticipate this feature being used to generate complex stereotyped trajectories or to model processes underpinning complex spatial behaviour (as we demonstrate in Figure 3b, e).

Any single toolkit cannot contain all possible neural representations of interest. Besides, static cell types (e.g. PlaceCells, GridCells etc.) which have fixed receptive fields are limiting if the goal is to study how representations and/or behaviour are learned. RatInABox provides two solutions: Firstly, being open-source, users can write and contribute their own bespoke Neurons (instructions and examples are provided) with arbitrarily complicated rate functions. Secondly, two types of function-approximator Neurons are provided which map inputs (the firing rate of other Neurons) to outputs (their own firing rate) through a parameterised function which can be hand-tuned or trained to represent an endless variety of receptive field functions including those which are mixed selective, non-linear, dynamic, and non-stationary.

• FeedForwardLayer: Calculates a weighted linear combination of the input Neurons with optional bias and non-linear activation function.

• NeuralNetworkNeurons: Inputs are passed through a user-provided artificial neural network.

Naturally, function-approximator Neurons can be used to model how neural populations in the brain communicate, how neural representations are learned or, in certain cases, neural dynamics. In an online demo, we show how GridCells and HeadDirectionCells can be easily combined using a FeedForwardLayer to create head-direction selective grid cells (aka. conjunctive grid cells Sargolini et al., 2006). In Figure 3d and associated demo GridCells provide input to a NeuralNetworkNeurons class which is then trained, on data generated during exploration, to have a highly complex and non-linear receptive field. Function-approximator Neurons can themselves be used as inputs to other function-approximator Neurons allowing multi-layer and/or recurrent networks to be constructed and studied.

Efficiently encoding what an Agent can ‘see’ in its local vicinity, aka. its field of view, is crucial for many modelling studies. A common approach is to use a convolutional neural network (CNN) to process a rendered image of the nearby environment and extract activations from the final layer. However, this method is computationally expensive and necessitates training the CNN on a large dataset of visual images.

RatInABox offers a more efficient alternative through the use of VectorCells. Three variants – FieldOfViewBVCs, FieldOfViewOVCs, and FieldOfViewAVCs – comprise populations of egocentric Boundary-, Object-, and AgentVectorCells with angular and distance preferences specifically set to tile the Agent’s field of view. Being egocentric means that the cells remained fixed in the reference frame of the Agent as it navigates the Environment. Users define the range and resolution of this field of view. Plotting functions for visualising the field of view cells, as shown in Figure 3c, are provided.

In RatInABox, PlaceCells and VectorCells are sensitive to walls in the Environment. Three distance geometries are supported: ‘euclidean’ geometry calculates the Euclidean distance to a place field centre and so cell activity will ‘bleed’ through boundaries as if they weren’t there. ‘line_of_sight’ geometry allows a place cell to fire only if there is direct line-of-sight to the place field centre from the current location. Finally ‘geodesic’ geometry (default) calculates distance according to the shortest boundary-avoiding path to the cell centre (notice smooth wrapping of the third place field around the wall in Figure 1b). The latter two geometries respect the observation that place fields don’t typical pass through walls, an observation which is thought to support efficient generalisation in spatial reinforcement learning (Gustafson and Daw, 2011). Boundary conditions can be periodic or solid. In the former case, place fields near the boundaries of the environment will wrap around.

RatInABox simplifies the calculation and visualization of rate maps through built-in protocols and plotting functions. Rate maps can be derived explicitly from their known analytic firing functions or implicitly from simulation data. The explicit method computes rate maps by querying neuron firing rates at all positions simultaneously, utilizing ’array programming’ to rapidly compute the rate map. In the implicit approach, rate maps are created by plotting a smoothed histogram of positions visited by the Agent, weighted by observed firing rates (a continuous equivalent of a smoothed spike raster plot). Additionally, the tool offers the option to visualize spikes through raster plots.

Our random motion model is based on the Ornstein Uhlenbeck (OU) process, ${\displaystyle {\mathsf{X}}_{\theta ,\lambda ,\mu }(t)}$, a stochastic process satisfying the Langevin differential equation

where ${\displaystyle \eta (t)\sim \mathcal{N}(0,1)}$ is Gaussian white noise and ${\displaystyle \theta }$, ${\displaystyle \lambda }$ and μ are constants. The first term in the update equation drives decay of ${\displaystyle {\mathsf{X}}_{\theta ,\lambda ,\mu }(t)}$ towards the mean μ. The second term is a stochastic forcing term, driving randomness. These stochastic processes are well studied; their unconditioned covariance across time is

Thus ${\displaystyle {\mathsf{X}}_{\theta ,\lambda ,\mu }(t)}$ decorrelates smoothly over a timescale of ${\displaystyle \tau =1/\theta }$. Over long periods ${\displaystyle {\mathsf{X}}_{\theta ,\lambda ,\mu }(t)}$ is stochastic and therefore unpredictable. Its long-run stationary probability distribution is a Gaussian with mean μ and standard deviation ${\displaystyle \sigma =\sqrt{{\lambda }^2/2\theta }}$. We can re-parameterise the Ornstein Uhlenbeck process in terms of these more intuitive parameters (the decoherence timescale ${\displaystyle \tau }$ and the long-run standard deviation ${\displaystyle \sigma }$) using the transformations

to give

Ornstein Uhlenbeck processes have the appealing property that they are temporally continuous (their statistics are independent of ${\displaystyle dt}$) and allow for easy control of the long-run standard deviation and the decoherence timescale of the stochastic variable. For these reasons, we use use them to model rotational and linear velocities within RatInABox.

To avoid travelling through walls, if a collision is detected the velocity is elastically reflected off the wall (normal component is flipped). The speed is then scaled to one half the average motion speed, ${\displaystyle {\mathsf{v}}_{2D}(t)=0.5{\sigma }_{\mathsf{v}}}$.

Spring-deceleration model. In order to slow down before colliding with a wall the Agent feels an acceleration, perpendicular to the wall, whenever it is within a small distance, ${\displaystyle d_{\text{wall}}}$, of the wall.

${\displaystyle d_{\perp ,j}(t)}$ is the perpendicular distance from the Agent to the ${\displaystyle j^{\text{th}}}$ wall, ${\displaystyle {\mathbf{n}}_j}$ is the perpendicular norm of the ${\displaystyle j^{\text{th}}}$ wall (the norm pointing towards the Agent) and ${\displaystyle k_1}$ & ${\displaystyle s}$ are constants (explained later). ${\displaystyle d_{\text{wall}}}$ is the distance from the wall at which the Agent starts to feel the deceleration, defaulting to ${\displaystyle d_{\text{wall}}=0.1}$ m.

Note that this acceleration is identical to that of an oscillating spring-mass where the base of the spring is attached a distance ${\displaystyle d_{\text{wall}}}$ from the wall on a perpendicular passing through the Agent. The spring constant is tuned such that a mass starting with initial velocity towards the wall of ${\displaystyle -s{\sigma }_{\mathsf{v}}{\mathbf{n}}_j}$ would stop just before the wall. In summary, for ${\displaystyle k_1=1}$, if the Agent approaches the wall head-on at speed of ${\displaystyle s{\sigma }_{\mathsf{v}}}$ (${\displaystyle s}$ times its mean speed) this deceleration will just be enough to avoid a collision.

${\displaystyle s}$ is the unitless wall repel strength parameter (default ${\displaystyle s=1}$). When it is high, walls repel the agent strongly (only fast initial speeds will result in the agent reaching the wall) and when it is low, walls repel weakly (even very slow initial speeds will not be slowed done by the spring dynamics). When ${\displaystyle s=0}$ wall repulsion is turned off entirely.

Conveyor-belt model. A second (similar, but not exactly equivalent) way to slow down motion near a wall is to consider a hypothetical conveyor belt near the wall. This conveyor belt has a non-uniform velocity pointing away from the wall of

When the Agent is close to the wall the hypothetical conveyor-belt moves it backwards on each time step, effectively slowing it down. Note that this velocity is identical to that of a spring-mass attached to the wall with initial velocity ${\displaystyle s{\sigma }_{\mathsf{v}}{\mathbf{n}}_j}$ away from the wall and spring constant tuned to stop the mass just before it reaches a distance ${\displaystyle d_{\text{wall}}}$. In summary, for ${\displaystyle k_2=1}$, if the Agent approaches the wall head-on at speed of ${\displaystyle s{\sigma }_{\mathsf{v}}}$ the conveyor belt will just be fast enough to bring it to a halt at the location of the wall.

Wall attraction (thigmotaxis). Although similar, there is an exploitable difference between the ‘spring-deceleration’ and ‘conveyor-belt’ models: the ‘conveyor-belt’ changes the Agents position, ${\displaystyle \mathbf{x}(t)}$, on each step but not its internal velocity variable ${\displaystyle \mathbf{v}(t)}$. As as result (and as the conveyor-belt intuition suggests) it will slow down the Agent’s approach towards the wall without causing it to turn around. This creates a ‘lingering’ or ‘thigmotactic’ effect whereby whenever the Agent heads towards a wall it may carry on doing so, without collision, for some time until the stochastic processes governing its motion (section ‘Temporally continuous random motion’) cause it to turn. Conversely the ‘spring-deceleration’ model has no ‘thigmotactic’ effect since it actively changes the internal velocity variable causing the Agent to turn around or ‘bounce’ off the walls.

The relative strengths of these two effects, ${\displaystyle k_1}$ and ${\displaystyle k_2}$, are controlled by a single thigmotaxis parameter, ${\displaystyle {\lambda }_{\text{thig}}\in [0,1]}$ which governs the trade-off between these two models.

When ${\displaystyle {\lambda }_{\text{thig}}=1}$ only the conveyor belt model is active giving a strong thigmotactic effects. When ${\displaystyle {\lambda }_{\text{thig}}=0}$ only the spring-deceleration model is active giving no thigmotactic effect. By default ${\displaystyle {\lambda }_{\text{thig}}=0.5}$. The constants 3 and 6 are tuning parameters chosen by hand in order that direct collisions with the walls are rare but not impossible.

Although this procedure, intended to smoothly slow the Agent near a wall, may seem complex, it has a two advantages: Firstly, deceleration near walls is smooth, becoming stronger as the Agent gets nearer and so induces no physically implausible discontinuities in the velocity. Secondly, it provides a tunable way by which to control the amount of thigmotaxis (evidenced in Figure 2c, d). Recall that these equations only apply to motion very near the wall (${\displaystyle {\displaystyle <d_{\text{wall}}}}$) and they can be turned off entirely (${\displaystyle s=0}$) (see Table 1, below).

A set of locations (the centre of the place fields), ${\displaystyle \{{\mathbf{x}}_i^{\text{PC}}\}}$, is randomly sampled from the Environment. By default these locations sit on a grid uniformly spanning the Environment to which a small amount of random jitter, half the scale of the sampled grid, is added. Thus, place cell locations appear ‘random’ but initialising in this way ensures all parts of the Environment are approximately evenly covered with the same density of place fields.

The environmental-distance from the Agent to the place field centres is calculated (${\displaystyle {\displaystyle d_i(t)=d({\mathbf{x}}_i^{\text{PC}},\mathbf{x}(t))}}$). The firing rate is then determined by one of the following functions (defaulting to ${\displaystyle F^{\text{{gaussian}}}}$):

Where used, ${\displaystyle w_i}$ is the user-provided radius (aka. width) of the place cells (defaulting to 0.2 m).

Each grid cell is assigned a random wave direction ${\displaystyle {\theta }_i\sim {\mathcal{U}}_{[0,2\pi ]}}$, gridscale ${\displaystyle {\lambda }_i\sim {\mathcal{U}}_{[0.5\text{m},1.0\text{m}]}}$ and phase offset ${\displaystyle {\varphi }_i\sim {\mathcal{U}}_{[0,2\pi ]}}$. The firing rate of each grid cell is given by the thresholded sum of three cosines

${\displaystyle {\mathbf{e}}_{\theta }}$ is the unit vector pointing in the direction ${\displaystyle \theta }$. We also provide a shifted (as opposed to rectified) sum of three cosines grid cell resulting in softer grid fields

${\displaystyle {\mathbf{e}}_{\theta }}$ is the unit vector pointing in the direction ${\displaystyle \theta }$.

VectorCells subclasses include BoundaryVectorCells, ObjectVectorCells and AgentVectorCells as well as FieldOfView versions of these three classes. The common trait amongst all types of VectorCell is that each cell is responsive to a feature of the environment (boundary segments, objects, other agents) at a preferred distance and angle. The firing rate of each vector cell is given by the product of two functions; a Gaussian radial function and a von Mises angular function. When the agent is a euclidean distance ${\displaystyle d(t)}$ from the feature, at an angle ${\displaystyle \varphi (t)}$ the contribution of that feature to the total firing rate is given by

where ${\displaystyle f_{\text{VM}}}$ is the radial von Mises distribution (a generalisation of a Gaussian for periodic variables)

Total firing rate is calculated by summing/integrating these contributions over all features in the Environment as described in the following sections. Distance and angular tuning parameters and defined/sampled as follows:

• ${\displaystyle d_i}$ is the distance tuning of the vector cell. By default ${\displaystyle d_i\sim {\mathcal{U}}_{[0.05\text{m},0.3\text{m}]}}$

• ${\displaystyle {\sigma }_{d,i}}$ is the distance tuning width. By default this increases linearly as a function of ${\displaystyle d_i}$: ${\displaystyle {\sigma }_{d,i}=d_i/\beta +\xi }$ for constants ${\displaystyle \beta }$ and ${\displaystyle \xi }$ but can be set otherwise. See Table 1, below for the values which are chosen to match those used by Cothi and Barry (de Cothi and Barry, 2020).

• ${\displaystyle {\varphi }_i}$ is the angular tuning of the vector cell. By default ${\displaystyle {\varphi }_i\sim {\mathcal{U}}_{[0^{\circ },{360}^{\circ }]}}$.

• ${\displaystyle {\sigma }_{\varphi ,i}}$ (which defines the von Mises concentration measure ${\displaystyle {\kappa }_i:=1/\sqrt{{\sigma }_{\varphi ,i}}}$) is the angular tuning width of the vector cell. By default ${\displaystyle {\sigma }_{\varphi ,i}\sim {\mathcal{U}}_{[{10}^{\circ },{30}^{\circ }]}}$.

The asymptotic equivalence between a Gaussian and a von Mises distribution (true for small angular tunings whereby von Mises distributions of concentration parameter ${\displaystyle \kappa }$ approach Gaussian distributions of variance ${\displaystyle {\sigma }^2=1/\kappa }$) means this model is effectively identical to the original boundary vector cell model proposed by Hartley et al., 2000 but with the difference that our vector cells (BVCs included) will not show discontinuities if they have wide angular tunings of order ${\displaystyle {360}^{\circ }}$.

All vector cells can be either

• allocentric (default): ${\displaystyle \varphi (t)}$ is the angle subtended between the x-direction vector ${\displaystyle {\mathbf{e}}_x=[1,0]}$, and the line between the Agent and the feature.

• egocentric: ${\displaystyle \varphi (t)}$ is the angle subtended between the heading direction of the agent ${\displaystyle \hat{\mathbf{h}}(t)}$, and the line between the Agent and the feature.

The environmental features which BoundaryVectorCells (BVCs) respond to are the boundary segments (walls) of the Environment. The total firing rate of of each cell is given by integrating (computationally we use a default value of ${\displaystyle d\theta =2^{\circ }}$ to numerically approximate this integral) the contributions from the nearest line-of-sight boundary segments (walls occluded by other walls are not considered) around the full ${\displaystyle 2\pi }$ field-of-view;

(computationally we use a default value of ${\displaystyle d\theta =2^{\circ }}$ to numerically approximate this integral). ${\displaystyle K_i=1/\underset{\mathbf{x}}{max}{F}_i(\mathbf{x})}$ is a normalisation constant calculated empirically at initialisation such that each BVC has a maximum firing rate (before scaling) of 1.0 Hz.

ObjectVectorCells (OVCs) respond to objects in the Environment. Objects are zero-dimensional and can be added anywhere within the Environment, each object, ${\displaystyle j}$, comes with a “type” attribute, ${\displaystyle t_j}$. Each object vector cell has a tuning type, ${\displaystyle t_i}$, and is only responsive to objects of this type. The total firing rate of of each cell is given by the sum of the contributions from all objects of the correct type in the Environment;

Since Equation 33 has a maximum value of 1 by definition the maximum firing rate of an object vector cell is also 1 Hz (unless multiple objects are closeby) and no normalisation is required.

AgentVectorCells respond to other Agents in the Environment. All cells in a given class are selective to the same Agent, index ${\displaystyle j}$. The firing rate of each cell is then given by;

FieldOfViewBVCs/OVCs/AVCs are a special case of the above vector cells where the tuning parameters (${\displaystyle d_i}$, ${\displaystyle {\sigma }_{d,i}}$, ${\displaystyle {\varphi }_i}$, ${\displaystyle {\sigma }_{\varphi ,i}}$) for a set of VectorCells are carefully set so that cells tile a predefined ‘field of view’. By default these cells are egocentric and so the field of view (as the name implies) is defined relative to the heading direction of the Agent; if the Agent turns the field of view turns with it.

Users define the angular and radial extent of the field of view as well as the resolution of the cells which tile it. There is some flexiblity for users to construct complex fields of view but baic API simplifies this process, exposing a few key parameters:

• ${\displaystyle r_{\text{fov}}=[r_{\text{fov}}^{\text{min}},r_{\text{fov}}^{\text{\% max}}]}$ (default [0.02 m, 0.2 m]): the radial extent of the field of view.

• ${\displaystyle {\theta }_{\text{fov}}}$ (default [0°, 75°]): the angular extend of the field of view (measured from the forward heading direction, symmetric left and right).

• ${\displaystyle {\delta }_{\text{fov}}^0}$ (default 0.02 m): FieldOfView VectorCells all have approximately circular receptive fields (i.e. the radial Gaussian and angular von Mises in Equation 33 have matched variances which depend on their tuning distance; ${\displaystyle {\sigma }_{d,i}=d_i\cdot {\sigma }_{\varphi ,i}:={\delta }_{\text{fov}}(d_i)}$). ${\displaystyle {\delta }_{\text{fov}}^0}$ sets the resolution of the inner-most row of cells in the field of view, ${\displaystyle {\delta }_{\text{fov}}^0={\delta }_{\text{fov}}(d_i=r_{\text{fov}}^{\text{min}})}$.

• Manifold type: For “diverging” manifolds (default) cells further away from the Agent have larger receptive fields ${\displaystyle {\delta }_{\text{fov}}(d_i)={\xi }_0+d_i/\beta }$ for user-defined ${\displaystyle \beta }$ (default ${\displaystyle \beta =5}$) and ${\displaystyle {\xi }_0:={\delta }_{\text{fov}}^0-r_{\text{fov}}^{\text{min}}/\beta }$. For “uniform” manifold all cells have the same sized receptive fields, ${\displaystyle {\delta }_{\text{fov}}(d_i)={\delta }_{\text{fov}}^0}$.

More complex field of views can be constructed and a tutorial is provided to show how.

In 2D Environments each head direction cell has an angular tuning mean ${\displaystyle {\theta }_i}$ and width ${\displaystyle {\sigma }_i:=1/\sqrt{{\kappa }_i}}$. The response function is then a von Mises in the head direction of the Agent:

By default all cells have the same angular tuning width of 3° and tuning means even spaced from 0° to ${\displaystyle {360}^{\circ }}$.

In 1D Environments there is always and only exactly ${\displaystyle n=2}$ HeadDirectionCells; one for leftward motion and one for rightward motion.

VelocityCells are a subclass of HeadDirectionCells which encode the full velocity vector rather than the (normalised) head direction. In this sense they are similar to HeadDirectionCells but their firing rate will increase with the speed of the Agent.

In 2D their firing rate is given by:

where ${\displaystyle {\theta }_{\mathsf{v}}(t)}$ is the angle of the velocity vector ${\displaystyle \mathbf{v}(t)}$ anticlockwise from the x-direction and ${\displaystyle {\sigma }_{\mathsf{v}}}$ is the likely speed scale of the Agent moving under random motion (this is chosen so the firing rate of the velocity cell before scaling is approximately ${\displaystyle \mathsf{O}(1)}$ Hz).

In 1D environments:

where the addition of ${\displaystyle {\mu }_{\mathsf{v}}}$ accounts for any bias in the motion.

A single cell encodes the scaled speed of the Agent

where, same as with the VelocityCells, ${\displaystyle {\sigma }_{\mathsf{v}}}$ (or ${\displaystyle {\sigma }_{\mathsf{v}}+{\mu }_{\mathsf{v}}}$ in 1D) is the typical speed scale of the Agent moving under random motion giving these cells ad pre-scaled maximum firing rate of ${\displaystyle \mathsf{O}(1)}$ Hz.

PhasePrecessingPlaceCells (a subclass of PlaceCells) display a phenomena known as phase precession with respect to an underlying theta oscillation; within each theta cycle the firing rate of a place cell peaks at a phase dependent on how far through the place field the Agent has travelled. Specifically, as the Agent enters the receptive field the firing rate peaks at a late phase in the cycle and as the Agent leaves the receptive field the firing rate peaks at an early phase in the cycle, hence the name phase precession. Phase precession is implemented by modulating the spatial firing rate of PlaceCells with a phase precession factor, ${\displaystyle F_i^{\theta }(t)}$,

which rises and falls each theta cycle, according to:

This is a von Mises factor where ${\displaystyle {\varphi }_{\theta }(t)=2\pi {\nu }_{\theta }t\mathrm{mod}2\pi }$ is the current phase of the ${\displaystyle {\nu }_{\theta }}$ Hz theta-rhythm and ${\displaystyle {\varphi }_i^{\ast }(\mathbf{x}(t),\widehat{\dot{\mathbf{x}}}(t))}$ is the current ‘preferred’ theta phase of a cell which is a function of it’s position ${\displaystyle \mathbf{x}(t)}$ and direction of motion ${\displaystyle \widehat{\dot{\mathbf{x}}}(t)}$. This preferred phase is calculated by first establishing how far through a cells spatial receptive field the Agent has travelled along its current direction of motion;

and then mapping this to a uniform fraction ${\displaystyle {\beta }_{\theta }}$ of the range ${\displaystyle [0,2\pi ]}$;

${\displaystyle {\sigma }_i}$ is the width of the cell at its boundary, typically defined as ${\displaystyle {\sigma }_i=w_i}$, except for gaussian place cells where the boundary is arbitrarily drawn at two standard deviations ${\displaystyle {\sigma }_i=2w_i}$.

The intuition for this formula can be found by considering an Agent travelling straight through the midline of a circular 2D place field. As the Agent enters into the receptive field (at which point ${\displaystyle (\mathbf{x}(t)-{\mathbf{x}}_i)\cdot \widehat{\dot{\mathbf{x}}}(t)=-{\sigma }_i}$) the firing rate will peak at a theta phase of ${\displaystyle \pi +\beta \pi }$. This then precesses backwards as it passes through the field until the moment it leaves (${\displaystyle (\mathbf{x}(t)-{\mathbf{x}}_i)\cdot \widehat{\dot{\mathbf{x}}}(t)={\sigma }_i}$) when the firing rate peaks at a phase of ${\displaystyle \pi -\beta \pi }$. This generalises to arbitrary curved paths through 2D receptive fields. This model has been used and validated before by Jeewajee et al., 2014 . ${\displaystyle {\kappa }_{\theta }}$ determines the spread of the von Mises, i.e. how far from the preferred phase the cell is likely to fire.

RandomSpatialNeurons provide spatially ‘tuned’ inputs for use in instances where PlaceCells, GridCells, BoundaryVectorCells etc. These neurons have smooth but, over long distances, random receptive fields (approximately) generated by sampling from a Gaussian process with a radial basis function kernel of lengthscale ${\displaystyle l}$ (default ${\displaystyle l=0.1}$ m). The kernel is given by:

where ${\displaystyle d(\mathbf{x},{\mathbf{x}}^{\mathrm{\prime }})}$ is the environmental-distance between two points in the environment. This distance measure (same as used for PlaceCells, and VectorCells etc.) accounts for walls in the environment and so the receptive fields of these neurons are smooth everywhere except across walls (see Section ‘Distance measures’).

Firing rates are calculated as follows: At initialisation an array of target locations, at least as dense as the lengthscale, is sampled across the environment ${\displaystyle \{{\mathbf{x}}_j\}}$. For each neuron, ${\displaystyle i}$, ${\displaystyle j}$ target values, ${\displaystyle [{\tilde{F}}_i{]}_{:}}$, is sampled from the multivariate Normal distribution

where ${\displaystyle \mathbf{K}}$ is the covariance matrix with elements ${\displaystyle K_{lm}=k({\mathbf{x}}_l,{\mathbf{x}}_m)}$. This creates a sparse set of locations, ${\displaystyle \{{\mathbf{x}}_j\}}$, and targets, ${\displaystyle {\tilde{F}}_{ij}}$, across the Environment: locations close to each other are likely to have similar targets (and hence similar firing rates) whereas locations far apart will be uncorrelated.

At inference time the firing rate at an arbitrary position in the Environment, ${\displaystyle \mathbf{x}(t)}$ (which will not neccesarily be one of the pre-sampled targets) is estimated by taking the mean of the targets weighted by the kernel function between the position and the target location:

This weighted average is a cheap and fast approximation to the true Bayesian Gaussian process which would require the inversion of the covariance matrix ${\displaystyle \mathbf{K}}$ at each time-step and which we find to be numerically unstable around exposed walls.

FeedForwardLayer and NeuralNetworkNeurons are different from other RatInABox classes; their firing rates are not textitexplicitly determined by properties (position, velocity, head direction etc.) of their Agent but by the firing rates of a set of input layers (other ratinabox.Neurons). They allow users to create arbitrary and trainable ‘function approximator’ Neurons with receptive fields depending non-trivially on the states of one or many Agent(s).

Each FeedForwardLayer has a list of inputs ${\displaystyle \{{\mathsf{L}}_j{\}}_{j=1}^N}$ which must be other ratinabox.Neurons subclasses (e.g. PlaceCells, BoundaryVectorCells, FeedForwardLayer). For input layer ${\displaystyle j}$ with ${\displaystyle n_j}$ neurons of firing rates ${\displaystyle F_k^{{\mathsf{L}}_j}(t)}$ for ${\displaystyle k\in [1,n_j]}$, a weight matrix is initialised by drawing weights randomly ${\displaystyle {\mathsf{w}}_{ik}^{{\mathsf{L}}_j}\sim \mathcal{N}(0,g/\sqrt{n_j})}$ (for default weight intialisation scale ${\displaystyle g=1}$). The firing rate of the ${\displaystyle i^{\text{th}}}$ FeedForwardLayer neuron is given by weighted summation of the inputs from all layers plus a bias term:

where ${\displaystyle \varphi (x)}$ is a potentially non-linear activation function defaulting to a linear identity function of unit gain. ${\displaystyle b_i}$ is a constant bias (default zero). A full list of available activations and their defining parameters can be found in the utils.py file; these include ReLU, sigmoid, tanh, Retanh, softmax and linear (the default) functions or users can pass their own bespoke activation function.

Alongside ${\displaystyle \varphi (r_i(t))}$ this layer also calculates and saves ${\displaystyle {\varphi }^{\mathrm{\prime }}(r_i(t))}$ where ${\displaystyle {\varphi }^{\mathrm{\prime }}}$ is the derivative of the activation function, a necessary quantity for many learning rules and training algorithms.

NeuralNetworkNeurons are a generalisation of FeedForwardLayer. Like FeedForwardLayer they are initialised with a list of inputs ${\displaystyle i^{\text{th}}}$. This class also recieves, at the point of initialisation, a neural network, NN. This can be any pytorch.nn.module. To calculate teh firing rate this class takes the firing rates of all input layers, concatenates them, and passes them through the neural network. The firing rate of the NeuralNetworkNeurons neuron is given by the activity of the neuron in the output layer of neural network:

If no neural network is provided by the user a default network with two hidden ReLU layers of size 20 is used.

In order to be compatible with the rest of the RatInABox API the firing rate returned by this class is a numpy array, however, on each update the output of the pytorch neural network is additionally saved as a torch tensor. By accessing this tensor, users can take gradients back through the embedded neural network and train is as we demonstrate in Figure 3e.

In Figure 3e and an associated demo script a NeuralNetworkNeurons layer is initialised with ${\displaystyle N=1}$ neuron/output. The inputs to the network come from a layer of 200 GridCells, ranging in grid scale from 0.2 m to 0.5 m. These are passed through a neural network with three hidden ReLU layers of size 100 and a linear readout. As the Agent randomly explores its Environment the network is trained with gradient descent to reduce the L2 error between the firing rate of the network and that of a ‘target’ rate map (a vector image of the letters ‘RIAB’). We use gradient descent with momentum and a learning rate of ${\displaystyle \eta =0.002\cdot {\text{dt}}^2}$ (which makes the total rate of learning time-step independent). Momentum is set to ${\displaystyle \mu =(1-\frac{dt}{{\tau }_{\text{et}}})}$ where ${\displaystyle {\tau }_{\text{et}}}$ is the eligibility trace timescale of 10 s which smoothes the gradient descent, improving convergence. We find learning converges after approximately 2 hr and a good approximation of the target function is achieved.

See the GitHub repository for instructions on how to install RatInABox. The following is a Python script demonstrating a very basic use-case.

Import RatInABox and necessary classes. Initialise a 2D Environment. Initialise an Agent in the Environment. Initialise some PlaceCells. Simulate for 20 s. Print table of times, position and firing rates. Plot the motion trajectory, the firing rate timeseries’ and place cell rate maps.

# Import RatInABox 
import ratinabox from ratinabox.Environment import Environment 
from ratinabox.Agent import Agent 
from ratinabox.Neurons import PlaceCells 
import pandas as pd

# Run a very simple simulation 
Env = Environment() 
Ag = Agent(Env) 
PCs = PlaceCells(Ag) 
for i in range(int(20/Ag.dt)): 
Ag.update() 
PCs.update()

# Export data into a dataframe 
pd.DataFrame(Ag.history)

# Plot data 
Ag.plot_trajectory() 
PCs.plot_rate_timeseries() 
PCs.plot_rate_map()

Jupyter script: https://github.com/RatInABox-Lab/RatInABox/blob/main/demos/paper_figures.ipynb.

In this demonstration we study, using RatInABox, which type of spatially modulated cell type is best for decoding position.

First we generate training and testing datasets. A set of Neurons (${\displaystyle n=N_{\text{cells}}=20}$) is initialised in a 1 m square Environment containing a small barrier (Appendix 1—figure 1, top). A six minute trajectory is simulated using the RatInABox random motion model to produce a dataset of inputs ${\displaystyle \{{\mathbf{x}}_i{\}}_{i=1}^{N_T}}$ and targets ${\displaystyle \{{\mathbf{y}}_i{\}}_{i=1}^{N_T}}$:

where ${\displaystyle \mathbf{x}(t_i)}$ is the position of the Agent at time ${\displaystyle t_i}$ and ${\displaystyle \overrightarrow{F}}$ is the firing rate of the neuronal population. These data are split into training (${\displaystyle {\displaystyle 0<t_i<5}}$ mins, Appendix 1—figure 1a purple) and testing (${\displaystyle {\displaystyle 5<t_i<6}}$ mins, Appendix 1—figure 1a black) fractions. The goal of the decoder is to learn a mapping ${\displaystyle G:\mathcal{X}\to \mathcal{Y}}$, from firing rates to positions.

To do this we use Gaussian Process Regression (GPR). GPR is a form of non-parameteric regression where a prior is placed over the infinite-dimensional function space ${\displaystyle P(G(\mathbf{x}))}$ in the form of its covariance kernel ${\displaystyle C(\mathbf{x},{\mathbf{x}}^{\mathrm{\prime }})}$ and mean ${\displaystyle \mu (\mathbf{x})}$ (typically zero). This defines a prior on the targets in the training set ${\displaystyle \mathbf{Y}=({\mathbf{y}}_1,{\mathbf{y}}_2,{\mathbf{y}}_3,...{)}^{\mathsf{T}}}$,

where ${\displaystyle {\mathsf{C}}_{ij}=C({\mathbf{x}}_i,{\mathbf{x}}_j)+{\sigma }_{\eta }{\delta }_{ij}}$ is a covariance matrix established over the data points. The second term accounts for additive noise in the data function. This can be used to make an inference on the posterior of the target for an unseen testing data point, ${\displaystyle P({\mathbf{y}}_{\text{test}}|\{{\mathbf{x}}_i{\}}^{\text{train}},\{{\mathbf{y}}_i{\}}^{\text{train}},{\mathbf{x}}_{\text{test}})}$ – itself a Gaussian – the mean of which is taken as the “prediction”. A more comprehensive reference/tutorial on Gaussian Process Regression is given by MacKay, 2003.

We use a radial basis function (aka “squared exponential”) kernel with width ${\displaystyle l=l_0\sqrt{N_{\text{cells}}}}$ which scales with the expected size of the population vector (${\displaystyle \sim \sqrt{N_{\text{cells}}}}$, we set ${\displaystyle l_0=\sqrt{20}}$)

and a small amount of target noise ${\displaystyle {\sigma }_{\eta }=1e-10}$. Note that the closest ‘parameterised’ analog to GPR with an RBF kernel is linear regression against Gaussian basis features of length scale ${\displaystyle l}$. Since the Gaussian is a non-linear function this means our regression prior is also non-linear function of firing rate (and therefore potential non-biologically plausible). We choose to optimise with the sklearn.gaussian_process.GaussianProcessRegressor package. Note we do not attempt to optimise the hyperparameters ${\displaystyle l_0}$ or ${\displaystyle {\sigma }_{\eta }}$ which one would probably do in a more rigorous experiment. RatInABox parameters are all default with the exception that the place cells are of type gaussian_threshold and width ${\displaystyle w_i=0.4}$ m and the timestep is set to ${\displaystyle dt=50}$ ms.

Appendix 1—figure 1b (lower) show the results over comparable sets of PlaceCells, GridCells and BoundaryVectorCells. Coloured dots show the prediction – mean of the posterior – of the GPR model “trained” on all points in the training dataset for that particular cell type. This is plotted on top of the true trajectory, shown in black. PlaceCells perform best achieving under 1 cm average decoding error, followed by BoundaryVectorCells then GridCells where the decoded position is visibly noisy.

Place cells outperform grid cells which outperform BVCs irrespective of how many cells are using in the basis feature set. More cells gives lower decoding error. Decoding errors in Appendix 1—figure 1c are smaller than would be expected if one decoded from equivalently sized populations of real hippocampal neurons. There are likely many reasons for this. Real neurons are noisy, communicate sparsely through spikes rather than rates and, most likely, jointly encode position and many other behaviourally relevant (or irrelevant) variables simultaneously. All of these factors could be straightforwardly incorporated into this analysis using existing RatInABox functionality.

Appendix 1—figure 1

Download asset Open asset

Jupyter script: https://github.com/RatInABox-Lab/RatInABox/blob/main/demos/reinforcement_learning_example.ipynb.

In this example we demonstrate how RatInABox can be used in a reinforcement learning (RL) study. The goal is as follows: train an artificial Agent to explore a 2D Environment where a reward is hidden behind a wall. The Agent should become proficient at navigating around the wall and towards the reward from all locations within the Environment.

The core of our approach will rest on model-free RL where an Agent first learns a value function for a policy (a process known as “policy evaluation”) and then uses this value function to define a new, improved policy (“policy improvement”). Iterating between these two procedures (“policy iteration”) can result in convergence towards an optimal or near-optimal policy.

A core pillar of RatInABox is its continuous approach to modelling time and space. This continuity will require revising typical approaches to how the value function if defined, approximated and then learned, as well as how motion control (aka action selection, in discrete space) is performed. This is not a weakness, in fact we would argue it is one of the strengths. Once we are done, we are left with a formulation of model-free RL which bears much higher resemblance to biological navigation. Furthermore, since most of the complexities of feature encoding and motion control in continuous time and space are handled by RatInABox innately this “upgrade” comes almost for free.

The value of a motion policy, ${\displaystyle \pi }$, is defined as the decaying sum (or integral in continuous time) of expected future rewards

where the expectation is taken over any stochasticity present in the current policy (i.e. how the Agent moves) and Environment/reward (although in this case both will be deterministic). This definition of value is temporally continuous. The key differences compared to the more common form – where value is written as a discrete sum of rewards over future timesteps – is that it is now a continuous integral over a reward density function and temporal discounting is done by exponentially decaying future reward over a time period ${\displaystyle \tau }$. The prefactor of ${\displaystyle 1/\tau }$ is an optional constant of normalisation.

In order to learn the value function we define a new ValueNeuron class. The ValueNeuron, which is a subclass of FeedForwardLayer, recieves feedforward input from a set of features corresponding to PlaceCells scattered across the Environment with firing rates ${\displaystyle \{{\varphi }_i{\}}_{i=1}^{N_{\varphi }=1000}}$ where ${\displaystyle {\varphi }_i(\mathbf{x})=F_i(\mathbf{x})}$ is the firing rate of the ${\displaystyle i^{\text{th}}}$ place cell at location ${\displaystyle \mathbf{x}}$. This linear approximation to the value function can be written as