Dynamical mean-field theory for a highly heterogeneous neural population with graded persistent activity of the entorhinal cortex

Abstract

The entorhinal cortex serves as a major gateway connecting the hippocampus and neocortex, playing a pivotal role in episodic memory formation. Neurons in the entorhinal cortex exhibit two notable features associated with temporal information processing: a population-level ability to encode long temporal signals and a single-cell characteristic known as graded-persistent activity, where some neurons maintain activity for extended periods even without external inputs. However, the relationship between these single-cell characteristics and population dynamics has remained unclear, largely due to the absence of a framework to describe the dynamics of neural populations with highly heterogeneous time scales. To address this gap, we extend the dynamical mean field theory, a powerful framework for analyzing large-scale population dynamics, to study the dynamics of heterogeneous neural populations. By proposing an analytically tractable model of graded-persistent activity, we demonstrate that the introduction of graded-persistent neurons shifts the chaos-order phase transition point and expands the network’s dynamical region, a preferable region for temporal information computation. Furthermore, we validate our framework by applying it to a system with heterogeneous adaptation, demonstrating that such heterogeneity can reduce the dynamical regime, contrary to previous simplified approximations. These findings establish a theoretical foundation for understanding the functional advantages of diversity in biological systems and offer insights applicable to a wide range of heterogeneous networks beyond neural populations.

Citation:Tomita F, Teramae J-n (2025) Dynamical mean-field theory for a highly heterogeneous neural population with graded persistent activity of the entorhinal cortex. PLoS Comput Biol 21(9): e1013484. https://doi.org/10.1371/journal.pcbi.1013484

Editor:Daniel Bush, University College London, UNITED KINGDOM OF GREAT BRITAIN AND NORTHERN IRELAND

Received:April 2, 2025; Accepted:September 3, 2025; Published: September 16, 2025

Copyright: © 2025 Tomita, Teramae. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability:The code used to run numerical simulations and analysis is available at https://github.com/fuuta/multiD_hetero_DMFT.

Funding:F. T. was supported by Japan Science and Technology Agency (JST) (https://www.jst.go.jp/EN/), the establishment of university fellowships towards the creation of science technology innovation, Grant Number JPMJFS2123. J. T. is supported by Japan Society for the Promotion of Science (JSPS) (https://www.jsps.go.jp/english/) KAKENHI Grant Number 24K15104, 22K12186, and 23K21352, and Japan Agency for Medical Research and Development (AMED, https://www.amed.go.jp/en/index.html) AMED-CREST 23gm1510005h0003. The funders have no roles in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

The entorhinal cortex, as the major gateway for information entering the hippocampus from various brain regions, plays an essential role in coding long temporal information, a crucial component of episodic memories [1–3]. Episodic memory is the ability of animals to store specific events they have experienced, along with their order and contextual details. The formation of such memories requires the brain to represent temporal information across various time scales. Recent experimental findings indicate that neurons in the lateral entorhinal cortex, along with hippocampal time cells [4–8], offer a fundamental mechanism for representing long temporal information required for episodic memory formation [9]. Experiments in rodents and humans have reported that neurons in the lateral entorhinal cortex exhibit activity characterized by gradually rising or decaying activity across various time scales [7,9]. A recent experiment also identified ‘temporal context cells’ in the entorhinal cortex, which, in contrast to time cells, respond rapidly following stimulus onset and then gradually return to baseline activity over extended temporal timescales, exhibiting a broad range of decay rates [10]. These activities may provide a suitable temporal code for the formation of episodic memory, capturing the different scales of time over which an animal’s experiences occur [9,10].

In addition to the above population activity, the entorhinal cortex is known for a distinctive single-neuron activity, which is also associated with representing long-term information [11]. This activity is termed graded-persistent activity (GPA). The firing rate of an isolated neuron generally decays rapidly in the absence of external input. However, unlike the general response, some isolated entorhinal cortex neurons can sustain their firing activity for several minutes even after external inputs have ended [11]. Moreover, the sustained firing rate can exhibit various graded values, reflecting the input history of the neuron. Because this single-cell property suggests that a subset of neurons in the entorhinal cortex has much longer characteristic time scales compared to other typical neurons, it has been suggested that these neurons may be involved in functions requiring long temporal information, such as working memory [12–18].

Considering the above, one might expect that the existence of neurons with GPA could influence the population dynamics of the entorhinal cortex in encoding long-term temporal information. However, it remains unclear how the partial introduction of the GPA neurons modulates the dynamical properties of the neuronal population. The main cause of this difficulty is the lack of a theory to describe the dynamics of a highly heterogeneous population of neurons [19,20]. The fact that only a subset of entorhinal neurons exhibits extremely long-term graded-persistent activity implies that the intrinsic time scales of each neuron are largely diverse across the population, beyond the range where the heterogeneity can be considered negligible.

To address the problem and reveal how the partial introduction of GPA neurons modulates the population dynamics of neurons, we first propose an analytically tractable model of neurons showing this characteristic activity. Then, we analyze their population dynamics by extending the dynamical mean-field theory (DMFT). DMFT is a comprehensive and efficient framework for analyzing the population dynamics of both biological and artificial recurrent neural networks [19,21–33]. The theory allows us to reduce the generally high-dimensional dynamics of large populations of neurons to an effective low-dimensional equation. This protocol is rigorously justified by a path integral approach. Specifically, the theory provides a theoretical basis for neural computation of temporal information, including temporal information processing in echo state networks and reservoir computing [34–36], by clarifying the onset of chaotic states in these networks.

DMFT enables us to marginalize the heterogeneous connection strengths between neurons into the effective mean field. However, we will see that we cannot simply average out the intrinsic properties of each neuron, including its characteristic time scale, by similarly incorporating them into the mean field. This is because of the difference in their dependency on the network size. To solve the problem, we extend the DMFT framework to networks consisting of heterogeneous neurons. Unlike conventional DMFT, which provides a single mean-field equation, we will obtain a set of mean-field equations reflecting the intrinsic heterogeneity of each neuron. Nevertheless, we will see that this set of equations can provide a single analytical expression to determine the critical coupling strength of the network. Results of the analysis show that the partial introduction of GPA neurons shifts the transition point to extend the dynamical region of the network. We confirm the validity of the theoretical prediction by comparing this with the results of numerical simulations with various network conditions.

To demonstrate the applicability of the approach, we will test the theory using networks with another type of heterogeneity, specifically the heterogeneity of adaptation in each neuron in the network, as discussed in a previous study [27]. The analysis in the previous work, based on the conventional treatment of heterogeneity, predicts that this heterogeneity would move the transition point in a way that expands the dynamical regime of the network. Contrary to this prediction, we will see that this heterogeneity can shift the network to stabilize the steady state, thus shrinking the dynamical region. This tendency is precisely described by the novel approach proposed here.

The organization of the paper is as follows. In the first subsection of the Results, an analytically tractable model of the GPA neuron is provided. Unlike previous models that have revealed underlying possible biological mechanisms of this characteristic activity [14,37–39], the model consists of only two variables. This model can describe the activity of a neuron with and without GPA by modulating a model parameter. Typical dynamics of networks of these neurons are also given in this subsection. DMFT of the heterogeneous network is developed in the next subsection. By deriving the equation determining the transition point, we discuss how adding GPA neurons to the network modulates its dynamics. Finally, we apply the proposed method to a network of neurons with heterogeneous adaptation. Possible extensions and remaining future subjects are examined in the Discussion. Details of the theory are given in the Methods section.

Result

A two-dimensional model for a neuron with and without graded-persistent activity

Single-neuron model.

A subset of neurons in the entorhinal cortex shows graded-persistent activity (GAP) characterized by prolonged sustained spike firings even after the input current to these neurons is terminated [11]. The sustained firing rate can be continuously increased by repeated additions of depolarizing input current to the neuron and, conversely, can be gradually decreased by inductions of hyperpolarizing input. This unique feature is attributed to the single-cell property rather than network effects, as it remains even with the blockage of synaptic connections [11]. Previous studies have suggested the contribution of intracellular calcium density and membrane channels, such as the activation of calcium-dependent nonspecific cation channels, to the single-cell property, and several detailed computational models for this have been proposed [14,37–39].

We use a simple, analytically tractable two-dimensional model to study networks consisting of neurons with and without graded-persistent activity. The model consists of the variable x representing the neural activity and an auxiliary variable a whose time scale can be very slow. The firing rate of the neuron is given by an activation function that is an increasing function of x. The auxiliary variable a may correspond to some slow dynamics of the cell, such as intracellular calcium concentration. Depending on the value of a model parameter, the model can qualitatively describe both neurons with and without the GPA.

The single-cell model is given by

(1)

where the dot indicates the temporal derivative and I(t) is the external input to the neuron. Coefficient γ represents the decay rate of the auxiliary variable a, and β is the feedback strength from the auxiliary variable a to neural activity x. When the value of β is positive, the auxiliary variable a works as positive feedback to the neural activity x. (Note that these values should satisfy because x will diverge otherwise.)

When the decay rate γ is large, i.e., the characteristic time scale of a is small, the single-cell dynamics is effectively described only by the variable x because the external input does not largely increase the value of a. It immediately vanishes with the termination of the input. Thus, the model behaves as a normal neuron without graded persistency, whose activity promptly decays without the external input (Fig 1, left panels). On the contrary, if the decay rate is small (Fig 1, right panels), the value of a is almost kept constant even without external inputs. Due to the finite support of the auxiliary variable, the neural activity can be sustained, avoiding decay to zero, even without external input, which quantitatively reproduces the graded-persistent activity. The sustained activity gradually increases or decreases, reflecting the input history of the neuron, which agrees with experimental findings [11].

[Figure omitted. See PDF.]

Fig 1. Typical dynamics of the model neuron, Eq (1), with and without the graded-persistent activity (GPA) when successive pulse external stimulus (top panels) are given.

For a normal neuron (left panels, ), the cell activity (left middle panel) decays rapidly when the pulse input is terminated as its auxiliary variable (left bottom panel) is kept small. By contrast, for a GPA neuron (right panels, ), the cell activity (right middle panel) is kept almost constant, even during intervals between pulse external inputs. The sustained cell activity gradually increases in response to repeated depolarizing pulse inputs and decreases in response to subsequent hyperpolarizing pulse inputs. These sustained activities are supported by the slow decay of the auxiliary variable (right bottom panel).

https://doi.org/10.1371/journal.pcbi.1013484.g001

Heterogeneous network with a subset of the GPA neurons.

Using the single neuron model, define the heterogeneous network whose subset neurons have the graded persistency

(2)

Here, denotes the index of neurons, N is the number of neurons in the network, and J = (Jij) is the connection weight matrix where Jij represents the synaptic connection strength from the jth to the ith neuron. To avoid self-connection, we set the diagonal components of J as Jii = 0 for all i. The values of the off-diagonal components Jij are independently chosen from a Gaussian distribution with variance g2/N, i.e., . Here, the parameter g controls coupling strength. In all numerical simulations in this paper, we use , while the results of the theoretical analysis are not restricted to this particular choice of the activation function.

To make the neurons in the network heterogeneous, we randomly and independently chose the decay rates from a two-point distribution:

(3)

where . Because a normal neuron is modeled by a high value of the decay rate, the parameter p, determining the ratio of the neurons with the low value of γ, controls the ratio of GPA neurons in the network. For the feedback strength parameter , we use a uniform value unless stated otherwise. (Note that, while we use a two-valued distribution for simplicity, the theory developed here is not limited to discrete parameter distributions. For details, please see the Methods section, S1 Appendix, S1 and S2 Figs.)

Fig 2 shows typical temporal dynamics of the network for values of model parameters. The increase in coupling strength g moves the network from a quiescent state (from the left to the right panels in Fig 2), where the activities promptly decay to zero, to a seemingly chaotic state where irregular neural activities are sustained. This result agrees with previous studies that reported the existence of the chaos-order transition at g = gc, where gc is the transition point [21,28].

[Figure omitted. See PDF.]

Fig 2. Typical dynamics of neurons in the random network, Eq (2), for values of the coupling strength g and the ratio of the GPA neurons p.

Panels are arranged such that the coupling strength g increases from left to right (), and the ratio of the GPA neurons increases from top to bottom (). For the larger value of the ratio of GPA, p, the transition from the silent to the chaotic state occurs at the smaller value of g. Other parameters are N = 3000, , , and .

https://doi.org/10.1371/journal.pcbi.1013484.g002

When the ratio p of graded-persistent neurons in the network is increased (from the top to the bottom panels in Fig 2), the transition from the quiescent state to the chaotic state occurs at smaller values of g. This result suggests that the transition point itself is shifted toward expanding the chaotic regime by introducing GPA neurons to the network. The observation is of particular interest because the results of numerical simulations show that graded persistent neurons and other normal neurons behave almost similarly in the network (Fig 3). In the next section, we will see that this is actually the case and derive an equation determining gc as a function of the model parameters, including the ratio p.

[Figure omitted. See PDF.]

Fig 3. Behavior of graded persistent neurons and other normal neurons in a network, Eq (2), operating in the chaotic regime.

(a) Temporal profiles of the activity of neurons in each population. (b) Autocorrelation functions of the time series for individual neurons. (c) Autocorrelation functions averaged over neurons in each population. (d) Density distributions of the correlation times of neurons in each population, obtained directly from the autocorrelation functions in (b). Parameters N = 5000 p = 0.5 and g = 2.0 are used. All results indicate that graded persistent neurons and other normal neurons behave similarly in the network under chaotic dynamics.

https://doi.org/10.1371/journal.pcbi.1013484.g003

Dynamical Mean Field Theory for the network with highly heterogeneous neurons

To analyze the shift of the transition point induced by the intrinsic heterogeneity of neurons, we use Dynamical Mean Field Theory (DMFT) [19,21,22,24–33]. DMFT mitigates the difficulty of studying the population dynamics of high-dimensional nonlinear systems by allowing us to replace the interaction term in each neuron’s dynamics, i.e., the third term of Eq (2), with an effective dynamical mean field. The properties of this mean-field are determined by the neural dynamics evolving under the mean field, leading to a self-consistent equation for the order parameters that characterize the macroscopic network dynamics.

In the conventional procedure of DMFT, one can obtain a single equation of motion driven by the mean-field to describe the dynamics of all neurons in the network due to their homogeneity. Unlike these conventional approaches, the heterogeneity of neurons cannot be averaged out from the network (See Method for and following paragraphs details). Instead, applying the DMFT framework leads to a set of N differential equations governing neural dynamics, with the heterogeneity preserved. The difference between the two arises from the system size dependence of the terms in the generating functional of mean-field theory: connection heterogeneity appears as a sum of N2 terms, while neuron heterogeneity appears as a sum of N terms; thus, one cannot treat them the same way. However, we will see that these equations can be averaged in Fourier space, resulting in a single equation that characterizes the chaos-order transition of the heterogeneous network.

Following previous works [21,27,28], let us assume that the external input Ii to each neuron is an independent realization of a Gaussian process with zero mean. Then, by adapting the path integral method to the dynamical equations, Eq (2), one can replace the interaction term in each neuron’s dynamics with another Gaussian process (See Method for details):

(4)

where the mean of is zero, and its autocorrelation should be determined self-consistently by the constraint:

(5)

If the network consisted of homogeneous neurons, the above equations would be the same for all neurons, and we could remove the subscript i representing the neural index. However, we must keep the index because the decay rates differ across neurons. Note that a previous study proposed eliminating similar intrinsic heterogeneity by averaging the above equations over all neurons [27]. However, we will demonstrate in the next section that such a naive averaging fails to deliver accurate results.

Whereas we cannot reduce the set of equations Eq (4) to a single one, it is still possible to obtain a single equation determining the transition point of the network from the set of equations (See Method for details). The Fourier transform of the above equations gives

(6)

where , , and are the normalized Fourier transforms of xi, ai, and , respectively, that is, , , and , where . Eliminating Ai from the above and multiplying the resulting equation by its complex conjugate, we have

(7)(8)

where is the power spectral density of the activity of the ith neuron xi, is the average power spectral density of the activity via activation function . One can safely average Eq (7) over all neurons because, unlike Eqs (2) and (4), it does not have higher-order terms of heterogeneity, such as , which would result in additional unknown terms. The average gives

(9)(10)(11)

where in the third line, we have taken the limit of and used the law of large numbers, defining .

Now, let us assume that the activation function satisfies the condition

(12)

which is certainly satisfied by most of the generally used activation functions, including we used. This condition leads to an inequality between the power spectral densities:

(13)

which, when combining with equation Eq (11), gives the required equation determining the transition point gc of the heterogeneous network:

(14)

Introduction of GPA neurons expands the dynamical regime of the network

When the decay rate of each neuron is independently chosen from the two-point distribution, Eq (3), one can obtain an analytical expression of and the transition point gc explicitly:

(15)

If the network consists of homogeneous neurons (i.e., p = 0), this expression simplifies to , which reproduces the conventional result gc = 1 in the limit as infinite decay rate () or no feedback (), where the single neuron model reduces to a conventional one-dimensional dynamical system without the auxiliary slow variable.

If one increases the heterogeneity of the network, the transition point gc generally decreases from one, as gc is a monotonically decreasing function of p since . Therefore, introducing GPA neurons with slow dynamics to the network consistently shifts the network toward a more dynamic regime, i.e., a state preferable for temporal information encoding.

To validate the theoretical prediction and to examine how the shift of the transition point is affected by model parameters, we performed numerical simulations of the network dynamics in Eq (2) and compared them with the theoretical prediction Eq (15) for various values of the model parameters. To numerically identify the transition point, we calculated the maximum value of the power spectrum from the numerically obtained time series. Since the power spectrum is nonnegative by definition, a positive maximum indicates that the network is in a dynamic regime. In contrast, the zero maximum implies the network has converged to a steady state.

Fig 4 shows the results. We can see that the theoretical predictions (red lines) agree well with the results of numerical simulations in all cases, and the dynamical regime consistently expands as the heterogeneity p increases from 0. This expansion of the dynamical regime is particularly pronounced when the decay rate of GPA neurons is small (the far left panel of Fig 4a) and the feedback strength β is large (the far right panel of Fig 4c), where GPA neurons exhibit slower auxiliary dynamics (small ) and have a greater influence on the neural activity x (large β). This is because the slow auxiliary variable of GPA neurons facilitates the nonzero firing activity of other neurons in the network and helps to prevent the population dynamics from converging to a quiescent state.

[Figure omitted. See PDF.]

Fig 4. Maximum power spectrum of the network dynamics for values of the ratio of the GPA neuron p (horizontal axis of each panel) and coupling strength g (vertical axis of each panel).

Red curves are the theoretical prediction of the transition point, gc, given by Eq (15). a: and are fixed and increases from left to right (). b: and are fixed and increases from left to right (). c: and are fixed and β increases from left to right (). Numerical simulations are performed with N = 3000 for all panels.

https://doi.org/10.1371/journal.pcbi.1013484.g004

Analysis of the transition induced by the introduction of GPA neurons

To more precisely examine the transition of the network from the quiescent state to the dynamical state, we numerically calculated the eigenspectrum (Fig 5) and the maximum Lyapunov exponent (Fig 6) of the network for values of the ratio of GPA neurons p and the coupling strength g.

[Figure omitted. See PDF.]

Fig 5. Eigen spectrum in the complex plane for the Jacobian at the trivial fixed point of the network dynamics.

Dots in each panel indicate the eigenvalue of the linear stability matrix on the complex plane. Blue and red dots indicate eigenvalues with negative and positive real parts, respectively. Panels are arranged such that the GPA neuron ratio, p, increases from left to right, and coupling strength, g, increases from top to bottom. The thick vertical line in each panel is the imaginary axis. Panels shaded by red background mean the coupling strength of the network is above the transition point predicted theoretically, i.e., g > gc. Other model parameters are and .

https://doi.org/10.1371/journal.pcbi.1013484.g005

[Figure omitted. See PDF.]

Fig 6. The maximum Lyapunov exponent of the network dynamics.

a: The numerically estimated maximum Lyapunov exponent, averaged over ten realizations of random initial conditions, for the model network as a function of the coupling strength g. Line colors indicate the value of p (from the far-right purple curve for p = 0.0 to the far-left yellow curve for p = 1.0). Dashed vertical lines indicate gc. The number of neurons in the network is fixed at N = 3000. b: The same Lyapunov exponent for different values of N with p = 0.8. Blue, orange, green, and red lines correspond to networks of N = 100, 500, 1000, and 3000, respectively. Error bars show standard deviation.

https://doi.org/10.1371/journal.pcbi.1013484.g006

The eigenspectrum (Fig 5) of a fixed point in the network characterizes the system’s linear stability around that point. For recurrent neural networks with random connections, the instability of the trivial fixed point often corresponds to the chaos-order transition point, as the loss of this stability may lead to chaotic behavior due to the network’s intrinsic randomness [21]. We linearized the model system, Eq (2), around its trivial fixed point, or the origin, and numerically computed the eigenspectrum of the Jacobian matrix. The panels in Fig 5 are arranged such that the coupling strength, g, increases from top to bottom, and the ratio of the GPA neuron, p, increases from left to right. The trivial fixed point is stable if the real parts of all eigenvalues are negative and unstable if at least one has a positive real part. The eigenvalue distributions exhibit complex, nontrivial shapes on the complex plane. However, it is evident that instability, indicated by the appearance of a positive eigenvalue, occurs at smaller values of g for larger values of p, as predicted by the theoretical results.

The maximum Lyapunov exponent is used to characterize the complexity of nonlinear systems, particularly their chaotic behavior, as a positive exponent indicates the network’s chaotic state. To compute the maximum Lyapunov exponent, we employed the method proposed by Benettin et al. [40]. This involves numerically solving the model equations, Eq (2), to obtain the trajectory and estimation of the exponent from the growth rate of the tangent vector along the trajectory. Fig 6a shows that the maximum Lyapunov exponent changes sign from negative to positive at points close to the theoretical predictions. While small discrepancies exist between the numerical simulations and theoretical results, we confirmed that these mismatches arise from finite-size effects. As the number of neurons in the network increases, the discrepancies decrease (Fig 6b), which is consistent with the fact that DMFT assumes an infinite number of neurons.

Gaussian heterogeneity of adaptation of neurons in random networks

To demonstrate the applicability of the developed theory, we apply it to a different type of neural heterogeneity: Gaussian-distributed adaptation of each neuron in the network. In a pioneering theoretical work studying how intrinsic properties of single neurons modulate population dynamics in random neural networks, Muscinelli et al. proposed a two-dimensional neuron model with strong adaptation [27]. Remarkably, this adaptation neuron model is equivalent to the GPA neuron model introduced here but with the opposite sign of a in the equation for neural activity. Specifically, the sign of the second term in Eqs (1) and (2) is negative, rather than positive, in the adaptation neuron model. Due to the equivalence, it is straightforward to apply our method to a heterogeneous network of the adaptation neuron model.

In the previous study, the authors introduced heterogeneity in adaptation for each neuron by randomly choosing the feedback strength parameter from a Gaussian distribution with mean and small variance , i.e., , using the variance to characterize the heterogeneity. (The variance was kept small to ensure that the sampled values of remained negative.) Unlike the theory developed here, the authors treated the heterogeneity similarly to the coupling strengths between neurons. They naively assumed that the heterogeneity could be effectively expressed by an additional Gaussian mean field in the DMFT. Based on this assumption and following the conventional DMFT procedure, they derived a single mean-field equation driven by two Gaussian processes, which gives the equation determining the transition point of the heterogeneous network:

(16)(17)

where the function G is given by Eq (8). To emphasize the difference between this transition point and the one obtained in the next section, we denoted the transition point obtained here with hat as .

However, the assumption of the previous study that the intrinsic properties of neurons can be expressed by an additional Gaussian process was not entirely correct. (Technically speaking, averaging the generating functional of the model equation is still possible even with this type of heterogeneity. However, the equation cannot be reduced to a single expression. This is because the system-size dependence of the term arising from heterogeneity scales differently from that of the connection heterogeneity (see Method for details)). Instead, we need to deal with a set of N differential equations, which must be averaged in Fourier space to determine the transition point of the network, as described in the previous section. Following the same procedure as in the previous section, we have the equation for the transition point:

(18)

where the function is given by Eq (8) and represents the average of the function over sampled values of in the network.

Interestingly, the equations for the transition point derived in the previous study and ours can predict opposite tendencies regarding how this heterogeneity shifts the chaos-order transition point in certain cases. The equation of the previous study predicts that the transition point decreases with increasing . Oppositely, the equation derived here predicts that the transition point will increase, meaning the dynamical regime will shrink due to this type of heterogeneity. To test these predictions, we numerically simulated the population dynamics for various model parameters. Fig 7 shows that, as predicted by the current theory, the transition point increases as increases, whereas shows the opposite trend. Further exploration of the model parameters causing this discrepancy between the two theories is an important future subject.

[Figure omitted. See PDF.]

Fig 7. Maximum power spectrum of the network dynamics with Gaussian heterogeneous adaptation for values of model parameters.

The horizontal and vertical axes are and g, respectively. The red curves show theoretical prediction gc obtained using the method developed in this study The white dashed lines indicate that are derived from a naïve approximation of the heterogeneity.

https://doi.org/10.1371/journal.pcbi.1013484.g007

Discussion

In this study, considering the properties of the entorhinal cortex, we develop a theory to describe the population dynamics of random neural networks consisting of highly heterogeneous neurons. We propose a simple two-dimensional neuron model representing both normal neurons and neurons exhibiting graded persistent activity (GPA). We extend the well-established theoretical tool, DMFT, which provides an effective mean-field equation for network dynamics, to cases where the neurons in the network are highly heterogeneous. Unlike conventional DMFT, the derived mean-field model consists of N-dimensional stochastic equations rather than a single equation. However, we showed that averaging these equations is still possible in Fourier space. This allows us to theoretically determine the transition point of the network by focusing on the average power spectrum of the network dynamics. Since recent experiments have revealed that cortical neurons are highly heterogeneous in their intrinsic properties, the theory developed here will be an important tool, alongside conventional DMFT, for clarifying the properties of recurrent networks of realistic neuron models.

The results of the theory suggest that the stable region of the network shrinks, and the network becomes more dynamic as the ratio of GPA neurons increases, regardless of other model parameters. Since network activity triggered by external inputs does not decay near the boundary between the stable and dynamical regimes, lowering the coupling strengths required to reach this boundary may facilitate the encoding of long-timescale information in the entorhinal cortex. This boundary is widely known as the “edge of chaos” in the context of reservoir computing and is associated with optimal computational capabilities.

To validate this concept within our framework, we numerically calculated the memory capacity, a common metric in reservoir computing [34] (See S1 Appendix for details). As demonstrated S4 Fig, the memory capacity is indeed maximized at the onset of the chaotic regime, just below the critical transition point predicted by our theory. This finding not only supports the computational advantages of operating at the edge of chaos but also implies a novel tuning mechanism: when the available range of synaptic coupling strengths is limited, a network may still be able to achieve its maximal memory capacity by modulating its intrinsic heterogeneity (See S1 Appendix for details).

Of particular interest is the relationship between our study and the computational model of temporal representation proposed by Shankar and Howard [41], which employs a set of leaky integrators with a spectrum of decay timescales to encode stimulus history. Our work and theirs address distinct, though related, questions. While our study investigates how the existence of neurons with diverse timescales impacts the collective dynamics of a network, Shankar and Howard aimed to answer the more specific computational question of how these diverse timescales can be utilized to construct a scale-invariant representation of past events. This divergence in research goals is reflected in their respective methodologies. In contrast to our use of an extended DMFT to analyze a random network’s statistical properties, Shankar and Howard proposed a specifically structured, layered network architecture designed to approximate an inverse Laplace transform of the stimulus history. Consequently, the functional role of heterogeneity is also interpreted differently. In our framework, the diversity of neural timescales acts as a modulator that alters the network’s global dynamical regime, whereas in their model, this same diversity is the essential substrate for building the temporal representation itself. Despite these differences, the two frameworks are not contradictory but are better viewed as complementary. The population of leaky integrators posited by Shankar and Howard could, for instance, be biologically realized by the heterogeneous neural population that includes GPA neurons, as modeled in our study. An integration of these two perspectives could therefore provide a more comprehensive understanding that bridges the gap from the role of diversity in single elements and population dynamics to the emergence of higher-order cognitive functions like temporal representation.

The networks considered here belong to a class of systems that maintain memory of past inputs through transient trajectories, rather than storing information in stable attractors like traditional Hopfield networks. The capacity for such temporal coding is thought to depend critically on non-normal dynamics, which can transiently amplify specific inputs [42,43]. Our model, based on random recurrent interactions, naturally gives rise to this non-normal property. This raises the intriguing possibility that by tuning the nature of the heterogeneity, for instance, by altering the specific variables or the shape of their statistical distributions, our framework could account for specific, dynamic coding phenomena. A compelling target for such investigation would be the transient amplification of odor representations observed experimentally in the insect antennal lobe, a phenomenon that has also been linked to non-normal network dynamics [44,45].

Experimental findings suggest that characteristic neural dynamics encoding time arises from the interaction between the network’s internal dynamics and external signals to the network. In general, when external signals are applied to a nonlinear system, they tend to entrain the system and increase its stability [22,46–54]. Therefore, if the internally generated intrinsic dynamics of the network are in a stable and quiescent state, the system struggles to integrate external signals. This is because the signals further stabilize the network, causing input information to decay too quickly to be retained.

In contrast, if the network operates in a dynamic regime near the transition point, external signals can entrain the network and drive its dynamics closer to or just below the transition point. This is the ideal state for integrating external inputs and maintaining their information over a long time. Thus, an important direction for future research will be to extend the theory to networks influenced by external signals. Additionally, it is also an important future task to investigate the relationship between the network’s intrinsic dynamics and the characteristic population activity of the entorhinal cortex, such as the activity of time cells and the ramping activity. In this regard, we have demonstrated that neurons in our heterogeneous recurrent network can reproduce response patterns similar to those of ’temporal context cells’ recently observed in the entorhinal cortex, when intrinsic time scales, in particular the feedback strengths , follow continuous distributions (see S1 Appendix and S3 Fig for details).

We also applied the theory to a different type of heterogeneity studied in previous work: neuronal adaptation. The theory developed here accurately describes the numerical results, demonstrating that heterogeneity in adaptation can reduce the dynamical regime of the network. However, we show that the predictions of the previous study do not always capture the direction of the shift in the transition point. This result highlights that intrinsic heterogeneity should be treated carefully and differently from the heterogeneity of coupling strengths. Another important message from this result comes from the comparison between Figs 4 and 7. Although both results show that introducing heterogeneity to the intrinsic properties of neurons modulates population dynamics and shifts the transition point, the directions of the shifts are entirely opposite, indicating that one cannot generally conclude whether heterogeneity expands or shrinks the dynamical regime of the network without knowing the specific details of the heterogeneity.

The chaos-order transition studied here has been recognized as playing an important role in temporal information processing, especially within the framework known as reservoir computing, where input sequences are stored in randomly connected recurrent networks. For instance, several studies have shown that the memory capacity of a neural population for input sequences is maximized near the transition point [55]. Therefore, an interesting future direction would be to evaluate the memory capacity of heterogeneous networks of GPA neurons by extending the theory developed here.

This study is closely related to recent important work by Stern et al. [19]. In their study, Stern et al. demonstrated that the heterogeneous distribution of cell assembly sizes can explain the experimentally observed heterogeneous time scales in neuronal circuits. By extending the DMFT, they showed that the heterogeneity of time scales gives rise to a novel chaotic regime characterized by bistable activity. To reveal this chaotic regime, their analysis employed random matrix theory, which can provide the eigenspectrum distribution of large random matrices [56]. In contrast, our approach demonstrates that the transition point of a heterogeneous network can be explicitly determined by a single equation, Eq (14), derived by directly averaging the heterogeneity in Fourier space. Given the difficulty in deriving the eigenspectrum distribution for most matrix ensembles, our results complement their work by offering an alternative framework to elucidate the significant role of intrinsic heterogeneity in large-scale systems. Further studies of the transition mechanism could deepen the connections between their theoretical advancements and ours.

Related to the above, another promising theoretical direction is the study of the eigenspectrum distribution of the model. As shown in Fig 5, the eigenvalues of the model network exhibit a highly intricate yet well-organized nontrivial structure. Several bulks or clusters of eigenvalues are observed, which gradually merge as the model parameters vary. Since the connection matrix of the network is randomly generated, such eigenspectrum distributions may be analyzed within the framework of random matrix theory [31,56–61]. Given the critical role of heterogeneity in a wide range of applications, understanding the mechanisms underlying the modulation of the eigenspectrum could provide deeper insights into the dynamical behavior of large, generally heterogeneous systems.

Methods

Dynamical mean field theory for the neural networks with heterogeneous parameters

In this section, within the framework of dynamical mean-field theory (DMFT) [19,21,22,24–33], we derive effective dynamics (Eq (4) in the main text) of each neuron from the original dynamics, i.e., Eq (2) in the main text:

(19)

By formally adding white Gaussian noise to each term in the above equations and defining the right-hand sides as and , we obtain a set of stochastic equations:

(20)

Considering this as the Ito stochastic equation and discretizing it in time with interval yields

(21)

where , , , and . Since Eq (21) give two-point relationships of variables, we can obtain the high-dimensional joint probability density function of marginalized over from this as

(22)(23)

where we have shortened the expression on the left-hand side of Eq (21) by defining and . In the second line, we used the Fourier representation of the Dirac delta function, and in the third line, we introduced the moment generating function of the stochastic variable , defined by .

Let us consider the behavior of the stochastic distribution, Eq (23), in the limit of , or . The first and the second terms of the exponential function converge as

(24)

and

(25)

where we used the moment generating function of Gaussian distribution and we denoted the variance of the noise as . By putting them together, we have

(26)

Eliminating the formally introduced Gaussian noise by putting , we obtain

(27)

where, , , , , and . Here, and B are the diagonal matrices whose ith diagonal component is given by and , respectively. Bold lowercase letters mean vectors whose inner product shall be simultaneously in the spatial and temporal directions, i.e., .

The moment generating functional of the stochastic process given by the probability distribution Eq (27) is defined as

(28)

where and are the auxiliary variables corresponding to x(t) and a(t), respectively. Other auxiliary variables and are introduced to represent arbitrariness of their initial values [28]. We also introduced vector notations for integrals, and .

Then, by putting Eq (27) to Eq (28), we have

(29)(30)

Since we want to know the network’s behavior independent of each realization of the coupling matrix J, let us average the moment generating functional over the distribution of the coupling matrix (quenched average). The coupling strengths independently follow the same normal distribution with mean 0 and variance g2/N, and J appears in Eq (29) in the form of . So, we can average this as:

(31)(32)

where we used the identity in the second line. In the last line, we omitted the diagonal term, i.e., the second term of the exponential function, because this term scales in the order of N, which is sufficiently small compared to the off-diagonal first term, which scales with N2 in the limit of . (Note that this difference in system size dependence will be important when discussing each neuron’s intrinsic heterogeneity.) Using Eq (32), we obtain

(33)

where, for simplicity, we assumed that there is no external input, i.e., I = 0. Note that each second-order term of in the exponential of the above equation, i.e., the term of , shares the common factor , which is independent of the subscript i. As we will see below, this independence from the neuron’s index is key to deriving a single effective equation.

To proceed, let us define the auxiliary field Q1:

(34)

and represent this relationship by using the delta functional with an additional auxiliary field Q2:

(35)

Putting this delta functional into Eq (33) allows us to treat Q1 as an independent variable in the integral. Then, by using the Fourier representation of the delta functional, we have

(36)

We can see that interactions among neurons are replaced by the effective interaction between each neuron and the auxiliary field. This result allows us to decouple neural dynamics, except for the effective interaction via the auxiliary field.

To have the explicit expression of the decoupled neural dynamics, formally rewrite the above expression by removing vector response terms , , , and and introducing two scalar variables and to incorporate terms for the auxiliary field into the generating functional [28]:

(37)

where

(38)

Here, we introduced notations and

.

The exponent of the exponential function of Eq (37) increases with the order of N as we increase the number of neurons in the network. Using this property, one can evaluate the above integral by using the saddle point approximation that allows us to replace the integral of the exponential function with the exponential function itself of the maximum values for its variables and . The maximum condition, i.e., the saddle-point equation, is given by

(39)

where denotes the derivative of a functional in Q1 and Q2. Because these saddle-point equations are explicitly written as

we can solve them as

(40)

where means the average with respect to the probability distribution given by the solution of the saddle-point equation. Then, using the solution, we have the final expression of the averaged moment generating functional:

(41)

The final expression is same to the moment generating functional with noise and no coupling in Eq (26). Therefore it represents that each variable xi is commonly driven by the Gaussian process with mean 0 and the autocorrelation . Thus, pulling back the expression to differential equations, we have the effective dynamics of the population of neurons driven by the Gaussian mean field:

(42)

which is equal to Eq (4) of the main text.

What happens if we naively average out intrinsic heterogeneity of neurons

In the previous section, we derived Eq (33) by averaging the heterogeneity in the coupling strengths . However, we left the intrinsic heterogeneity, namely the variability in each neuron’s parameter , untouched. In this subsection, we will illustrate how the analysis becomes problematic if we average the intrinsic heterogeneity in the same way as the coupling strengths.

Let us revisit Eqs (29) and (30), which provide the moment generating functional of the stochastic process before the averaging. Assuming that each neuron’s parameter independently follows a Gaussian distribution , we have the functional averaged over the intrinsic parameter as:

(43)

where is the functional given by Eq (28). Since is in only Sxa of the functional, as described in Eq (30), performing the integral requires evaluation of

(44)

where is equivalent to Z except Sxa in it is replaced by . We can perform the integral and obtain

(45)

where we used the identity in the second line and the notation .

The equation above corresponds to Eq (31) in the previous subsection. In that subsection, one could omit the order N term because the leading order of the exponent in the exponential function was N2. However, the order N term cannot be neglected here because no higher-order terms are present. Consequently, instead of , the factor determines the variance of the effective Gaussian noise in the decoupled equation of neural dynamics. Unlike , however, this factor depends on the neuron index i, making the resulting one-body equation heterogeneous. As a result, a single effective equation cannot be derived, and the effective equation must be given as a set of equations, still with the heterogeneity of each neuron. (Note that a single equation could be obtained if we were to forcibly replace with , ignoring the i-dependence. However, as shown in Sect 3, this approach fails to produce accurate predictions.

Derivation of the relationship between power spectrums

In this subsection, we explain details of the derivation of the Eq (11) of the main text.

Let us define the Fourier transform and the normalized Fourier transform by

(46)(47)

Then, the effective network dynamics, Eq (4) of the main text, without external input

(48)

are transformed to

(49)

From the above expression, we have

(50)

where we have omitted ω from the expression for simplicity. Substituting this to the first line of Eq (49) gives

(51)

Then, by multiplying complex conjugates of themselves to both sides of this expression, we have

(52)

where, in the last line, we used and that directly followed from the definition of the power-spectrum density.

For the right-hand side of the above equation, we have from the Wiener-Khinchin theorem, and, using them, we can show

(53)

with defining and .

Substituting Eq (53) to Eq (52) and averaging both hand sides of the equation over the index i gives the desired relationship between power-spectrums:

(54)(55)

In the last line, we replace the arithmetic average of G with its mean, that is, .

Derivation of the equation for the transition point of the heterogeneous network

In this subsection, we will derive the equation, Eq (14), that determines the transition point gc of the network from the power-spectrum equation derived in the above subsection, Eq (55), or Eq (11) of the main text.

Assume that the activation function satisfies the condition , which is the condition that is satisfied by most of the standard activation functions, including and ReLU. Then, it directly follows that

(56)(57)

Because Parseval’s theorem of the power spectrum gives

(58)

Substituting them into Eq (57) and taking the mean for both hand sides of the equation gives

(59)(60)(61)

Then, by substituting Eq (61) to Eq (55), we arrive

(62)

Note that Eq (62) is an absolute inequality that must be satisfied by any values of , g, and x, as far as the activation function satisfies the above-introduced condition . Now, let assume that the coupling strength g satisfies an inequality . Then, one can immediately see that must satisfies to fulfill the absolute inequality Eq (62). Because the power spectrum for all ω means the neural activity is , we can conclude the network must be in the silent state when . On the other hand, if the coupling strength satisfies , we can say that the network is allowed to be in an active state, while we cannot say that the network must be in an active state. Thus, we can conclude that the condition

(63)

gives the critical coupling strength. Note that we can rewrite this condition as because takes its maximum value at under the usual condition of where the activity of the GPA neurons does not diverse.

Analytical expression of the critical coupling strength for the two-point distribution of

In this subsection, we will give the analytical expression of the critical coupling strength gc when the decay rate follows the two-point distribution:

(64)

Let us assume that the model parameters fulfill the condition to each single neural activity does not diverge. by differentiating Eq (8) of the main text, we have

(65)

Under the condition of , has a unique real solution at . Therefore, since and , it follows that provides the maximum value of G given by .

Then, as is just an average of over the given distribution of γ and gives the maximum of G regardless of values of γ, we can conclude that also takes its maximum value at . Therefore, combining this result with Eq (14) of the main text, we have the explicit expression of the transition point gc as follows:

(66)(67)

Supporting information

S1 Appendix. Impact of continuous heterogeneity in intrinsic parameters on network dynamics and computational capacity.

https://doi.org/10.1371/journal.pcbi.1013484.s001

(PDF)

S1 Fig. Network dynamics and temporal properties of neuron subgroups in a network with different decay rates.

(a) Temporal profiles of the activity of neurons with small (left) and large (right) decay rate parameters, corresponding to slower and faster time scales, respectively. (b) Autocorrelation functions of the same neuron groups are shown in (a). (c) Average autocorrelation functions within subgroups of neurons with ranges of decay rate parameters. (d) Distributions of relaxation times computed from the autocorrelation functions of individual neurons prior to averaging. Decay rates follow a uniform distribution U[1,10] across the network. Parameters N = 5000 and g = 2.0 are used.

https://doi.org/10.1371/journal.pcbi.1013484.s002

(TIF)

S2 Fig. Maximum power spectrum of network dynamics with continuous heterogeneity.

The left panel shows the result for a network in which the decay rates follow a uniform distribution which center is 5. The right panel shows the result for a network in which the feedback strengths follow a truncated normal distribution with mean 0 and range [–2,2]. Red lines indicate the theoretical predictions.

https://doi.org/10.1371/journal.pcbi.1013484.s003

(TIF)

S3 Fig. Amplitudes of neuronal responses to a transient impulse input in a heterogeneous network.

The upper heatmap shows the response amplitudes of individual neurons to the same pulse input applied to the network. Neurons are sorted from top to bottom in ascending order of their feedback strength . The lower panel shows the external input, which is shared by all neurons. The network operates in a stable regime where g = 0.15 < gc = 0.183, and the feedback strengths follow a uniform distribution U[0,2.9]. Other parameters are N = 3000 and .

https://doi.org/10.1371/journal.pcbi.1013484.s004

(TIF)

S4 Fig. Memory capacity (MC) of heterogeneous networks for various levels of heterogeneity.

Vertical dashed lines indicate the critical coupling strengths for each network. The input signal u(t) is modeled as Gaussian white noise, and the time evolution of the network is computed using the Euler–Maruyama method for stochastic differential equations. The other parameters are N = 3000, , , , and dt = 0.02.

https://doi.org/10.1371/journal.pcbi.1013484.s005

(TIF)

References

1. 1. Squire LR, Wixted JT. The cognitive neuroscience of human memory since H.M.. Annu Rev Neurosci. 2011;34(1):259–88.

* View Article

* Google Scholar

2. 2. Garcia AD, Buffalo EA. Anatomy and function of the primate entorhinal cortex. Annu Rev Vis Sci. 2020;6:411–32. pmid:32580662

* View Article

* PubMed/NCBI

* Google Scholar

3. 3. Persson BM, Ambrozova V, Duncan S, Wood ER, O’Connor AR, Ainge JA. Lateral entorhinal cortex lesions impair odor-context associative memory in male rats. J Neurosci Res. 2022;100(4):1030–46. pmid:35187710

* View Article

* PubMed/NCBI

* Google Scholar

4. 4. MacDonald CJ, Lepage KQ, Eden UT, Eichenbaum H. Hippocampal “time cells” bridge the gap in memory for discontiguous events. Neuron. 2011;71(4):737–49. pmid:21867888

* View Article

* PubMed/NCBI

* Google Scholar

5. 5. Eichenbaum H. Time cells in the hippocampus: a new dimension for mapping memories. Nat Rev Neurosci. 2014;15(11):732–44.

* View Article

* Google Scholar

6. 6. Sugar J, Moser M-B. Episodic memory: neuronal codes for what, where, and when. Hippocampus. 2019;29(12):1190–205. pmid:31334573

* View Article

* PubMed/NCBI

* Google Scholar

7. 7. Umbach G, Kantak P, Jacobs J, Kahana M, Pfeiffer BE, Sperling M, et al. Time cells in the human hippocampus and entorhinal cortex support episodic memory. Proc Natl Acad Sci U S A. 2020;117(45):28463–74. pmid:33109718

* View Article

* PubMed/NCBI

* Google Scholar

8. 8. Lin D, Huang AZ, Richards BA. Temporal encoding in deep reinforcement learning agents. Sci Rep. 2023;13(1):22335. pmid:38102369

* View Article

* PubMed/NCBI

* Google Scholar

9. 9. Tsao A, Sugar J, Lu L, Wang C, Knierim JJ, Moser M-B, et al. Integrating time from experience in the lateral entorhinal cortex. Nature. 2018;561(7721):57–62. pmid:30158699

* View Article

* PubMed/NCBI

* Google Scholar

10. 10. Bright IM, Meister MLR, Cruzado NA, Tiganj Z, Buffalo EA, Howard MW. A temporal record of the past with a spectrum of time constants in the monkey entorhinal cortex. Proc Natl Acad Sci U S A. 2020;117(33):20274–83. pmid:32747574

* View Article

* PubMed/NCBI

* Google Scholar

11. 11. Egorov AV, Hamam BN, Fransén E, Hasselmo ME, Alonso AA. Graded persistent activity in entorhinal cortex neurons. Nature. 2002;420(6912):173–8. pmid:12432392

* View Article

* PubMed/NCBI

* Google Scholar

12. 12. Frank LM, Brown EN. Persistent activity and memory in the entorhinal cortex. Trends Neurosci. 2003;26(8):400–1. pmid:12900167

* View Article

* PubMed/NCBI

* Google Scholar

13. 13. Okamoto H, Tsuboshita Y, Fukai T. New horizons for associative memory broaden by neurophysiological and computational findings of graded persistent activity. Brain & Neural Networks. 2005;12(4):235–48.

* View Article

* Google Scholar

14. 14. Teramae J-N, Fukai T. A cellular mechanism for graded persistent activity in a model neuron and its implications in working memory. J Comput Neurosci. 2005;18(1):105–21. pmid:15789172

* View Article

* PubMed/NCBI

* Google Scholar

15. 15. Hasselmo ME. Grid cell mechanisms and function: contributions of entorhinal persistent spiking and phase resetting. Hippocampus. 2008;18(12):1213–29. pmid:19021258

* View Article

* PubMed/NCBI

* Google Scholar

16. 16. Yoshida M, Hasselmo ME. Persistent firing supported by an intrinsic cellular mechanism in a component of the head direction system. J Neurosci. 2009;29(15):4945–52. pmid:19369563

* View Article

* PubMed/NCBI

* Google Scholar

17. 17. Koyluoglu OO, Pertzov Y, Manohar S, Husain M, Fiete IR. Fundamental bound on the persistence and capacity of short-term memory stored as graded persistent activity. Elife. 2017;6:e22225. pmid:28879851

* View Article

* PubMed/NCBI

* Google Scholar

18. 18. Ebato Y, Nobukawa S, Sakemi Y, Nishimura H, Kanamaru T, Sviridova N, et al. Impact of time-history terms on reservoir dynamics and prediction accuracy in echo state networks. Sci Rep. 2024;14(1):8631. pmid:38622178

* View Article

* PubMed/NCBI

* Google Scholar

19. 19. Stern M, Istrate N, Mazzucato L. A reservoir of timescales emerges in recurrent circuits with heterogeneous neural assemblies. Elife. 2023;12:e86552. pmid:38084779

* View Article

* PubMed/NCBI

* Google Scholar

20. 20. Park JI, Lee DS, Lee SH, Park HJ. Incorporating heterogeneous interactions for ecological biodiversity. arXiv preprint 2024. https://arxiv.org/abs/physicsbio-ph

21. 21. Sompolinsky H, Crisanti A, Sommers H. Chaos in random neural networks. Phys Rev Lett. 1988;61(3):259–62. pmid:10039285

* View Article

* PubMed/NCBI

* Google Scholar

22. 22. Molgedey L, Schuchhardt J, Schuster H. Suppressing chaos in neural networks by noise. Phys Rev Lett. 1992;69(26):3717–9. pmid:10046895

* View Article

* PubMed/NCBI

* Google Scholar

23. 23. Aljadeff J, Stern M, Sharpee T. Transition to chaos in random networks with cell-type-specific connectivity. Phys Rev Lett. 2015;114(8):088101. pmid:25768781

* View Article

* PubMed/NCBI

* Google Scholar

24. 24. Harish O, Hansel D. Asynchronous rate chaos in spiking neuronal circuits. PLoS Comput Biol. 2015;11(7):e1004266.

* View Article

* Google Scholar

25. 25. Kadmon J, Sompolinsky H. Transition to Chaos in random neuronal networks. Physical Review X. 2015;5(4):041030.

* View Article

* Google Scholar

26. 26. Schuecker J, Goedeke S, Helias M. Optimal sequence memory in driven random networks. Physical Review X. 2018;8(4):041029.

* View Article

* Google Scholar

27. 27. Muscinelli SP, Gerstner W, Schwalger T. How single neuron properties shape chaotic dynamics and signal transmission in random neural networks. PLoS Comput Biol. 2019;15(6):e1007122. pmid:31181063

* View Article

* PubMed/NCBI

* Google Scholar

28. 28. Helias M, Dahmen D. Statistical field theory for neural networks. 1st ed. Cham: Springer; 2020.

29. 29. Keup C, Kühn T, Dahmen D, Helias M. Transient chaotic dimensionality expansion by recurrent networks. Physical Review X. 2021;11(2):021064.

* View Article

* Google Scholar

30. 30. Bordelon B, Pehlevan C. Self-consistent dynamical field theory of kernel evolution in wide neural networks. arXiv preprint 32240–56.

* View Article

* Google Scholar

31. 31. Krishnamurthy K, Can T, Schwab DJ. Theory of gating in recurrent neural networks. Phys Rev X. 2022;12(1):011011. pmid:36545030

* View Article

* PubMed/NCBI

* Google Scholar

32. 32. Pereira-Obilinovic U, Aljadeff J, Brunel N. Forgetting leads to chaos in attractor networks. Phys Rev X. 2023;13(1):011009.

* View Article

* Google Scholar

33. 33. Engelken R, Wolf F, Abbott LF. Lyapunov spectra of chaotic recurrent neural networks. Phys Rev Res. 2023;5(4):043044.

* View Article

* Google Scholar

34. 34. Jaeger H. The “echo state” approach to analysing and training recurrent neural networks-with an erratum note. 148. Bonn, Germany: German National Research Center for Information Technology; 2001.

35. 35. Maass W, Natschläger T, Markram H. Real-time computing without stable states: a new framework for neural computation based on perturbations. Neural Comput. 2002;14(11):2531–60. pmid:12433288

* View Article

* PubMed/NCBI

* Google Scholar

36. 36. Lukoševičius M, Jaeger H. Reservoir computing approaches to recurrent neural network training. Computer Science Review. 2009;3(3):127–49.

* View Article

* Google Scholar

37. 37. Brody CD, Romo R, Kepecs A. Basic mechanisms for graded persistent activity: discrete attractors, continuous attractors, and dynamic representations. Curr Opin Neurobiol. 2003;13(2):204–11. pmid:12744975

* View Article

* PubMed/NCBI

* Google Scholar

38. 38. Loewenstein Y, Sompolinsky H. Temporal integration by calcium dynamics in a model neuron. Nat Neurosci. 2003;6(9):961–7. pmid:12937421

* View Article

* PubMed/NCBI

* Google Scholar

39. 39. Fransén E, Tahvildari B, Egorov AV, Hasselmo ME, Alonso AA. Mechanism of graded persistent cellular activity of entorhinal cortex layer v neurons. Neuron. 2006;49(5):735–46. pmid:16504948

* View Article

* PubMed/NCBI

* Google Scholar

40. 40. Benettin G, Galgani L, Giorgilli A, Strelcyn JM. Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 2: Numerical application. Meccanica. 1980;15(1):21–30.

* View Article

* Google Scholar

41. 41. Shankar KH, Howard MW. A scale-invariant internal representation of time. Neural Comput. 2012;24(1):134–93.

* View Article

* Google Scholar

42. 42. Ganguli S, Huh D, Sompolinsky H. Memory traces in dynamical systems. Proc Natl Acad Sci U S A. 2008;105(48):18970–5. pmid:19020074

* View Article

* PubMed/NCBI

* Google Scholar

43. 43. Goldman MS. Memory without feedback in a neural network. Neuron. 2009;61(4):621–34. pmid:19249281

* View Article

* PubMed/NCBI

* Google Scholar

44. 44. Mazor O, Laurent G. Transient dynamics versus fixed points in odor representations by locust antennal lobe projection neurons. Neuron. 2005;48(4):661–73. pmid:16301181

* View Article

* PubMed/NCBI

* Google Scholar

45. 45. Bondanelli G, Ostojic S. Coding with transient trajectories in recurrent neural networks. PLoS Comput Biol. 2020;16(2):e1007655. pmid:32053594

* View Article

* PubMed/NCBI

* Google Scholar

46. 46. Pikovsky AS, Rosenblum MG, Osipov GV, Kurths J. Phase synchronization of chaotic oscillators by external driving. Physica D. 1997;104(3–4):219–38.

* View Article

* Google Scholar

47. 47. Rosenblum M, Pikovsky A, Kurths J. Phase synchronization in driven and coupled chaotic oscillators. IEEE Transactions on Circuits and Systems I-regular Papers. 1997;44(10):874–81.

* View Article

* Google Scholar

48. 48. Toral R, Mirasso CR, Hernandez-Garcia E, Piro O. Analytical and numerical studies of noise-induced synchronization of chaotic systems. Chaos. 2001;11(3):665–73. pmid:12779505

* View Article

* PubMed/NCBI

* Google Scholar

49. 49. Teramae J-N, Tanaka D. Robustness of the noise-induced phase synchronization in a general class of limit cycle oscillators. Phys Rev Lett. 2004;93(20):204103. pmid:15600929

* View Article

* PubMed/NCBI

* Google Scholar

50. 50. Nakao H, Arai K, Kawamura Y. Noise-induced synchronization and clustering in ensembles of uncoupled limit-cycle oscillators. Phys Rev Lett. 2007;98(18):184101. pmid:17501578

* View Article

* PubMed/NCBI

* Google Scholar

51. 51. Sussillo D, Abbott LF. Generating coherent patterns of activity from chaotic neural networks. Neuron. 2009;63(4):544–57. pmid:19709635

* View Article

* PubMed/NCBI

* Google Scholar

52. 52. Kanno K, Uchida A. Generalized synchronization and complexity in unidirectionally coupled dynamical systems. Nonlinear Theory Appl IEICE. 2012;3(2):143–54.

* View Article

* Google Scholar

53. 53. Nicola W, Clopath C. Supervised learning in spiking neural networks with FORCE training. Nat Commun. 2017;8(1):2208. pmid:29263361

* View Article

* PubMed/NCBI

* Google Scholar

54. 54. Takasu S, Aoyagi T. Suppression of chaos in a partially driven recurrent neural network. Phys Rev Res. 2024;6(1):013172.

* View Article

* Google Scholar

55. 55. Barančok P, Farkaš I. Memory capacity of input-driven echo state networks at the edge of chaos. In: Artificial Neural Networks and Machine Learning – ICANN 2014. 2014. p. 41–8.

56. 56. Ahmadian Y, Fumarola F, Miller KD. Properties of networks with partially structured and partially random connectivity. Phys Rev E Stat Nonlin Soft Matter Phys. 2015;91(1):012820. pmid:25679669

* View Article

* PubMed/NCBI

* Google Scholar

57. 57. Sommers H, Crisanti A, Sompolinsky H, Stein Y. Spectrum of large random asymmetric matrices. Phys Rev Lett. 1988;60(19):1895–8. pmid:10038170

* View Article

* PubMed/NCBI

* Google Scholar

58. 58. Rajan K, Abbott LF. Eigenvalue spectra of random matrices for neural networks. Phys Rev Lett. 2006;97(18):188104. pmid:17155583

* View Article

* PubMed/NCBI

* Google Scholar

59. 59. Tao T. Topics in random matrix theory. Providence, RI: American Mathematical Society; 2012.

60. 60. Baron JW, Galla T. Dispersal-induced instability in complex ecosystems. Nat Commun. 2020;11(1):6032. pmid:33247107

* View Article

* PubMed/NCBI

* Google Scholar

61. 61. Baron JW, Jewell TJ, Ryder C, Galla T. Eigenvalues of random matrices with generalized correlations: a path integral approach. Phys Rev Lett. 2022;128(12):120601. pmid:35394295

* View Article

* PubMed/NCBI

* Google Scholar

Citation: Tomita F, Teramae J-n (2025) Dynamical mean-field theory for a highly heterogeneous neural population with graded persistent activity of the entorhinal cortex. PLoS Comput Biol 21(9): e1013484. https://doi.org/10.1371/journal.pcbi.1013484

References

1. Squire LR, Wixted JT. The cognitive neuroscience of human memory since H.M.. Annu Rev Neurosci. 2011;34(1):259–88.

2. Garcia AD, Buffalo EA. Anatomy and function of the primate entorhinal cortex. Annu Rev Vis Sci. 2020;6:411–32. pmid:32580662

3. Persson BM, Ambrozova V, Duncan S, Wood ER, O’Connor AR, Ainge JA. Lateral entorhinal cortex lesions impair odor-context associative memory in male rats. J Neurosci Res. 2022;100(4):1030–46. pmid:35187710

4. MacDonald CJ, Lepage KQ, Eden UT, Eichenbaum H. Hippocampal “time cells” bridge the gap in memory for discontiguous events. Neuron. 2011;71(4):737–49. pmid:21867888

5. Eichenbaum H. Time cells in the hippocampus: a new dimension for mapping memories. Nat Rev Neurosci. 2014;15(11):732–44.

6. Sugar J, Moser M-B. Episodic memory: neuronal codes for what, where, and when. Hippocampus. 2019;29(12):1190–205. pmid:31334573

7. Umbach G, Kantak P, Jacobs J, Kahana M, Pfeiffer BE, Sperling M, et al. Time cells in the human hippocampus and entorhinal cortex support episodic memory. Proc Natl Acad Sci U S A. 2020;117(45):28463–74. pmid:33109718

8. Lin D, Huang AZ, Richards BA. Temporal encoding in deep reinforcement learning agents. Sci Rep. 2023;13(1):22335. pmid:38102369

9. Tsao A, Sugar J, Lu L, Wang C, Knierim JJ, Moser M-B, et al. Integrating time from experience in the lateral entorhinal cortex. Nature. 2018;561(7721):57–62. pmid:30158699

10. Bright IM, Meister MLR, Cruzado NA, Tiganj Z, Buffalo EA, Howard MW. A temporal record of the past with a spectrum of time constants in the monkey entorhinal cortex. Proc Natl Acad Sci U S A. 2020;117(33):20274–83. pmid:32747574

11. Egorov AV, Hamam BN, Fransén E, Hasselmo ME, Alonso AA. Graded persistent activity in entorhinal cortex neurons. Nature. 2002;420(6912):173–8. pmid:12432392

12. Frank LM, Brown EN. Persistent activity and memory in the entorhinal cortex. Trends Neurosci. 2003;26(8):400–1. pmid:12900167

13. Okamoto H, Tsuboshita Y, Fukai T. New horizons for associative memory broaden by neurophysiological and computational findings of graded persistent activity. Brain & Neural Networks. 2005;12(4):235–48.

14. Teramae J-N, Fukai T. A cellular mechanism for graded persistent activity in a model neuron and its implications in working memory. J Comput Neurosci. 2005;18(1):105–21. pmid:15789172

15. Hasselmo ME. Grid cell mechanisms and function: contributions of entorhinal persistent spiking and phase resetting. Hippocampus. 2008;18(12):1213–29. pmid:19021258

16. Yoshida M, Hasselmo ME. Persistent firing supported by an intrinsic cellular mechanism in a component of the head direction system. J Neurosci. 2009;29(15):4945–52. pmid:19369563

17. Koyluoglu OO, Pertzov Y, Manohar S, Husain M, Fiete IR. Fundamental bound on the persistence and capacity of short-term memory stored as graded persistent activity. Elife. 2017;6:e22225. pmid:28879851

18. Ebato Y, Nobukawa S, Sakemi Y, Nishimura H, Kanamaru T, Sviridova N, et al. Impact of time-history terms on reservoir dynamics and prediction accuracy in echo state networks. Sci Rep. 2024;14(1):8631. pmid:38622178

19. Stern M, Istrate N, Mazzucato L. A reservoir of timescales emerges in recurrent circuits with heterogeneous neural assemblies. Elife. 2023;12:e86552. pmid:38084779

20. Park JI, Lee DS, Lee SH, Park HJ. Incorporating heterogeneous interactions for ecological biodiversity. arXiv preprint 2024. https://arxiv.org/abs/physicsbio-ph

21. Sompolinsky H, Crisanti A, Sommers H. Chaos in random neural networks. Phys Rev Lett. 1988;61(3):259–62. pmid:10039285

22. Molgedey L, Schuchhardt J, Schuster H. Suppressing chaos in neural networks by noise. Phys Rev Lett. 1992;69(26):3717–9. pmid:10046895

23. Aljadeff J, Stern M, Sharpee T. Transition to chaos in random networks with cell-type-specific connectivity. Phys Rev Lett. 2015;114(8):088101. pmid:25768781

24. Harish O, Hansel D. Asynchronous rate chaos in spiking neuronal circuits. PLoS Comput Biol. 2015;11(7):e1004266.

25. Kadmon J, Sompolinsky H. Transition to Chaos in random neuronal networks. Physical Review X. 2015;5(4):041030.

26. Schuecker J, Goedeke S, Helias M. Optimal sequence memory in driven random networks. Physical Review X. 2018;8(4):041029.

27. Muscinelli SP, Gerstner W, Schwalger T. How single neuron properties shape chaotic dynamics and signal transmission in random neural networks. PLoS Comput Biol. 2019;15(6):e1007122. pmid:31181063

28. Helias M, Dahmen D. Statistical field theory for neural networks. 1st ed. Cham: Springer; 2020.

29. Keup C, Kühn T, Dahmen D, Helias M. Transient chaotic dimensionality expansion by recurrent networks. Physical Review X. 2021;11(2):021064.

30. Bordelon B, Pehlevan C. Self-consistent dynamical field theory of kernel evolution in wide neural networks. arXiv preprint 32240–56.

31. Krishnamurthy K, Can T, Schwab DJ. Theory of gating in recurrent neural networks. Phys Rev X. 2022;12(1):011011. pmid:36545030

32. Pereira-Obilinovic U, Aljadeff J, Brunel N. Forgetting leads to chaos in attractor networks. Phys Rev X. 2023;13(1):011009.

33. Engelken R, Wolf F, Abbott LF. Lyapunov spectra of chaotic recurrent neural networks. Phys Rev Res. 2023;5(4):043044.

34. Jaeger H. The “echo state” approach to analysing and training recurrent neural networks-with an erratum note. 148. Bonn, Germany: German National Research Center for Information Technology; 2001.

35. Maass W, Natschläger T, Markram H. Real-time computing without stable states: a new framework for neural computation based on perturbations. Neural Comput. 2002;14(11):2531–60. pmid:12433288

36. Lukoševičius M, Jaeger H. Reservoir computing approaches to recurrent neural network training. Computer Science Review. 2009;3(3):127–49.

37. Brody CD, Romo R, Kepecs A. Basic mechanisms for graded persistent activity: discrete attractors, continuous attractors, and dynamic representations. Curr Opin Neurobiol. 2003;13(2):204–11. pmid:12744975

38. Loewenstein Y, Sompolinsky H. Temporal integration by calcium dynamics in a model neuron. Nat Neurosci. 2003;6(9):961–7. pmid:12937421

39. Fransén E, Tahvildari B, Egorov AV, Hasselmo ME, Alonso AA. Mechanism of graded persistent cellular activity of entorhinal cortex layer v neurons. Neuron. 2006;49(5):735–46. pmid:16504948

40. Benettin G, Galgani L, Giorgilli A, Strelcyn JM. Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 2: Numerical application. Meccanica. 1980;15(1):21–30.

41. Shankar KH, Howard MW. A scale-invariant internal representation of time. Neural Comput. 2012;24(1):134–93.

42. Ganguli S, Huh D, Sompolinsky H. Memory traces in dynamical systems. Proc Natl Acad Sci U S A. 2008;105(48):18970–5. pmid:19020074

43. Goldman MS. Memory without feedback in a neural network. Neuron. 2009;61(4):621–34. pmid:19249281

44. Mazor O, Laurent G. Transient dynamics versus fixed points in odor representations by locust antennal lobe projection neurons. Neuron. 2005;48(4):661–73. pmid:16301181

45. Bondanelli G, Ostojic S. Coding with transient trajectories in recurrent neural networks. PLoS Comput Biol. 2020;16(2):e1007655. pmid:32053594

46. Pikovsky AS, Rosenblum MG, Osipov GV, Kurths J. Phase synchronization of chaotic oscillators by external driving. Physica D. 1997;104(3–4):219–38.

47. Rosenblum M, Pikovsky A, Kurths J. Phase synchronization in driven and coupled chaotic oscillators. IEEE Transactions on Circuits and Systems I-regular Papers. 1997;44(10):874–81.

48. Toral R, Mirasso CR, Hernandez-Garcia E, Piro O. Analytical and numerical studies of noise-induced synchronization of chaotic systems. Chaos. 2001;11(3):665–73. pmid:12779505

49. Teramae J-N, Tanaka D. Robustness of the noise-induced phase synchronization in a general class of limit cycle oscillators. Phys Rev Lett. 2004;93(20):204103. pmid:15600929

50. Nakao H, Arai K, Kawamura Y. Noise-induced synchronization and clustering in ensembles of uncoupled limit-cycle oscillators. Phys Rev Lett. 2007;98(18):184101. pmid:17501578

51. Sussillo D, Abbott LF. Generating coherent patterns of activity from chaotic neural networks. Neuron. 2009;63(4):544–57. pmid:19709635

52. Kanno K, Uchida A. Generalized synchronization and complexity in unidirectionally coupled dynamical systems. Nonlinear Theory Appl IEICE. 2012;3(2):143–54.

53. Nicola W, Clopath C. Supervised learning in spiking neural networks with FORCE training. Nat Commun. 2017;8(1):2208. pmid:29263361

54. Takasu S, Aoyagi T. Suppression of chaos in a partially driven recurrent neural network. Phys Rev Res. 2024;6(1):013172.

55. Barančok P, Farkaš I. Memory capacity of input-driven echo state networks at the edge of chaos. In: Artificial Neural Networks and Machine Learning – ICANN 2014. 2014. p. 41–8.

56. Ahmadian Y, Fumarola F, Miller KD. Properties of networks with partially structured and partially random connectivity. Phys Rev E Stat Nonlin Soft Matter Phys. 2015;91(1):012820. pmid:25679669

57. Sommers H, Crisanti A, Sompolinsky H, Stein Y. Spectrum of large random asymmetric matrices. Phys Rev Lett. 1988;60(19):1895–8. pmid:10038170

58. Rajan K, Abbott LF. Eigenvalue spectra of random matrices for neural networks. Phys Rev Lett. 2006;97(18):188104. pmid:17155583

59. Tao T. Topics in random matrix theory. Providence, RI: American Mathematical Society; 2012.

60. Baron JW, Galla T. Dispersal-induced instability in complex ecosystems. Nat Commun. 2020;11(1):6032. pmid:33247107

61. Baron JW, Jewell TJ, Ryder C, Galla T. Eigenvalues of random matrices with generalized correlations: a path integral approach. Phys Rev Lett. 2022;128(12):120601. pmid:35394295