Phase relations of interneuronal activity relative to theta rhythm

The theta rhythm plays a crucial role in synchronizing neural activity during attention and memory processes. However, the mechanisms behind the formation of neural activity during theta rhythm generation remain unknown. To address this, we propose a mathematical model that explains the distribution of interneurons in the CA1 field during the theta rhythm phase. Our model consists of a network of seven types of interneurons in the CA1 field that receive inputs from the CA3 field, entorhinal cortex, and local pyramidal neurons in the CA1 field. By adjusting the parameters of the connections in the model. We demonstrate that it is possible to replicate the experimentally observed phase relations between interneurons and the theta rhythm. Our model predicts that populations of interneurons receive unimodal excitation and inhibition with coinciding peaks, and that excitation dominates to determine the firing dynamics of interneurons.

1. Introduction

The hippocampus is a brain structure that plays a key role in the processes of attention and memory. To process information, neural ensembles in the hippocampus need to be synchronized with rhythms. The main rhythm that organizes the neural activity of the hippocampus during cognitive tasks is the theta rhythm (4–12 Hz) (Vinogradova, 1995; Buzsáki, 2002; Buzsáki and Moser, 2013). Almost all hippocampal neurons are modulated by theta rhythm (Mizuseki et al., 2009). This is expressed in the fact that each population has a theta rhythm phase, in which the probability of discharges of its neurons is maximal. Modulation of neuronal activity by rhythm makes it possible to synchronize different areas of the hippocampal formation during information processing (Fries, 2015; Nuñez and Buño, 2021; Mysin and Shubina, 2022). In addition, the theta rhythm provides an ordered structure of the place cell activity due to phase precession (Burgess and O'Keefe, 2011; Jaramillo and Kempter, 2017). Thus, understanding the mechanisms of theta rhythm formation is the most important problem of neuroscience.

Our study has two aims. The first aim is to explain the effect of modulation of the firing rate of the interneurons of the CA1 field by the theta rhythm. The second aim is to adapt the optimization methods of spike network models to describe experimental data.

The mechanisms of formation of phase relations of interneuronal activity and theta rhythm are unknown. Several theoretical studies have investigated the formation of phase relations between neurons of the CA1 field and theta rhythm (Bezaire et al., 2016; Mysin, 2021). Despite the detailed nature of these models, the authors were unable to reproduce the form of distribution of most types of interneurons in the theta rhythm phase. The main problem in constructing a model of the distribution of neurons by theta rhythm phases is to identify the mechanism of stabilization of antagonistic relationships between different populations of interneurons. For example, parvalbumin-containing (PV) and cholecystokinin-containing (CCK) basket neurons inhibit each other and fire at opposite phases of the theta rhythm (Lasztóczi et al., 2011; Royer et al., 2012; Varga et al., 2012; Bezaire and Soltesz, 2013; Somogyi et al., 2014). It is difficult to choose the parameters of connections when these neurons form a stable antiphase activity. Most often, one of the populations completely inhibit the other. When considering a larger number of interneuron populations, the problem becomes much more complicated (Bezaire et al., 2016; Mysin, 2021).

We hypothesized that short-term synaptic plasticity stabilizes the structure of interneuronal activity. In a recent study, the authors of the project Hippocampome.org provide a large meta-analysis of data on short-term synaptic plasticity in the hippocampus (Moradi et al., 2022). We used the estimates of this study as the basis of our model. We have considered the mechanism of formation of phase relations for populations of interneurons of the CA1 field. The CA1 field was chosen because there is the greatest amount of data on interneuron activity for this region of the hippocampus. We modeled seven populations of interneurons: PV (pvbas) and CCK (cckbas) basket, axo-axonal (aac), OLM (oriens-lacunosum moleculare) (olm), Ivy (ivy), and neurogliaform (ngf) and bistratified (bis) neurons. These populations were selected because the phase relationships for them are described (Somogyi and Klausberger, 2005; Fuentealba et al., 2010; Lapray et al., 2012; Varga et al., 2012). Neurons in these populations make up ~70% of all interneurons in the CA1 field (Bezaire and Soltesz, 2013). We also took into account inputs from pyramidal neurons of the CA3 and CA1 fields, as well as from neurons of the third layer of the entorhinal cortex. We were able to select the parameters of neurons and connections in the model and found a good description of the mechanisms of phase relations.

Considering a network of seven populations requires setting up several hundred parameters. This makes it impossible to apply gradient-free optimization methods. We were able to adapt the model to apply gradient descent optimization. This approach can be applied to modeling problems of other experimental phenomena in neural networks of the brain.

2. Materials and methods

2.1. Models of neurons

We consider a network of seven populations (Figure 1). The neurons of each population are described by the leaky integrate-and-fire (LIF) model. We used simulations using the conduction-based refractory density (CBRD) approach (Section 2.3) and direct stimulation of neurons with noise input (Monte Carlo). In this section, we will only describe the equations used. In Section 3.1, we will discuss in detail the motivation for using each approach for modeling. Monte Carlo simulations were provided with equations as follows:

$\begin{array}{l}{C}_m\frac{dV}{dt}=g_L(E_L-V)+I_{ext}+\sum I_{syn}+{\sigma }_m·\eta (t)\\ \text{if}V>V^T:V=V_{reset}& (1)\end{array}$

Here, C_m is the membrane capacitance, V is the membrane potential, g_L is the leak conductance, E_L is the leak reverse potential, I_ext is the external current, and I_syn is the synaptic current. σ_m is the standard deviation of noise $\eta ~\mathcal{N}(0,1)$. For all neurons in all simulations, we used following values: $C_m=1\mu F/c{m}^2,g_L=0.1mS/c{m}^2,E_L=-60mV,V^T=-50mV,V_{reset}=-90mV,{\sigma }_m=0.3\mu A/c{m}^2/\sqrt{ms}$. Value I_ext is optimized for each population. After generating the action potential, the refractory = 3 ms. We used typical values used in simulations. The exact values of the parameters can be viewed in the Hippocampome database (Venkadesh et al., 2019)

FIGURE 1

Figure 1. Model architecture. Hippocampal CA1 microcircuit is simulated. Cell types and their connectivity are shown. ca3pyr, pyramidal neurons of the CA3 field; ca1pyr, pyramidal neurons of the CA1 field; ec3, principal neurons of the 3 layer of the entorhinal cortex; pvbas, PV basket cells; olm, Oriens-lacunosum molecular cells; cckbas, CCK basket cells; ivy, Ivy cells; ngf, neurogliaform cells; bis, bistratified cells; aac, axo-axonal cells. Neurons of the third layer of the entorhinal cortex and pyramidal neurons of the CA3 and CA1 fields are not simulated explicitly. They were taken into account as external inputs to the model. All interneurons are simulated with the LIF model.

Population frequency for Monte Carlo simulations as follows:

$\begin{array}{ll}\nu =\frac{n_{fired}}{N·\text{Δ}t}& (2)\end{array}$

where n_fired is the number of fired cell at each time step, N = 4, 000 is the full number of neurons in population, and Δt = 0.1 ms is the integration step. We used the stochastic Heun method from the brian2 package for Monte Carlo simulations (Stimberg et al., 2019).

2.2. Synapse model

Synapses were simulated with the Tsodyks—Markram model (Tsodyks and Markram, 1997; Tsodyks et al., 1998; Moradi et al., 2022) as follows:

$\begin{array}{ll}\frac{du}{dt}=-\frac{u}{{\tau }_f}+U_{inc}·(1-u_{-})·w·{\nu }_{pre}·dt& (3)\end{array}$

$\begin{array}{ll}\frac{dx}{dt}=\frac{1-x-y}{{\tau }_r}-u_{+}{x}_{-}·w·{\nu }_{pre}·dt& (4)\end{array}$

$\begin{array}{ll}\frac{dy}{dt}=-\frac{y}{{\tau }_d}+u_{+}{x}_{-}·w·{\nu }_{pre}·dt& (5)\end{array}$

U_inc is the increment of u produced by a spike. u₋, x₋ are variables just before the arrival of the spike, and u₊ refers to the moment just after the spike. From the first equation (Equation 3), u₊ = u₋ + U_inc·(1 − u₋). Variables y and x mean the activated, deactivated, and recovered states, respectively, u is a fraction of resource ready for use. After each presynaptic spike, an instantaneous shift occurs from recovered to activated state. The amount of shift is determined by u. The active resources, then, decay to the deactivated state by the decay time constant τ_d. Since synaptic resources are limited, the more resources stay in the deactivated state, the more a synapse is depressed. Synaptic resources exponentially recover from depression with the recovery time constant τ_r. w is normalization coefficient, it makes sense of the density of connections. ν_pre is firing rate of presynaptic population. In Monte Carlo simulations ν_pre is defined by Equation 2, and in CBRD simulations, ν_pre is defined by Equation 9. For each connection parameters U_inc, τ_d, τ_r, τ_f, w, and g_{syn, max} are optimized. We consider infinitely large populations of presynaptic and postsynaptic neurons. In this context, variables mean the averages of all synapses between two populations [for detail, see Section 3 in Tsodyks et al. (1998)].

The synaptic current as follows:

$\begin{array}{ll}{I}_{syn}=g_{syn,max}·y·(E_{syn}-V)& (6)\end{array}$

where g_{syn, max} is maximal synaptic conductance and E_syn is reversal potential for synaptic current. The connection probability in Monte Carlo simulations is 0.5.

2.3. Conductance-based refractory density approach

Conductance-based refractory density (CBRD) approach simulates the distribution (ρ) of neurons in space-times after spike generation (t^*) (Chizhov and Graham, 2007, 2008; Chizhov, 2017). The function ρ(t, t^*) is described by the transfer equation as follows:

$\begin{array}{ll}\frac{\partial \rho }{\partial t}+\frac{\partial \rho }{\partial t^{*}}=\rho H& (7)\end{array}$

$\begin{array}{ll}{\displaystyle {\int }_0^{\infty }\rho ·d{t}^{*}}=1& (8)\end{array}$

The function H is defined as the probability of a single neuron to generate a spike when the mathematical expectation of the neuronal state variables is known (Equation 11).

Population firing rate is defined as follows:

$\begin{array}{ll}\nu (t)\equiv \rho (t,0)={\displaystyle {\int }_0^{\infty }\rho ·H·d{t}^{*}}& (9)\end{array}$

Population firing rate is determined by the probability of finding a neuron in the state of spike generation (t^* = 0).

The CBRD approach consider neurons only according to t^* variable. The state of each neuron is parameterized by this phase variable. For LIF neuron, its only state variable is V. Neuron model modified as follows:

$\begin{array}{ll}{C}_m(\frac{\partial V}{\partial t}+\frac{\partial V}{\partial t^{*}})=g_L(E_L-V)+I_{ext}+{\displaystyle \sum I_{syn}}& (10)\end{array}$

The function H has two components A and B:

$\begin{array}{ll}H=\frac{A+B}{{\tau }_M}& (11)\end{array}$

A is the probability for a neuron to cross the threshold because of noise, derived analytically in the study mentioned in Chizhov and Graham (2007) and approximated by exponential and polynomial functions for convenience. B is the probability of dischange of a neuron during depolarization due to external input, i.e., the hazard from drift in the voltage phase space. τ_M is time of membrane:

$\begin{array}{ll}{\tau }_M=\frac{g_{tot}}{C_m}& (12)\end{array}$

C_m is capacity of membrane, g_tot is total conductance of all channels, and it is sum of g_L and all synaptic conductances.

$\begin{array}{ll}A=exp(0.0061-1.12T-0.257T^2-0.072T^3-0.0117T^4)& (13)\end{array}$

$\begin{array}{ll}B=-\sqrt{2}·{\left[\frac{dT}{dt}\right]}_{+}·F_T·{\tau }_M& (14)\end{array}$

[z]₊ is heaviside function:

$\begin{array}{ll}{[z]}_{+}:=\{\begin{array}{ll}0\hfill & \text{if }z\le 0;\hfill \\ 1\hfill & \text{if }z>1\hfill \end{array}& (15)\end{array}$

$\begin{array}{ll}T=\frac{V^T-V}{\sqrt{2}{\sigma }_m}& (16)\end{array}$

where V is membrane potential, σ_m is the noise amplitude, and V^T is the threshold for spike generation. T is the membrane potential relative to the threshold, scaled by noise amplitude.

$\begin{array}{ll}\frac{dT}{dt}=\frac{-1}{\sqrt{2}{\sigma }_m}\frac{dV}{dt}& (17)\end{array}$

$\begin{array}{ll}{F}_T=\sqrt{2/\pi }·exp(-(T^2))/(1+erf(T))& (18)\end{array}$

where erf is Gauss error function. The function F_T characterizes the firing probability in the regime of fast changes of T.

All values for model are same as for Monte Carlo simulations.

2.4. Numerical methods for CBRD simulations

In general, the equations can be written as follows:

$\begin{array}{ll}\frac{\partial z}{\partial t}+\frac{\partial z}{\partial t^{*}}=S_z(t,t^{*})& (19)\end{array}$

where z(t, t^*) is one of the functions ρ or V, and $S_z(t,t^{*})$ is the source term of the equation.

Numerical scheme was taken from Harten and Osher (1987). Numerical scheme is as follows:

$\begin{array}{ll}{z}_i^{n+1}=z_i^n-\frac{\text{Δ}t}{\text{Δ}{t}^{*}}\left(z_i^n-z_{i-1}^n+\frac{1}{2}\left(1-\frac{\text{Δ}t}{\text{Δ}{t}^{*}}\right)(w_i^z-w_{i-1}^z)\right)+\text{Δ}t·{(S_z)}_i^n& (20)\end{array}$

n is index in time, and i is index in space.

$\begin{array}{ll}{w}_i^z=limiter(z_{i+1}-z_i,z_i-z_{i-1})& (21)\end{array}$

$\begin{array}{ll}limiter(a,b)=\{\begin{array}{ll}0\hfill & \text{if }ab<=0;\hfill \\ -min(0.5\text{ |}ab\text{|},2min(\text{|}a\text{|},\text{|}b\text{|}))\hfill & \text{if }a<0,ab>0;\hfill \\ min(0.5\text{ |}ab\text{|},2min(\text{|}a\text{|},\text{|}b\text{|}))\hfill & \text{otherwise}\hfill \end{array}& (22)\end{array}$

2.4.1. Bound conditions

Bound values $w_0^z=0$, $w_N^z=0$. Firing rate (left bound for ρ) is calculated as follows:

$\begin{array}{ll}{\rho }_0^{n+1}={\rho }_0^n-\frac{\text{Δ}t}{\text{Δ}{t}^{*}}{\rho }_0^n-\text{Δ}t{\displaystyle \sum _{i=0}^N}{(S_{\rho })}_i^n& (23)\end{array}$

Left bound for V:

$\begin{array}{ll}{V}_0^{n+1}=V_{reset}& (24)\end{array}$

Right bound for ρ:

$\begin{array}{ll}{\rho }_N^{n+1}={\rho }_N^n+\frac{\text{Δ}t}{\text{Δ}{t}^{*}}({\rho }_{N-1}^n+\frac{1}{2}(1-\frac{\text{Δ}t}{\text{Δ}{t}^{*}}\left)w_{N-1}^{\rho }\right)+\text{Δ}t·{(S_{\rho })}_N^n& (25)\end{array}$

Right bound for V:

$\begin{array}{ll}{V}_N^{n+1}=V_N^n+\text{Δ}t{(S_V)}_N^n& (26)\end{array}$

Δt = 0.1 ms, Δt^* = 0.5 ms, and N = 400 is the number of spatial states. All simulations using the CBRD approach were performed with the TensorFlow package (version 2.10.0). All gradients were calculated using automatic differentiation (Abadi et al., 2015).

2.5. Target firing rates, inputs, and loss function

Firing rate of inputs and target firing rate for simulated populations as follows:

$\begin{array}{ll}F{R}_{pop}(t)=\frac{F{R}_{mean,pop}}{I_0({\kappa }_{pop})}·exp({\kappa }_{pop}·(cos(2\pi {\omega }_{\theta }t-{\varphi }_{pop})))& (27)\end{array}$

where FR_pop(t) is population firing rate in time in spikes/second. FR_mean,pop is mean population firing rate (spikes/second). κ_pop is a measure of concentration of von Mises distribution, I₀(κ) is a zero order Bessel function. ω_θ = 7 Hz is frequency of theta rhythm, and ϕ_pop is peak phase of population firing. Parameters for all populations are presented in Table 1. Authors of experimental studies usually use R (ray length) as a measure of phase variation rather than κ (Mardia and Jupp, 1999). κ is calculated from R with following approximation (Mardia and Jupp, 1999):

$\begin{array}{ll}\kappa (R)=\{\begin{array}{ll}2·R+R^3+5/6·R^5\hfill & \text{if }R<0.53;\hfill \\ -0.4+1.39·R+0.43/(1-R)\hfill & \text{if }0.53\le R<0.85;\hfill \\ 1/(3·R-4·R^2+R^3)\hfill & \text{if }R\ge 0.85\hfill \end{array}& (28)\end{array}$

TABLE 1

Table 1. Parameters of the target function.

Artificial inputs in Monte Carlo simulations are modeled as Poisson generators with rate given by the Equation 27.

The loss function for estimating the discrepancy between the simulation and the target function as follows:

$\begin{array}{ll}{L}_{simulation}={\displaystyle \sum _{k=1}^K}{\displaystyle \sum _{t_n=0}^T}{(ln({\rho }_{target}(t_n)+1)-ln({\rho }_{simulated}(t_n)+1))}^2& (29)\end{array}$

Summation is carried out by time and all populations of neurons. K = 7 is the number of populations in the model. T is the number of time steps in the simulation.

The optimized parameters must be within certain bounds. In particular, all parameters of synapses should be positive, and U_inc should not exceed one. We have added barrier terms to the loss function.

$\begin{array}{l}{L}_{barrier}={\displaystyle \sum _{m=1}^M(\text{−}0.001·ln(100g_{syn,max}))}+\\ +{\displaystyle \sum _{m=1}^M(\text{−}0.001·ln(100{\tau }_r))}+\\ +{\displaystyle \sum _{m=1}^M(\text{−}0.001·ln(100{\tau }_f))}+\\ +{\displaystyle \sum _{m=1}^M(\text{−}0.001·ln(100{\tau }_d))}+\\ +{\displaystyle \sum _{m=1}^M(\text{−}0.001·ln(100U_{inc}))}+\\ +{\displaystyle \sum _{m=1}^M(\text{−}0.001·ln(100w))}+\\ +{\displaystyle \sum _{m=1}^M(\text{−}0.001·ln(100(1-U_{inc}))}& (30)\end{array}$

Summarization is carried out for all synapses in the model. M = 49 M is the number of connections in the model.

Full loss function for optimization as follows:

$\begin{array}{ll}{L}_{full}=L_{simulation}+L_{barrier}& (31)\end{array}$

We have used Adam optimizer with standard parameters: learning rate = 0.001, β₁ = 0.9, and β₂ = 0.999 (Kingma and Ba, 2014).

3. Results

3.1. Optimization of the model

We considered a network of biologically plausible networks consisting of seven populations of interneurons. The graph of connections is complex and contains 49 connections (Figure 1) (Moradi et al., 2022). Excitatory inputs were not modeled as neural populations but were taken into account artificial generators (functions from time). Each synapse is described by six parameters (Equations 3–5). We fit all these parameters. In addition, the I_ext parameter for each neuron population was fit in simulations (Equation 1). Thus, we needed to tune 301 model parameters. Our approach contains elements of novelty, so we will begin the description of the results with a brief description of the formulation and solution of the optimization problem.

We propose to fit the parameters by reducing the process of building a model to solving the optimization problem. In other words, it is proposed to find the optimal parameters of the model so that it describes the experimental data in the best way. The requirement to be able to configure numerous parameters makes it impossible to use zero-order optimization methods. Gradient estimation is possible by using automatic differentiation packages (Tensorflow, Torch) (Abadi et al., 2015) or the adjoint state method (Chen et al., 2019; Sun et al., 2020). However, these approaches require direct computation of variables by which differentiation will be carried out.

The activity of neurons in the real brain is noisy. The vast majority of experimental phenomena are obtained as a result of averaging the activity of one or more neurons depending on external stimuli or internal state. Modulation of neuronal activity by theta rhythm is just one example. In other words, the experimental phenomena of neural network dynamics are empirical distributions. Modeling consists of describing the parameters of distributions. In this context, direct simulation of distributions is necessary to be able to estimate gradients.

We have used the population approach for modeling, since it allows us to describe the population firing rate. As a population approach, we chose the CBRD (conductance-based refectory density) method (Chizhov and Graham, 2007, 2008; Chizhov, 2017). Direct simulation of the firing rate makes it possible to estimate the gradient of the loss function from the model parameters. All operators of the numerical scheme of the CBRD model are differentiable, so automatic differentiation can be applied (Section 2.4). Estimating the gradient of the loss function makes it possible to apply gradient descent methods to optimize model parameters (Equation 31). All simulations using the CBRD approach were performed with the TensorFlow package, and gradients were calculated using automatic differentiation (Abadi et al., 2015). As initial conditions for optimization, we used the average values of the parameters described in the study by Moradi et al. (2022). We ran Monte Carlo simulations with the parameters found using optimization. Monte Carlo simulations were used to control of the accuracy of the CBRD approach.

We have introduced a target activity function for each population (Equation 27). This function approximate experimental data about modulation of neuronal activity by the theta rhythm. It describes the population firing rate of each group of neurons over time. The target function takes into account the average firing rate and theta rhythm modulation. Theta rhythm modulation is determined by the mean phase (circular average of the phases of neuronal discharges) and the phase variation (ray length - R). The population firing rate of excitatory inputs was described by similarly functions (Equation 27).

This approach has advantages and disadvantages. On the one hand, population models are only applicable to simple models of neurons. In our study, we used the LIF (leaky integrate-and-fire) model (Equation 1). On the other hand, the use of gradient descent allows optimization over a much larger number of parameters than gradient-free methods.

3.2. Description of the model results

Figure 2 shows the optimization results. For each population, the target function, the results of simulation using the CBRD approach, and the results of Monte Carlo simulation with optimal parameters are presented. The model under consideration can approximate the target functions for all populations of interneurons with significant accuracy. The simulation results are in the confidence intervals of experimental estimates. The model describes all the required characteristics of the experimental data: the average spike rate, the average phase of the theta rhythm, and the ray length (R). We also note a good agreement between the results of Monte Carlo simulations and simulations using the CBRD approach. In the next step, we looked at how the model works. Figure 3 shows the dynamics of synaptic conductivities in each neuron population. Both excitatory and inhibitory conductivities have an unimodal distribution over the theta cycle for all neuron populations, except neurogliaform cells. An unexpected result is that the peaks of the excitatory and inhibitory inputs coincide. This effect is observed for all populations of neurons. The ratio of the total excitatory to inhibitory conductivity remains approximately constant over the theta cycle and ranges from 0.7 to 2.4 for different populations. Since the reversion potential for AMPA receptors (0 mV) is much further from the resting potential than the reversion potential of GABA-A receptors (–75 mV). Excitatory currents dominate over inhibitory ones and determine the dynamics of neuronal discharges. For all neuron populations except axo-axonal ones, the peak of excitatory input is located near the peak of firing rate. The dynamics of axo-axonal neurons in the model is determined by the external current (Supplementary Table S1). All optimal parameters are given in Supplementary material.

FIGURE 2

Figure 2. The results of model optimization. For each population, plots show the target function, and the population spike rate obtained with the CBRD approach using Monte Carlo simulation. In total, 1 second of simulation is shown after stabilization of the model dynamic mode. The notation of neurons is similar to Figure 1.

FIGURE 3

Figure 3. Synaptic conductivities. Two theta cycles of the simulation using the CBRD approach are shown. Excitatory (upper series of plots) and inhibitory (lower series of plots) conductivities are shown for each population. The color indicates presynaptic populations. Each plot shows the sum of conductivities, “sum exc" and “sum inh" notes sum of excitatory and inhibitory conductivities, respectively. The notation of neurons is similar to Figure 1.

We performed optimization for a similar network with synapses without short-term plasticity. The results are presented in the Supplementary material. For a part of neural populations (pvbas, olm, bis, and aac), it fits the parameters to describe phase relations. However, the remaining populations are completely silent (ivy, ngf, and cckbas). These results show that short-term plasticity is a necessary element of the model for reproducing phase relations.

3.3. Generalizing power of the model

At the last stage of research, we tested the stability of the found solution at other theta rhythm frequencies. We searched for optimal parameters at a theta rhythm frequency of 7 Hz. Figure 4 shows the results of simulations with varying theta rhythm frequency. Figure 5 shows raster plots for a Monte Carlo simulation with a theta rhythm frequency of 5 Hz. In this series of computational experiments, we changed the frequencies of the input populations. The model parameters were used that were found earlier. Over the entire frequency range (4–12 Hz), the shape of the distribution of neuronal activity over the theta rhythm phase is preserved. The theta rhythm frequency is a parameter that varies greatly in animal experiments, for example, it correlates with the animal's running speed (Hinman et al., 2011; Ledberg and Robbe, 2011; Long et al., 2014; Justus et al., 2017). The conservation of the phase relations of interneuron activity with a change in the theta rhythm frequency is an important property of our model. These experiments can be considered as analogous to cross-validation used in machine learning. In the language of machine learning, we can say that our model shows good quality on the “test set”, i.e., samples that did not participate in the “training”.

FIGURE 4

Figure 4. Phase relations. Each plot shows the distribution of activity of neuronal populations by theta rhythm phase. The frequency of the theta rhythm was set in the equations of the exciting inputs. All simulations were carried out with optimal parameters found at a theta rhythm frequency of 7 Hz. Simulations were carried out using the CBRD approach (thick lines) and Monte Carlo (dashed line). The notation of neurons is similar to Figure 1.

FIGURE 5

Figure 5. Raster plots for Monte Carlo simulation with the theta rhythm frequency of 5 Hz. For each population, only 1,000 neurons are shown, but 4,000 neurons are simulated. The dashed line is reference cos. The notation of neurons is similar to Figure 1.

4. Discussion

4.1. Assumptions and limitations of the model

Our results are based on several key assumptions. We used the von Mises function to describe the inputs and target activities for interneurons (Equation 27). These assumptions of the model are based on experimental data that more than 85% of pyramidal neurons in the CA3 and CA1, three fields of the EC layer are modulated by theta rhythm (Mizuseki et al., 2009). The proportion of interneurons of the CA1 field modulated by the theta rhythm is 95% (Mizuseki et al., 2009). The degree of phase modulation R for all hippocampal neurons is in the range of 0.2–0.4; we used R = 0.3 as the mean estimate. Binding phases for exciting inputs and local interneurons of the CA1 field (Table 1) have been measured in several experimental studies and are reliable (Klausberger et al., 2003; Somogyi and Klausberger, 2005; Mizuseki et al., 2009; Somogyi et al., 2014). Our model takes into account the average discharge frequencies of input populations and interneurons (Table 1). We selected the values for simulations as the mean and median values measured in the literature. The exact values of the mean firing rate are not important in our context, since the discharge frequency of the presynaptic population is linearly related to the weight of the synaptic connection with the postsynaptic population. In other words, a decrease in the firing rate of the presynaptic population can be compensated by an increase in synaptic conduction in postsynaptic neurons. For example, if a neuron receives inputs from N neurons with a connection density p, there are N·p synapses on it. If N is decreased but p is increased by the same number of times, the number of synapses will not change. There is a coefficient w in the synapse model (Equations 3–5). It takes into account the density of connections. We have optimized the coefficient w, so the density of connections is indirectly taken into account in the model. Thus, our model results describe the network of the CA1 field up to normalization.

We used a homogeneous representation of each population of interneurons. Recent studies show that interneurons are involved in the formation of neural ensembles, i.e., their activity is modulated by external stimuli (Geiller et al., 2022). However, a study of simultaneous registration of a large number of interneurons of different classes shows that cells of all classes are active simultaneously for times of hundreds of milliseconds (Geiller et al., 2020). Thus, it is reasonable to assume that a significant proportion of neurons in each population is active at the same time. Notably, the effect of simultaneous activity of interneurons of different populations is observed, despite the inhibitory connections of populations between each other (Bezaire and Soltesz, 2013). Our model reproduces this effect well. Due to the homogeneity of each population, we used the same number of neurons in direct simulations. Interneurons are unevenly distributed across classes. For example, population of Ivy neurons are six times larger than axo-axonal cells (Bezaire and Soltesz, 2013). A representative representation of neurons would not improve the accuracy of the model, but it would increase computational costs. We did not take into account the input to the model from the medial septum. Recent studies show that the strongest rhythmic input from the septum goes to the interneurons of the CA3 field and the entorhinal cortex (Joshi et al., 2017; Unal et al., 2018; Viney et al., 2018). Other data demonstrate the importance of input from the CA3 field and the entorhinal cortex to the CA1 field for theta rhythm in it (López-Madrona et al., 2020; Zutshi et al., 2022). Together, these results lead to the hypothesis that the theta rhythm in the CA1 field results from secondary rhythmic inputs rather than direct input from the medial septum. There are two significant limitations of our model. The first is usage of the LIF model of neurons. Most types of interneurons have slow potassium channels (Bezaire et al., 2016; Komendantov et al., 2019), which can make a significant contribution to the frequency of neuronal discharges due to spike-frequency adaptation. OLM neurons have pronounced H-currents (Saraga et al., 2003), which can also make a significant contribution to the impulsing of these neurons due to depolarization or resonant properties (Avella Gonzalez et al., 2015). The usage of the LIF model is dictated by a limitation of the population approach. We did not explicitly simulate pyramidal neurons of the CA1 field due to the limitations of the population method. Different classes of interneurons have entrances to different compartments of pyramidal cells. The correct inclusion of pyramid neurons in the network requires the use of multicompartment models to describe them. The problem of finding optimal synaptic inputs to pyramidal neurons can be solved separately.

The second limitation is that we did not take into account all populations of interneurons. This limitation is based on the weak knowledge of other populations of interneurons. The lack of data makes it impossible to include other populations in the model. Accounting for more populations can change the balance of inhibition and shift its peak.

4.2. Comparison with other models

The literature presents several attempts to explain the phase relations between neurons and theta rhythm. The key idea of the research is to determine the optimal structure of the input from the medial septum and excitatory inputs (Cutsuridis et al., 2010; Cutsuridis and Hasselmo, 2012; Cutsuridis and Poirazi, 2015; Bezaire et al., 2016; Mysin et al., 2019; Mysin, 2021). Our results do not negate the contribution of external inputs but rather dismantle an additional mechanism for stabilizing phase relations. Models with a few interneuron populations give a good approximation of the distribution of neurons in theta rhythm phases due to input from the medial septum (Cutsuridis et al., 2010; Mysin et al., 2019). However, with an increase in the number of interneuron classes in models, the quality of the approximation of phase relations decreases (Cutsuridis and Poirazi, 2015; Bezaire et al., 2016; Mysin, 2021). Although in the cited studies, the authors used a different input structure. We believe that taking into account the detailed structure of the inputs to the CA1 field in the model is not sufficient to explain the phase relations of interneurons. In this study, we have shown that the short-term plasticity of synapses between populations of interneurons may be the missing component in the stabilization of phase relationships.

4.3. Comparison of simulation results with experimental data

The main result of our study is an unimodal distribution of excitation and inhibition with coinciding peaks for most types of interneurons. Checking this fact is a complex experimental task. In the literature, there is a study of synaptic currents on PV basket and OLM neurons during theta rhythm generation in hippocampal slices (Huh et al., 2016). The authors also found that unimodal excitation and inhibition coincided in time. The peak of inhibition occurred 12 ms after the peak of excitation for PV basket neurons and 7 ms after the peak of arousal for OLM neurons (Huh et al., 2016). We used experimental data in vivo as the basis for our model; therefore, comparison with data on hippocampal preparations is limited.

4.4. Prospects for applying population model optimization

Several attempts have been made in the literature to optimize networks of spiking neurons (Lee et al., 2016; Wunderlich and Pehle, 2021; Zenke and Vogels, 2021; Ramezanian-Panahi et al., 2022). However, these studies are aimed at finding practical applications in the field of data analysis not brain modeling. These approaches are not directly applicable for building a model in neuroscience. The activity of real individual neurons is noisy. The vast majority of experimental phenomena are obtained as a result of averaging the activity of one or several neurons, depending on external stimuli or internal state. The phase relations of neuronal activity are one of the many examples. Another example is place cells. This phenomenon can be described as the distribution of neuronal activity depending on the position of the animal. Phase precession of place cells is the distribution of neuronal activity depending on the position of the animal and the phase of the theta rhythm. Population models describe physiological phenomena in the language of distributions. Modeling complex systems require numerous equations and parameters. Projects such as the Human Brain Project or the Hippocampus Project aim to model the brain by collecting and organizing parameters for equations (Wheeler et al., 2015; Bjerke et al., 2018). Optimization of population models can well complement the processes of collecting parameters. The development of the mathematical apparatus of population models may be one of the reasons for a breakthrough in the construction of large-scale models of the brain.

4.5. General remarks

Our model shows a possible mechanism for the formation of phase relations between interneurons and theta rhythm in the CA1 field. We assume that a similar mechanism operates in other areas of the hippocampal formation. The establishment of phase relationships among interneurons leads to the coupling of principal neurons to the phase of the theta rhythm. Rhythmic inputs to the principal neurons creates oscillations of local field potential (Buzsáki et al., 2012; Einevoll et al., 2013). This, in turn, synchronizes the regions of the hippocampus with each other through the projections of principal neurons, which is necessary for the transmission of information (Bastos et al., 2015; Mysin and Shubina, 2022). Synchronization of hippocampal regions creates a capability for the formation of phase precession (Burgess and O'Keefe, 2011; Fernández-Ruiz et al., 2017; Zutshi et al., 2022). Another result of the synchronization of different areas of the hippocampus is the emergence of coupling of theta and gamma rhythms. For example, a slow gamma rhythm is formed in a feedback loop between PV basket and principal neurons (Colgin and Moser, 2010; Buzsáki and Wang, 2012). In the CA1 field, excitation from the CA3 field is necessary for the formation of a slow gamma rhythm (Colgin, 2015; Fernández-Ruiz et al., 2017). PV basket neurons CA1 have a peak of discharges in the descending phase of theta rhythm, and the input from field CA3 falls into the same phase (Mizuseki et al., 2009; Belluscio et al., 2012). Therefore, a model that reproduces phase relationships with respect to the theta rhythm will have a good predictive ability with respect to other phenomena associated with the theta rhythm. Identification of the mechanisms of formation of the dynamics of the activity of interneurons during the generation of the theta rhythm is important for understanding the processing of information and the formation of hippocampal rhythms.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found at: https://static-content.springer.com/esm/art%3A10.1038%2Fs42003-022-03329-5/MediaObjects/42003_2022_3329_MOESM3_ESM.xlsx.

Author contributions

IM: main idea, coding, and article writing.

Funding

This study was supported by the Russian Science Foundation, Grant No. 20-71-10109.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncir.2023.1198573/full#supplementary-material

References

Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., et al (2015). TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems. Software Available online at: tensorflow.org (accessed May 6, 2023).

Google Scholar

Avella Gonzalez, O. J., Mansvelder, H. D., van Pelt, J., and van Ooyen, A. (2015). H-channels affect frequency, power and amplitude fluctuations of neuronal network oscillations. Front. Comput. Neurosci. 9, 141. doi: 10.3389/fncom.2015.00141

PubMed Abstract | CrossRef Full Text | Google Scholar

Bastos, A. M., Vezoli, J., and Fries, P. (2015). Communication through coherence with inter-areal delays. Curr. Opin. Neurobiol. 31, 173–180. doi: 10.1016/j.conb.2014.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Belluscio, M. A., Mizuseki, K., Schmidt, R., Kempter, R., and Buzsáki, G. (2012). Cross-frequency phase-phase coupling between θ and γ oscillations in the hippocampus. J. Neurosci. 32, 423–435. doi: 10.1523/JNEUROSCI.4122-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezaire, M. J., Raikov, I., Burk, K., Vyas, D., and Soltesz, I. (2016). Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent ca1 circuit. eLife 5, e18566. doi: 10.7554/eLife.18566.172

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezaire, M. J., and Soltesz, I. (2013). Quantitative assessment of ca1 local circuits: knowledge base for interneuron-pyramidal cell connectivity. Hippocampus 23, 751–785. doi: 10.1002/hipo.22141

PubMed Abstract | CrossRef Full Text | Google Scholar

Bjerke, I. E., Øvsthus, M., Papp, E. A., Yates, S. C., Silvestri, L., Fiorilli, J., et al. (2018). Data integration through brain atlasing: Human brain project tools and strategies. Eur. Psychiat. 50, 70–76. doi: 10.1016/j.eurpsy.2018.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Burgess, N., and O'Keefe, J. (2011). Models of place and grid cell firing and theta rhythmicity. Curr. Opin. Neurobiol. 21, 734–744. doi: 10.1016/j.conb.2011.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G. (2002). Theta oscillations in the hippocampus. Neuron 33, 325–340. doi: 10.1016/S0896-6273(02)00586-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., Anastassiou, C. A., and Koch, C. (2012). The origin of extracellular fields and currents-eeg, ecog, lfp and spikes. Nat. Rev. Neurosci. 13, 407–420. doi: 10.1038/nrn3241

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Moser, E. I. (2013). Memory, navigation and theta rhythm in the hippocampal-entorhinal system. Nat. Neurosci. 16, 130–138. doi: 10.1038/nn.3304

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Wang, X.-J. (2012). Mechanisms of gamma oscillations. Ann. Rev. Neurosci. 35, 203–225. doi: 10.1146/annurev-neuro-062111-150444

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, R. T. Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D. (2019). “Neural ordinary differential equations,” in Advances in Neural Information Processing Systems 31 (NeurIPS 2018).

Google Scholar

Chizhov, A. V. (2017). Conductance-based refractory density approach: comparison with experimental data and generalization to lognormal distribution of input current. Biol. Cyber. 111, 353–364. doi: 10.1007/s00422-017-0727-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Chizhov, A. V., and Graham, L. J. (2007). Population model of hippocampal pyramidal neurons, linking a refractory density approach to conductance-based neurons. Phys. Rev. 75, 011924. doi: 10.1103/PhysRevE.75.011924

PubMed Abstract | CrossRef Full Text | Google Scholar

Chizhov, A. V., and Graham, L. J. (2008). Efficient evaluation of neuron populations receiving colored-noise current based on a refractory density method. Phys. Rev. 77, 011910. doi: 10.1103/PhysRevE.77.011910

PubMed Abstract | CrossRef Full Text | Google Scholar

Colgin, L. L. (2015). Do slow and fast gamma rhythms correspond to distinct functional states in the hippocampal network? Brain Res. 1621, 309–315. doi: 10.1016/j.brainres.2015.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Colgin, L. L., and Moser, E. I. (2010). Gamma oscillations in the hippocampus. Physiology 25, 319–329. doi: 10.1152/physiol.00021.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cutsuridis, V., Cobb, S., and Graham, B. P. (2010). Encoding and retrieval in a model of the hippocampal ca1 microcircuit. Hippocampus 20, 423–446. doi: 10.1002/hipo.20661

PubMed Abstract | CrossRef Full Text | Google Scholar

Cutsuridis, V., and Hasselmo, M. (2012). Gabaergic contributions to gating, timing, and phase precession of hippocampal neuronal activity during theta oscillations. Hippocampus 22, 1597–1621. doi: 10.1002/hipo.21002

PubMed Abstract | CrossRef Full Text | Google Scholar

Cutsuridis, V., and Poirazi, P. (2015). A computational study on how theta modulated inhibition can account for the long temporal windows in the entorhinal-hippocampal loop. Neurobiol. Learn. Memory 120, 69–83. doi: 10.1016/j.nlm.2015.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Einevoll, G. T., Kayser, C., Logothetis, N. K., and Panzeri, S. (2013). Modelling and analysis of local field potentials for studying the function of cortical circuits. Nat. Rev. Neurosci. 14, 770–785. doi: 10.1038/nrn3599

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernández-Ruiz, A., Oliva, A., Nagy, G. A., Maurer, A. P., Berényi, A., and Buzsáki, G. (2017). Entorhinal-ca3 dual-input control of spike timing in the hippocampus by theta-gamma coupling. Neuron 93, 1213–1226.e5. doi: 10.1016/j.neuron.2017.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Fries, P. (2015). Rhythms for cognition: Communication through coherence. Neuron 88, 220–235. doi: 10.1016/j.neuron.2015.09.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuentealba, P., Klausberger, T., Karayannis, T., Suen, W. Y., Huck, J., Tomioka, R., et al. (2010). Expression of coup-tfii nuclear receptor in restricted gabaergic neuronal populations in the adult rat hippocampus. J. Neurosci. 30, 1595–1609. doi: 10.1523/JNEUROSCI.4199-09.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Geiller, T., Sadeh, S., Rolotti, S. V., Blockus, H., Vancura, B., Negrean, A., et al. (2022). Local circuit amplification of spatial selectivity in the hippocampus. Nature 601, 105–109. doi: 10.1038/s41586-021-04169-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Geiller, T., Vancura, B., Terada, S., Troullinou, E., Chavlis, S., Tsagkatakis, G., et al. (2020). Large-scale 3d two-photon imaging of molecularly identified ca1 interneuron dynamics in behaving mice. Neuron 108, 968–983.e9. doi: 10.1016/j.neuron.2020.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Harten, A., and Osher, S. (1987). Uniformly high-order accurate nonoscillatory schemes. SIAM J. Numer. Analy. 24, 279–309. doi: 10.1137/0724022

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinman, J. R., Penley, S. C., Long, L. L., Escab,í, M. A., and Chrobak, J. J. (2011). Septotemporal variation in dynamics of theta: speed and habituation. J. Neurophysiol. 105, 2675–2686. doi: 10.1152/jn.00837.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Huh, C. Y. L., Amilhon, B., Ferguson, K. A., Manseau, F., Torres-Platas, S. G., Peach, J. P., et al. (2016). Excitatory inputs determine phase-locking strength and spike-timing of ca1 stratum oriens/alveus parvalbumin and somatostatin interneurons during intrinsically generated hippocampal theta rhythm. J. Neurosci. 36, 6605–6622. doi: 10.1523/JNEUROSCI.3951-13.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaramillo, J., and Kempter, R. (2017). Phase precession: a neural code underlying episodic memory? Curr. Opin. Neurobiol. 43, 130–138. doi: 10.1016/j.conb.2017.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Joshi, A., Salib, M., Viney, T. J., Dupret, D., and Somogyi, P. (2017). Behavior-dependent activity and synaptic organization of septo-hippocampal gabaergic neurons selectively targeting the hippocampal ca3 area. Neuron 96, 1342–1357.e5. doi: 10.1016/j.neuron.2017.10.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Justus, D., Dalügge, D., Bothe, S., Fuhrmann, F., Hannes, C., Kaneko, H., et al. (2017). Glutamatergic synaptic integration of locomotion speed via septoentorhinal projections. Nat. Neurosci. 20, 16–19. doi: 10.1038/nn.4447

PubMed Abstract | CrossRef Full Text | Google Scholar

Katona, L., Lapray, D., Viney, T. J., Oulhaj, A., Borhegyi, Z., Micklem, B. R., et al. (2014). Sleep and movement differentiates actions of two types of somatostatin-expressing gabaergic interneuron in rat hippocampus. Neuron 82, 872–886. doi: 10.1016/j.neuron.2014.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Kingma, D. P., and Ba, J. (2014). Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.

Google Scholar

Klausberger, T., Magill, P. J., Márton, L. F., Roberts, J. D. B., Cobden, P. M., Buzsáki, G., et al. (2003). Brain-state- and cell-type-specific firing of hippocampal interneurons in vivo. Nature 421, 844–848. doi: 10.1038/nature01374

PubMed Abstract | CrossRef Full Text | Google Scholar

Klausberger, T., Marton, L. F., O'Neill, J., Huck, J. H. J., Dalezios, Y., Fuentealba, P., et al. (2005). Complementary roles of cholecystokinin- and parvalbumin-expressing gabaergic neurons in hippocampal network oscillations. J. Neurosci. 25, 9782–9793. doi: 10.1523/JNEUROSCI.3269-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Komendantov, A. O., Venkadesh, S., Rees, C. L., Wheeler, D. W., Hamilton, D. J., and Ascoli, G. A. (2019). Quantitative firing pattern phenotyping of hippocampal neuron types. Sci. Rep. 9, 17915. doi: 10.1038/s41598-019-52611-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Lapray, D., Lasztoczi, B., Lagler, M., Viney, T. J., Katona, L., Valenti, O., et al. (2012). Behavior-dependent specialization of identified hippocampal interneurons. Nat. Neurosci. 15, 1265–1271. doi: 10.1038/nn.3176

PubMed Abstract | CrossRef Full Text | Google Scholar

Lasztóczi, B., Tukker, J. J., Somogyi, P., and Klausberger, T. (2011). Terminal field and firing selectivity of cholecystokinin-expressing interneurons in the hippocampal ca3 area. J. Neurosci. 31, 18073–18093. doi: 10.1523/JNEUROSCI.3573-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ledberg, A., and Robbe, D. (2011). Locomotion-related oscillatory body movements at 6-12 hz modulate the hippocampal theta rhythm. PLoS ONE 6, e27575. doi: 10.1371/journal.pone.0027575

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. H., Delbruck, T., and Pfeiffer, M. (2016). Training deep spiking neural networks using backpropagation. Front. Neurosci. 10, 508. doi: 10.3389/fnins.2016.00508

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, L. L., Hinman, J. R., Chen, C.-M., Escabi, M. A., and Chrobak, J. J. (2014). Theta dynamics in rat: speed and acceleration across the septotemporal axis. PLoS ONE 9, e97987. doi: 10.1371/journal.pone.0097987

PubMed Abstract | CrossRef Full Text | Google Scholar

López-Madrona, V. J., Pérez-Montoyo, E., Álvarez Salvado, E., Moratal, D., Herreras, O., Pereda, E., et al. (2020). Different theta frameworks coexist in the rat hippocampus and are coordinated during memory-guided and novelty tasks. eLife 9, e57313. doi: 10.7554/eLife.57313.sa2

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardia, K. V., and Jupp, P. E. (1999). Directional Statistics. Wiley Series in Probability and Statistics. Hoboken, NJ, USA: John Wiley and Sons, Inc. doi: 10.1002/9780470316979

CrossRef Full Text | Google Scholar

Mizuseki, K., Sirota, A., Pastalkova, E., and Buzsáki, G. (2009). Theta oscillations provide temporal windows for local circuit computation in the entorhinal-hippocampal loop. Neuron 64, 267–280. doi: 10.1016/j.neuron.2009.08.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Moradi, K., Aldarraji, Z., Luthra, M., Madison, G. P., and Ascoli, G. A. (2022). Normalized unitary synaptic signaling of the hippocampus and entorhinal cortex predicted by deep learning of experimental recordings. Commun. Biol. 5, 418. doi: 10.1038/s42003-022-03329-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Mysin, I. (2021). A model of the ca1 field rhythms. eNeuro 8, ENEURO.0192-21.2021. doi: 10.1523/ENEURO.0192-21.2021

PubMed Abstract | CrossRef Full Text | Google Scholar

Mysin, I., and Shubina, L. (2022). From mechanisms to functions: The role of theta and gamma coherence in the intrahippocampal circuits. Hippocampus 32, 342–358. doi: 10.1002/hipo.23410

PubMed Abstract | CrossRef Full Text | Google Scholar

Mysin, I. E., Kitchigina, V. F., and Kazanovich, Y. B. (2019). Phase relations of theta oscillations in a computer model of the hippocampal ca1 field: Key role of schaffer collaterals. Neural Netw. 116, 119–138. doi: 10.1016/j.neunet.2019.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Nuñez, A., and Buño, W. (2021). The theta rhythm of the hippocampus: from neuronal and circuit mechanisms to behavior. Front. Cell. Neurosci. 15, 649262. doi: 10.3389/fncel.2021.649262

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramezanian-Panahi, M., Abrevaya, G., Gagnon-Audet, J.-C., Voleti, V., Rish, I., and Dumas, G. (2022). Generative models of brain dynamics. Front. Artif. Intell. 5, 807406. doi: 10.3389/frai.2022.807406

PubMed Abstract | CrossRef Full Text | Google Scholar

Royer, S., Zemelman, B. V., Losonczy, A., Kim, J., Chance, F., Magee, J. C., et al. (2012). Control of timing, rate and bursts of hippocampal place cells by dendritic and somatic inhibition. Nat. Neurosci. 15, 769–775. doi: 10.1038/nn.3077

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakalar, E., Klausberger, T., and Lasztóczi, B. (2022). Neurogliaform cells dynamically decouple neuronal synchrony between brain areas. Science 377, 324–328. doi: 10.1126/science.abo3355

PubMed Abstract | CrossRef Full Text | Google Scholar

Saraga, F., Wu, C. P., Zhang, L., and Skinner, F. K. (2003). Active dendrites and spike propagation in multi-compartment models of oriens-lacunosum/moleculare hippocampal interneurons. J. Physiol. 552, 673–689. doi: 10.1113/jphysiol.2003.046177

PubMed Abstract | CrossRef Full Text | Google Scholar

Somogyi, P., Katona, L., Klausberger, T., Lasztóczi, B., and Viney, T. J. (2014). Temporal redistribution of inhibition over neuronal subcellular domains underlies state-dependent rhythmic change of excitability in the hippocampus. Philos. Trans. R Soc. London. B, Biol. Sci. 369, 20120518. doi: 10.1098/rstb.2012.0518

PubMed Abstract | CrossRef Full Text | Google Scholar

Somogyi, P., and Klausberger, T. (2005). Defined types of cortical interneurone structure space and spike timing in the hippocampus. J. Physiol. 562, 9–26. doi: 10.1113/jphysiol.2004.078915

PubMed Abstract | CrossRef Full Text | Google Scholar

Stimberg, M., Brette, R., and Goodman, D. F. (2019). Brian 2, an intuitive and efficient neural simulator. eLife 8, e47314. doi: 10.7554/eLife.47314.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y., Zhang, L., and Schaeffer, H. (2020). “NeuPDE: Neural network based ordinary and partial differential equations for modeling time-dependent data,” in Mathematical and Scientific Machine Learning (PMLR) 352–372.

Google Scholar

Tsodyks, M., Pawelzik, K., and Markram, H. (1998). Neural networks with dynamic synapses. Neural Comput. 10, 821–835. doi: 10.1162/089976698300017502

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsodyks, M. V., and Markram, H. (1997). The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability. Proc. Nat. Acad. Sci. USA. 94, 719–723. doi: 10.1073/pnas.94.2.719

PubMed Abstract | CrossRef Full Text | Google Scholar

Unal, G., Crump, M. G., Viney, T. J., Éltes, T., Katona, L., Klausberger, T., et al. (2018). Spatio-temporal specialization of gabaergic septo-hippocampal neurons for rhythmic network activity. Brain Struc. Funct. 223, 2409–2432. doi: 10.1007/s00429-018-1626-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Varga, C., Golshani, P., and Soltesz, I. (2012). Frequency-invariant temporal ordering of interneuronal discharges during hippocampal oscillations in awake mice. Proc. Nat. Acad. Sci. USA. 109, E2726-E2734. doi: 10.1073/pnas.1210929109

PubMed Abstract | CrossRef Full Text | Google Scholar

Varga, C., Oijala, M., Lish, J., Szabo, G. G., Bezaire, M., Marchionni, I., et al. (2014). Functional fission of parvalbumin interneuron classes during fast network events. eLife 3, e04006. doi: 10.7554/eLife.04006.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Venkadesh, S., Komendantov, A. O., Wheeler, D. W., Hamilton, D. J., and Ascoli, G. A. (2019). Simple models of quantitative firing phenotypes in hippocampal neurons: Comprehensive coverage of intrinsic diversity. PLoS Comput. Biol. 15, e1007462. doi: 10.1371/journal.pcbi.1007462

PubMed Abstract | CrossRef Full Text | Google Scholar

Viney, T. J., Salib, M., Joshi, A., Unal, G., Berry, N., and Somogyi, P. (2018). Shared rhythmic subcortical gabaergic input to the entorhinal cortex and presubiculum. eLife 7, e34395. doi: 10.7554/eLife.34395.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinogradova, O. S. (1995). Expression, control, and probable functional significance of the neuronal theta-rhythm. Progr. Neurobiol. 45, 523–583. doi: 10.1016/0301-0082(94)00051-I

PubMed Abstract | CrossRef Full Text | Google Scholar

Wheeler, D. W., White, C. M., Rees, C. L., Komendantov, A. O., Hamilton, D. J., and Ascoli, G. A. (2015). Hippocampome.org: a knowledge base of neuron types in the rodent hippocampus. eLife 4, e09960. doi: 10.7554/eLife.09960

PubMed Abstract | CrossRef Full Text | Google Scholar

Wunderlich, T. C., and Pehle, C. (2021). Event-based backpropagation can compute exact gradients for spiking neural networks. Scient. Rep. 11, 12829. doi: 10.1038/s41598-021-91786-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Zenke, F., and Vogels, T. P. (2021). The remarkable robustness of surrogate gradient learning for instilling complex function in spiking neural networks. Neur. Comput. 33, 899–925. doi: 10.1162/neco_a_01367

PubMed Abstract | CrossRef Full Text | Google Scholar

Zutshi, I., Valero, M., Fernández-Ruiz, A., and Buzsáki, G. (2022). Extrinsic control and intrinsic computation in the hippocampal ca1 circuit. Neuron 110, 658–673.e5. doi: 10.1016/j.neuron.2021.11.015

PubMed Abstract | CrossRef Full Text | Google Scholar