Noise-driven bifurcations in a neural field system modelling networks of grid cells

The activity generated by an ensemble of neurons is affected by various noise sources. It is a well-recognised challenge to understand the effects of noise on the stability of such networks. We demonstrate that the patterns of activity generated by networks of grid cells emerge from the instability of homogeneous activity for small levels of noise. This is carried out by analysing the robustness of network activity patterns with respect to noise in an upscaled noisy grid cell model in the form of a system of partial differential equations. Inhomogeneous network patterns are numerically understood as branches bifurcating from unstable homogeneous states for small noise levels. We show that there is a phase transition occurring as the level of noise decreases. Our numerical study also indicates the presence of hysteresis phenomena close to the precise critical noise value.


Introduction
By now it is well established that grid cells, and the characteristic hexagonal firing patterns they create in physical space, play an important role in the navigational system of mammalian brains (McNaughton et al. 2017). Since grid cells were discovered in Hafting et al. (2005), there has been extensive activity in order to understand their precise behaviour, see (Rowland et al. 2016;McNaughton et al. 2017) and the references therein. The main challenges ten years after the discovery of grid cells, such as how grid cells are organised and how they are connected to other cell types in the brain, were highlighted in Rowland et al. (2016). In particular, as the brain is inherently noisy (Rolls and Deco 2010), the lack of understanding of the effect of noise on grid cells was emphasised as a challenge. The recent results in Gardner et al. (2022) has provided insight into the organisation of grid cells by showing that the activity of the network (called a module in Gardner et al. 2022) is arranged on a torus. The question regarding the effects of noise on grid cells, however, remains open.
In accordance with previous experimental studies and general belief in the field, the results in Gardner et al. (2022) provided further evidence in favour of describing the grid cell network by continuous attractor network dynamics through a system of neural field models (Ermentrout and Terman 2010). The first attractor network models for grid cells were presented in McNaughton et al. (2006), Burak and Fiete (2009), Couey et al. (2013), which were based on the classical papers Cowan (1972, 1973) and Amari (1977), see also Pinto et al. (1996). In Burak and Fiete (2009), Couey et al. (2013), the grid cells are assumed to have orientation preferences in four different directions. The hexagonal grid cell patterns are then generated by a system of 4N 2 neural field ordinary differential equations with β = 1, . . . , 4. Here s β i ≥ 0 represents the activity level of neuron i with orientation preference β, and τ β i is its relaxation time. The right-hand side of (1.1) represents the firing rate of the neuron, see Bressloff (2012). The function is a given activation function, often of the form of a ReLU or sigmoid function. The firing rate of neuron i depends on an external input B β i (t) and the response of the network. It is assumed that the neurons are arranged on a square, which we will denote and call the neural sheet, according to the strength of their pairwise connection. The position of neuron i is denoted by x i . The strength of the connectivity between neuron i of type β and j of type β is W β i j = W (x i − x j − r β ), where W (x), x ∈ , is assumed to be even in each coordinate, in . The connectivity is shifted in the direction of the orientation preference of neuron j of type β with r β which is given by shifts of equal length in the four cardinal directions: north, south, east and west. It has been commonly considered, and, as mentioned, recently shown in Gardner et al. (2022), that the network of neurons creates a torus connectivity. This is realised in the model by assuming that W is extended periodically outside .
By modelling the movement of a rat traversing physical space through the input B β i (t), (1.1) can recreate the hexagonal patterns in physical space produced by the firing of a single grid cell as observed in experiments (Couey et al. 2013). In this model, the patterns in physical space are a consequence of the patterns generated on the neural sheet . However, (1.1) being a deterministic model, it does not offer much insight into the effects of noise.
Works on understanding noisy neural fields have in general been lacking (Bressloff 2012, Sect. 6) until recently (Touboul 2012;Kilpatrick and Ermentrout 2013;Kilpatrick 2014;MacLaurin and Bressloff 2020;Touboul et al. 2012;Bressloff 2019;Byrne et al. 2019). In Burak and Fiete (2012) fundamental limits on how information dissipates in networks of noisy neurons were derived. The author in Kilpatrick (2014) presents a study of two coupled noisy neural field models with a focus on the consequences of the coupling on the neural activity waves, while Kilpatrick and Ermentrout (2013) and MacLaurin and Bressloff (2020) investigate the effect of noise on stationary bumps in one-dimensional spatially extended networks.
Taking a different perspective than Burak and Fiete (2012); Kilpatrick and Ermentrout (2013); Kilpatrick (2014) and MacLaurin and Bressloff (2020), by adding noise to the common model of a grid cell network (1.1), the main goal of this work is to analyse the robustness of the hexagonal patterns in the activity level (Ermentrout and Cowan 1979) observed in (1.1) with respect to noise strength. We show that the stationary spatial patterns of the activity level emerge from the instability of homogeneous brain activity as the noise level diminishes. By upscaling the model (1.1) with noise to a system of Fokker-Planck-like partial differential equations, our analysis gives an estimate on the noise strength above which there is no coherent activity pattern. We also numerically explore the different branches of inhomogeneous stationary patterns bifurcating from the homogeneous state depending on the noise for several activation functions indicating the presence of hysteresis phenomena.
Instabilities of homogeneous steady states of noisy neural fields were also investigated by Byrne et al. (2019) utilising a partial differential equation (PDE) description. However, the PDEs were of a very different form than the ones presented in this manuscript. A Fokker-Planck-like system describing a network of noisy neurons can be found in Bressloff (2019), where neural variability in a coupled ring network was studied.
It is classical to analyse the behavioural change of neural fields without noise in the form of ordinary differential equations (ODEs) by standard bifurcation analysis (Murray 2002;Bressloff 2012;Veltz et al. 2015;Kilpatrick and Poll 2017;Schmidt and Avitabile 2020). Finding noise-driven bifurcations is more challenging, and one has to rely on other technical tools unless the coupling of the network has a particular structure where closed ODEs for the mean and the variance are available Touboul 2012). The system (1.1) with noise has, in addition, a nonlinear coupling, and the activity levels must remain nonnegative due to their physical interpretation, leading to technical additional constraints on the stochastic processes involved. In the following we deal with these challenges by analysing a system of Fokker-Plancklike partial differential equations with boundary conditions describing the space-time evolution of the law of the stochastic processes with respect to the noise level.

The PDE system: derivation, main goal, and numerical experiment
For the sake of the reader, we start by discussing the simplest classical case of no spatial connectivity, see Hopfield (1984); Bressloff (2012) and the references therein. Let us consider the classical neural field stochastic dynamical system for a network of M coupled neurons given by Here, the neurons are considered indistinguishable and all-to-all coupled with equal strengths given by W 0 ∈ R whose sign depends on the type of neurons considered: inhibitory or excitatory. We also consider that the relaxation time for all neurons is the same and equal to τ . B(t) is the external input for this neural network and σ > 0 is the strength of the noise W k . We have considered independent Brownian motion for each neuron in the network. Classical stochastic analysis implies that we can derive a Fokker-Planck equation for the evolution of the probability density of neurons with activity level s at time t in the large population limit M → ∞, i.e., the law of the limiting stochastic process follows the PDE where f = f (t, s) denotes the probability to observe the activity s at time t, and f denotes the mean value of the activity level s Notice that the noise can drive the activity level to be negative in (2.1), which is clearly not desirable from the modelling viewpoint. In order to avoid this, it is common practise to consider the Fokker-Planck equation (2.2) on s ∈ [0, ∞) with no-flux boundary conditions This ensures that particles cannot escape from non-negative values of the activity level variable s at the PDE level while keeping an evolution of a probability density, see Carrillo et al. (2011), Carrillo et al. (2013) for instance.
Remark 2.1 (Microscopic Model) Reflective boundary conditions for stochastic processes have been incorporated at the stochastic differential equation (SDE) level in order to avoid particles to escape a fixed domain (Sznitman 1984;Lions and Sznitman 1984;Faugeras and Inglis 2015). One can produce a microscopic stochastic process by adding an additional process counting when particles touch the boundary of the domain. The law of the rigorous mean-field limit, M → ∞, of the following system . . , M, follows the evolution of (2.2)-(2.3) under suitable smoothness assumptions on , see Lions and Sznitman (1984), Sznitman (1991).
The next step in the modelling is to reinterpret M as the number of neurons in each of the cortical columns of a neural sheet of N columns. Given space points x 1 , . . . , x N in the region of the neural cortex, the interaction among N M neurons stacked in N columns at locations x i with M neurons each, where s β ik represents the activity level with orientation β of the k th neuron at location x i is given for i = 1, . . . , N and k = 1, . . . , M by Here, we consider the same periodic setting, imposed through the periodicity of the interactions W β for β = 1, . . . , 4, as in Burak and Fiete (2009), Couey et al. (2013). Moreover, the neurons are inhibitory (Couey et al. 2013) and the activity in the network is modulated by a time dependent external input as in (1.1). We are dealing with a population network of neurons structured by their orientation preference β = 1, . . . , 4 corresponding to the four cardinal points (north, west, south, east). The network population includes the localised in space cross inhibition of neurons with different orientations modulated by the shape function Following the approach outlined above in the case of one population, we can formally write a Fokker-Planck type equation, in the limit N , M → ∞, for the evolution of the probability density f β (t, x, s) of finding neurons of type β at position x on the neural sheet with activity level s ≥ 0 at time t ≥ 0. We refer to Cai et al. (2006) for a similar approach in conductance-voltage models. The system of equations reads where β (x) is given by periodic boundary conditions in x, and the no-flux boundary conditions at s = 0 given by for each position x in the square sheet . To realise the torus connectivity, we assume that W is periodic with respect to and even in each coordinate on . The function is typically a smooth approximation of the ReLU activation function (x) + = max{0, x} or a sigmoid function. The initial probability density of the system (2.7) is denoted by f β 0 .

Remark 2.2
The system of Fokker-Planck equations (2.7)-(2.8) can be rigorously derived from the microscopic stochastic processes (2.6) under suitable assumptions. The rigorous proof of this mean-field limit for the spatially extended system (2.7)-(2.8) has recently been obtained in Carrillo et al. (2021) by a generalisation of the coupling method of Sznitman (1991). This rigorous passage to the limit is a very interesting area of mathematical research on its own with a multitude of different models and limiting systems derived under different assumptions on the ingredients of the network. For instance, we refer to the works Moynot and Samuelides (2002), Faugeras et al. (2009), Faugeras andInglis (2015), Touboul (2012), Touboul et al. (2012), and Cabana and Touboul (2018) in which the authors deal with spatially extended systems of neural networks modelled by their voltage with random connectivity interactions using large deviation principles (Arous and Guionnet 1995;Guionnet 1997).
To summarise, the main goal of this work is to focus on the biological information carried by the system of PDEs (2.7)-(2.8). More precisely, we study how noise affects the dynamics of (2.7)-(2.8) under the following assumptions: (A1) The grid cells are arranged on a torus, realised by setting the neural sheet = [−0.5, 0.5] 2 and extending W periodically outside . (A2) The inhibitory (Couey et al. 2013) connectivity function W ≤ 0 is at least in L 2 ( ), and is an even function in each coordinate on . Furthermore, we define The modulation function is in C 1 (unless otherwise stated).
It is well-known that grid cell firing is strongly connected to mammals' navigation, but unknown exactly how the grid cell network communicates with other networks in the brain. We will therefore simply assume in the numerical experiments that the external input B β (t) in (2.7)-(2.8) depends on the velocity at time t, v(t), of a moving animal in the following manner (Burak and Fiete 2009;Couey et al. 2013): where B > 0 is a constant external excitatory input, assumed to be the same for different β, α the velocity modulation, θ(t) the orientation of the animal at time t according to the reference frame, and θ β the orientation preference of the neurons of type β (θ 1 = π 2 , θ 2 = π, θ 3 = 3 2 π, θ 4 = 2π ). This particular form of the input, together with the right set of parameters in (1.1), has been shown to enable single cells of the ODE system (1.1) to create hexagonal firing patterns in physical space, see Burak and Fiete (2009), Couey et al. (2013).

Numerical reproduction of the hexagonal patterns
We numerically demonstrate that the PDE system (2.7)-(2.8) with (2.9) is able to reproduce the characteristic single-cell hexagonal firing pattern as discovered by Hafting et al. (2005) for rats and see how this pattern depends on the noise strength σ > 0. For this, we use a numerical scheme that has been extensively utilised for Fokker-Planck like equations (Carrillo et al. 2015). For more details on the numerical approach and its validation, see Appendix 1. Before connecting the grid cell system (2.7)-(2.8) with the movement of a rat, we initialise the activity on the neuronal sheet by running the simulation with α = 0 until f β has numerically stabilised into stationary patterns equal to the ones in the top and middle rows of Fig. 1, modulo translations.
Then we connect with the rats movement by setting α = 0.3 in (2.9) as in Couey et al. (2013). The velocity v(t) and orientation θ(t) used in the numerical experiment are calculated using timestamped position data from the physical experiments in Hafting et al. (2005) where rats moved around in a circular enclosure with a radius of 80 cm. The shape of the enclosure can be seen in the plots in the bottom row of Fig. 1. The path of the rat after moving around for t = 5 minutes, generated with the position data, is visualised in grey.
The red coloured areas in each plot in the bottom row in Fig. 1 make up the firing field of a grid cell in the network. The firing field consists of smaller red circular areas, which again are made up by even smaller red dots. Each dot marks a firing of the grid cell placed at position (0, 0) on the neuronal sheet as the velocity and orientation data is fed into the numerical method of (2.7)-(2.8) through (2.9). For simplicity, we have assumed in the numerical experiment that a neuron at a particular position x ∈ fires as soon as the firing rate in (2.7)-(2.8) satisfies (x) > 0. Now, from left to right in the bottom row of Fig. 1, we see the pattern of the firing fields-the patterns created by the red dots over the path of the rat in the enclosure in physical space-for the grid cell at (0, 0) for increasing noise strength.
The top and middle rows in the figure display snapshots of the probability density f β at s = 0.5 and s = 0, respectively. As t increases, these patterns are translated in accordance with the movement of the rat.
As can be observed in the bottom row, the red fields generated with the PDE system (2.7)-(2.8) form hexagonal patterns similar to the ones observed in physical experiments. However, the distance between the activity bumps and the area they cover decreases as the noise strength increases. The second main observation is that by increasing the noise the firing becomes less and less localised. The numerical experiments support the existence of a critical value of the noise, σ c > 0, at which a single cell could fire no matter where the rat is on its path.
The question we will address in the following is how stable the patterns observed in Fig. 1 are with respect to the noise strength σ .

Stability of the neural field system
In this section, we start by studying the spatially homogeneous solutions. We would like to understand the pattern formation in the neural field system (2.7)-(2.8) as a byproduct of the instability of these homogeneous solutions. We assume that solutions to (2.7) are sufficiently smooth and decay fast enough as s → ∞. We further assume from now on that B β = B is constant, i.e., α = 0 in (2.9), in order to study the emergence of stationary network patterns of the system. Note that setting B β to different constant values depending on β could yield non-stationary network patterns: a rat running with constant speed in one direction would give v(t) = const > 0 and θ(t) = const in (2.9), which results in different constant values of B β . This would, with the right set of parameter values, consequently translate the network patterns in time. To avoid this technicality, we let B be identical for the four different direction preferences. Homogeneous with no-flux boundary conditions (2.8) at s = 0, i.e., and In order to find stationary spatially homogeneous states f ∞ we first assume that their mean f ∞ is given. Denote by 0 = W 0 f ∞ + B the corresponding firing rate for simplicity. Thus, by integrating (3.1) and using the boundary condition (3.2), the stationary spatially homogeneous states f ∞ (s) are given by with the mass normalisation factor Z such that where the error function has been defined as However, note that (3.3) is an implicit equation as 0 depends on the mean f ∞ . To show the existence of stationary solutions, we need to solve the consistency equation for the mean f ∞ given by We prove next that the stationary state exists and is unique by leveraging on (3.5) under suitable conditions on the firing rate function .
Proof Define for m ≥ 0 and σ > 0 the function It is not difficult to check that the supremum of g(η) over η ∈ [0, ∞) is given by Our assumptions on and W 0 ≤ 0 imply that ∂G ∂m (m, σ ) < 0, and then we obtain the desired unique zero of G defining our stationary state f ∞ through (3.3)-(3.5).

Remark 3.2
Notice that the previous proposition can also be applied for functions admitting negative values as exp(−η 2 ) 1+erf(η) behaves like − √ πη in the limit η → −∞, such that one can show for any η ∈ R. For instance, the theorem is valid for the ε-approximation of (x) = (x) + defined by for ε small enough such that ε (x) ≥ 1 W 0 . The smooth approximation˜ ε (x) = 0.5(x + √ x 2 + ε) of (x) = (x) + can also be used as˜ ε (x) trivially satisfies the hypotheses of Proposition 3.1 since˜ ε (x) > 0, and it is strictly increasing.
We have numerically analysed the stability of the stationary solutions obtained in the previous result among spatially homogeneous solutions of (3.1). In Fig. 2a, we illustrate that the computed stationary state and the numerical solution to the evolution problem after time t = 150ms are indistinguishable for the firing rate ε with ε = 0.01. In Fig. 2b, we observe the convergence in time towards the stationary state by computing the difference in L 1 and the difference in average between the stationary solution and the evolution problem. We conclude that the stationary state and the corresponding numerical solution to the evolution problem (3.1) are identical to machine precision after t = 150ms, and that the convergence in time is exponential with the rates computed by averaging over 100 runs.
Next, we focus on the linear stability of the spatially homogeneous solution as a solution of the nonlinear system (2.7)-(2.8). For comparison, we first find the condition for linear spatial stability in the case of zero noise (σ = 0) following the classical approach as in Murray (2002). LetŴ (k) be the two-dimensional Fourier transform of W restricted to ,Ŵ As an example, if W (x) is the characteristic of a ball with radius r fully contained in , then Remark 3.3 Given the assumptions on W in (A1) and (A2), the function F(k) is realvalued when the shifts r β , β = 1, 2, 3, 4, satisfy the assumptions in (A4). With shifts of equal size z, one can check that β exp − ik · r β = 2 cos(2π k 1 z) + 2 cos(2π k 2 z) with k = 2π k 1 2π k 2 using Euler's formula.
The following lemma presents a linear stability condition for the system (2.7) without noise.

Remark 3.4
The proof of the mean field limit in Carrillo et al. (2021) relies on σ > 0. The mean field limit in the case σ = 0 is easier to obtain and leads to a pure Vlasov equation. This is a very classical result in the smooth setting and is derived via estimates in transport distances, see Hauray and Jabin (2015), Cañizo et al. (2011), Jabin (2014), Hauray and Jabin (2007), Golse (2016) and the references therein. (3.7). Then the mean of the zero noise, i.e., σ = 0, spatially homogeneous, stationary solution of (2.7)-(2.8) is linearly asymptotically stable if F(k) < 1.

Lemma 3.5 Assume (A1)-(A4). Let F(k) be as in
Proof By taking the mean of (2.7)-(2.8) with σ = 0, we find that the mean at position , x), evolves according to (dropping the t dependence for ease of notation) where h β = s β − s ∞ and 0 = (W 0 s ∞ + B). Applying the ansatz h β (x, t) ∝ exp(λt + ik · x) to the equation, we find that each mode k has the characteristic polynomial We see that three of the eigenvalues are stable as long as τ > 0. The fourth eigenvalue is which determines whether the linear system is stable or not. The eigenvalue is negative if F(k) < 1.
We now turn to the case with noise, σ > 0. Let h β = f β − f ∞ , where f ∞ is defined through (3.3)-(3.5). Notice that for the perturbations to be admissible, we need to ensure that ∞ 0 f β ds = 1, and consequently ∞ 0 h β ds = 0.
In principle one needs f β ≥ 0 for it to be a probability density. However, we will prove below that linear stability holds without any assumption on the sign of f β .
After linearising (2.7) around the spatially homogeneous state f ∞ , we get with β = 1, . . . , 4. We now restrict the set of perturbations in L 2 ( × [0, ∞)) to the ones of the form where u β k (s, t) is sufficiently smooth. We can then reduce (3.10) to one Fourier mode. Dropping the subscript k we set v β (x, s, t) = exp(ik · x)u β (s, t), where u β may be complex-valued, and U = β exp(−ik · r β )u β , such that (3.10) turns into for β = 1, . . . , 4. An equation for the time evolution of U can also derived, namely where F(k) is defined in (3.7). Denote where the dependence on σ enters through f ∞ defined by (3.3)-(3.5). We remark that we cannot obtain closed ODE equations for the moments of the distribution in the activity variable s, and thus a similar analysis as in Touboul (2012); Touboul et al. (2012) is not possible here. Building on (3.12), we can prove the following result. Recall that U = β exp(−ik · r β )u β . To obtain a time decaying estimate from the bound, we need to separate the linear part of U from the nonlinear (Part II). Finally, we establish the stability of U , and consequently u β , which then yields the asymptotic stability of f ∞ (Part III).

Theorem 3.6 Assume (A1)-(A4). Let satisfy the assumptions in Proposition 3.1 and let F(k) be as in (3.7) and satisfy the condition in Remark 3.3. Then the spatially homogeneous steady solution f ∞ to (2.7)-(2.8) is linearly asymptotically stable in
Part I: Note that We multiply (3.13) withŪ (the complex conjugate of U ), and integrate over [0, ∞) with respect to f ∞ : After applying (3.17) to the first term on the right-hand side and integrating the last term by parts, we find that U satisfies (where (z) denotes the real part of the complex Using the boundary condition on the second term above and by again utilizing the equivalence (3.17), we get Performing an integration by parts on the second to last integral and then yet again using (3.17), From (3.12) and the definition of U , it can be shown that f ∞ U follows with F(k) as in (3.7). We add the two equalities, By setting α = 1 and β = −F(k), we get In the above derivation we have used that By applying the Cauchy-Schwarz inequality to ∞ U s , we see that the right-hand side of (3.19) is non-positive. However, it is not straightforward to determine whether the right-hand side is strictly negative or if there is a Grönwall type decay estimate to be obtained from this expression. In particular, the Cauchy-Schwarz inequality with the chosen functions whenever u 2 = c(t)u 1 for any function c(t), which here means that U = c(t)(s − f ∞ ) when we make sure that the requirement ∞ 0 U f ∞ ds = 0 holds. Part II: To separate the linear part from the nonlinear part, we split where M ∞ is defined in (3.14). We can summarise this as Inserting this into the first square in (3.21) and then expanding the square, (3.21) turns into We now apply the Poincaré inequality with respect to the measure f ∞ (s)ds (Muckenhoupt 1972;Roustant et al. 2017), given by In the above calculation, we have used the orthogonality of V and c(t)(s − f ∞ ) with respect to f ∞ ds. Part III: Defining due to (3.15). Note that also due to (3.15), D(t) ≥ 0. It can be checked after some tedious computations using the explicit expression of f ∞ in (3.3)-(3.5) that R < 0 We now choose α > 0 such that C 2 > 0, and then apply Poincaré's inequality to the integral, where the exponential decay (3.23) is applied in the last step. This leads to One can avoidC 2 = C 3 by choosing α appropriately.
The asymptotic stability of f ∞ in L 2 × [0, ∞) for the set of perturbations given by (3.11) now follows by an application of Parseval's identity with respect to x to the identity (3.11), and the estimate above. Notice that f ∞ (s) ≤ 1 Z from (3.4). Remark 3.8 In principle, the linear stability analysis is valid only for smooth . However, the stability condition (3.15) of the linearised problem does only depend on such that the result holds for the linearised system with (x) = (x) + . Notice that condition (3.15) is continuous with respect to the regularisation parameter ε in Remark 3.2 for which the linearisation is valid.
We also remark that the value of R in the proof above remains strictly negative as long as ε is small enough despite the fact that ε may give negative values.

Bifurcation diagrams and phase transitions
With our linear stability analysis at hand, we will now investigate how the stationary patterns of the nonlinear PDE system (2.7)-(2.8) change as we vary the noise parameter σ . We do this by numerically computing bifurcation branches from the spatially homogeneous solution f ∞ for various choices of the modulation function . The numerical procedure is described in more detail in Sect. 4.2.

Bifurcations and phase transitions of the nonlinear PDE system
We now continue with our examination of the stability of the spatial patterns, generated by the full system (2.7)-(2.8), with respect to σ . In Fig. 5, we numerically compute bifurcation branches from the spatially homogeneous solution f ∞ for different modulation functions with W (|x|) = −0.005 · 128 2 (1 + tanh(10 − 50|x|)) for the nonlinear problem (2.7)-(2.8). This is done by using a continuation based method on σ over an accurate numerical solver for the evolution in time of Fokker-Planck like equations developed in Carrillo et al. (2015); further details are given in Appendix 1. The continuation method starts either at the largest or the smallest noise value σ of the interval under consideration and it solves for the evolution in time of (2.7)-(2.8) up to stabilisation to a steady value. This allows for recursive computation of the stationary states for smaller or larger values of the noise by taking as initial data the already computed steady state. With this procedure we ensure, up to numerical accuracy, that we compute the stable stationary states, either by sweeping the noise values from left-to-right (l2r) or from right-to-left (r2l).
Each subplot in Fig. 5 shows the maximum and minimum over space x of the average activity rate f (x) = β f β (x) of the computed steady states for each noise value σ . We show both the spatial maximum and minimum of f (x) to illustrate the fact that the computed stationary states are not spatially homogeneous, in other words, that they lead to spatial patterns. We also plot the spatially homogeneous branch numerically solving the implicit expression (3.5) as reference. The red dots indicate the stability threshold in σ for the condition F(k) < 1, as in Lemma 3.5, to hold.
In Fig. 5a, we observe the bifurcation branches for the sigmoid function (x) = 1/(1 + exp(−15x)). All of them show a sharp discontinuity at different noise values. We restrict the discussion to the lines representing the spatial maximum. We first focus on the full line (l2r) and the dashed line (r2l) that connect two bifurcation branches at different noise values corresponding to a hexagonal-like pattern similar to Fig. 6a. This clearly indicates that there is a discontinuous phase transition near the noise value indicated by the arrow. The fact that the l2r and r2l curves do not coincide further indicates that there is a hysteresis phenomenon. This conclusion is supported by the fact that the blue dot, the minimum noise value for linear stability (3.15), is to the left of both branches. This allows the possibility of branches of dynamically unstable steady states bending backwards in noise at the phase transition point. Unstable branches are not computable with our numerical approach. Finally, we find a second bifurcation branch given by the dotted line (l2r-s) in Fig. 5c corresponding to a stripe-like pattern similar to Fig. 6b. This branch was found by imposing a particular symmetry on the initial data, i.e., enforcing a horizontal band with activity level one.
In Fig. 5b-d, we show analogous computations for the case of the modulation function given by (x) = (x) + and its regularisations (x) = ε (x) with ε = 0.1 and ε = 0.01. Similarly to Fig. 5a, we observe a discontinuous phase transition for the full line (l2r) and the dashed line (r2l), and the linear stability blue dot is also to the left of the phase transition point as above. We remark that the blue and the red dots may lie outside the noise intervals in Fig. 5b-d, but they follow the same order. The case of ε = 0.1 in Fig. 5b resembles the behaviour observed for the sigmoid function in Fig. 5a. The hexagonal-like patterns are the preferred stable configurations both for generic initial data, full line (l2r), and starting with small perturbations of the homogeneous stationary state, dashed line (r2l). Again stripe-like patterns are obtained by choosing specific initial data. Similar branches and the numerical observation that the hexagonal-like pattern is the most stable configuration has already been reported for a neural field model without noise (Veltz et al. 2015).
This behaviour changes in Fig. 5c, d. The hexagonal-like patterns are still the preferred stable configurations for generic initial data, full line (l2r). However, starting with small perturbations of the homogeneous stationary state, dashed line (r2l), we connect to the stripe-like bifurcation branch, dotted line (l2r).
The bifurcation branches and their dynamics gets richer as the regularisation parameter gets smaller. We observe that for ε = 0.01 in Fig. 5c, there is an additional branch, dark red line (l2r-e), leading to eye-like patterns as in Fig. 6c. This branch jumps to the stripe-like pattern for larger noise values. It is difficult to extract information on the range of noise values σ ∈ [0.025, 0.0265] since the branch, dark red line (l2r), does not show a sharp transition point while the dashed line (r2l) does. However, this becomes much clearer in the limiting case of the positive part in Fig. 5d. We observe two sharper discontinuous transition points in the stripe-like branch, lead- Fig. 6 Stationary patterns of f β (s = 0) at σ = 0.022 for (x) = ε (x), ε = 0.01 with W (|x|) = −0.005 · 128 2 (1 + tanh(10 − 50|x|)): hexagonal-like (a), stripe-like (b), eye-like (c). Left to right: black, dotted, and red line in Fig. 5c ing to an intermediate pure-stripe branch, σ ∈ [0.025, 0.028], before jumping to the homogeneous state, see the dashed line (r2l) and the dotted line (l2r).

Concluding remarks
The first conclusion of our analysis is that the mean-field limit of the grid cell model (1.1) with constant external input introduced by Burak and Fiete (2009), Couey et al. (2013), and Burak and Fiete (2012), presents phase transitions driven by the noise strength as demonstrated in Figs. 5 and 6. This behaviour resembles the phenomena appearing in the classical Kuramoto model (Kuramoto 1981;Sakaguchi et al. 1988;Acebrón et al. 2005;Carrillo et al. 2019Carrillo et al. , 2020 for synchronisation and other neural field models Touboul et al. (2012) in the computational neuroscience literature.
It is shown that the homogeneous in space stationary state is linearly unstable for small noise strength, similarly to basic ring and neural field models (MacLaurin and Bressloff 2020; Kilpatrick and Ermentrout 2013;Byrne et al. 2019). We numerically analysed the bifurcation diagram of stationary patterns showing the appearance of different branches identified by their symmetries, see Figs. 5 and 6. Our numerical experiments with random initial data demonstrate that the stationary hexagonal-like pattern in space of the activity level of neurons in Fig. 6a, leading to the solid black bifurcation branches in Fig. 5, has the largest basin of attraction. Moreover, the numerical simulations indicate that there is a sharp transition in the mean activity level together with a hysteresis phenomenon suggesting a discontinuous phase transition. Whether more stationary network patterns exist is another interesting topic.
The crucial implication of this phase transition on the rats' navigation path is that the larger the noise the less localised are the spatial firing fields of each neuron. This can be observed in the bottom row of Fig. 1 which shows that the firing field (coloured in red) gets denser as the noise increases. Moreover, there is a sharp value of the noise after which there is no localisation at all, leading the rats to not being able to orientate themselves in physical space. In other words, the point of transition from a homogeneous pattern (all neurons have the same mean activity level) to a nonhomogeneous pattern (neurons at different locations in the network have different mean activity levels) gives an upper bound for the noise strength for which single grid cells no longer can fire in a hexagonal pattern in physical space when connected with the rats movements through (2.9).
With a hexagonal network activity configuration as in Fig. 6a, a single neuron can create hexagonal neural field patterns in physical space as in the third row of Fig. 1. Exactly how the firing fields in physical space of a single neuron are affected by initial network activity patterns as the ones in Fig. 6b-c remains to be investigated.
From the methodological viewpoint, we remark that as the bifurcation branches are computed using a numerical approximation of the PDE (2.7)-(2.8), they can differ slightly from the actual branches of the PDE itself. To study bifurcations and phase transitions of the nonlinear PDE analytically will require sophisticated mathematical tools, and they will be investigated elsewhere.
From the computational neuroscience viewpoint, we expect that noise driven phase transitions will also naturally appear in related attractor dynamic models as the ones in Burak and Fiete (2012), Agamon and Burak (2020). Instability of homogeneous stationary network patterns should also play an important role therein. Additional investigations of more realistic models of coupled place and grid cells are needed. This will allow to connect with experiments and further contribute to the challenge of how noise affects network dynamics in Rowland et al. (2016, Future Issue 3). We carefully checked that the support of the distribution in s is essentially inside the interval [0, 1.3] for all times such that imposing no-flux boundary condition on the right end is a good approximation for s ∈ [0, ∞). The initial distribution is prescribed by randomly setting the activity at one percent of the grid points {x n } n , denoted {x n i } i , in the four sheets to 1, i.e.,  or the maximum time of 6000ms is reached. The respective noise intervals are divided into at least 100 points.

A.1 Grid refinement of the numerical method
To briefly check the robustness of the numerical method described in the manuscript, we performed a standard grid refinement analysis within a computationally feasible computational range from n = 64 to n = 256 cells in x and s. For n = 256, the initial data f β,0 (x)| s=1 is randomly set to 1/256 at 1% of the locations x, and f β,0 (x)| s=0 is set to 1/256 on the complement. The rest is set to zero. On the coarser grids, we have used piecewise averages of this initial data. In Fig. 7, the solutions on a grid consisting of n 3 cubes, with n = 32, 64, 128, and 256, are plotted. The corresponding L 1 and L 2 errors can be found in the table below. The coarser solutions (n = 32, 64, 128) are compared to the one computed on the grid consisting of 256 3 cubes (Table 1). It can be observed from the table and the figure that the numerical method used is stable when refining the grid.