
==== Front
eLife
Elife
eLife
eLife
2050-084X
eLife Sciences Publications, Ltd

38354040
87356
10.7554/eLife.87356
Research Article
Computational and Systems Biology
Neuroscience
A dynamical computational model of theta generation in hippocampal circuits to study theta-gamma oscillations during neurostimulation
Vardalakis Nikolaos https://orcid.org/0000-0002-5436-4091
12
Aussel Amélie https://orcid.org/0000-0003-0498-2905
123
Rougier Nicolas P https://orcid.org/0000-0002-6972-589X
123
Wagner Fabien B https://orcid.org/0000-0002-9582-6109
fabien.wagner@u-bordeaux.fr
1
1 https://ror.org/057qpr032 University of Bordeaux, CNRS, IMN Bordeaux France
2 https://ror.org/02kvxyf05 University of Bordeaux, INRIA, IMN Bordeaux France
3 https://ror.org/054qv7y42 University of Bordeaux, CNRS, Bordeaux INP Talence France
Toth Katalin Reviewing Editor https://ror.org/03c4mmv16 University of Ottawa Canada

Colgin Laura L Senior Editor https://ror.org/00hj54h04 University of Texas at Austin United States

14 2 2024
2024
12 RP8735624 3 2023
This manuscript was published as a preprint.27 3 2023

This manuscript was published as a reviewed preprint.26 6 2023

The reviewed preprint was revised.11 1 2024

© 2023, Vardalakis et al
2023
Vardalakis et al
https://creativecommons.org/licenses/by/4.0/ This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Neurostimulation of the hippocampal formation has shown promising results for modulating memory but the underlying mechanisms remain unclear. In particular, the effects on hippocampal theta-nested gamma oscillations and theta phase reset, which are both crucial for memory processes, are unknown. Moreover, these effects cannot be investigated using current computational models, which consider theta oscillations with a fixed amplitude and phase velocity. Here, we developed a novel computational model that includes the medial septum, represented as a set of abstract Kuramoto oscillators producing a dynamical theta rhythm with phase reset, and the hippocampal formation, composed of biophysically realistic neurons and able to generate theta-nested gamma oscillations under theta drive. We showed that, for theta inputs just below the threshold to induce self-sustained theta-nested gamma oscillations, a single stimulation pulse could switch the network behavior from non-oscillatory to a state producing sustained oscillations. Next, we demonstrated that, for a weaker theta input, pulse train stimulation at the theta frequency could transiently restore seemingly physiological oscillations. Importantly, the presence of phase reset influenced whether these two effects depended on the phase at which stimulation onset was delivered, which has practical implications for designing neurostimulation protocols that are triggered by the phase of ongoing theta oscillations. This novel model opens new avenues for studying the effects of neurostimulation on the hippocampal formation. Furthermore, our hybrid approach that combines different levels of abstraction could be extended in future work to other neural circuits that produce dynamical brain rhythms.

computational model
hippocampus
theta-gamma oscillations
medial septum
Kuramoto oscillators
Research organism

None
http://dx.doi.org/10.13039/501100009468 Conseil Régional Aquitaine Bordeaux Neurocampus chair of excellence Wagner Fabien B http://dx.doi.org/10.13039/501100006251 Université de Bordeaux Bordeaux Neurocampus chair of excellence Wagner Fabien B http://dx.doi.org/10.13039/501100000781 European Research Council ERC Starting Grant #101040391 (MEMOPROSTHETICS project) Wagner Fabien B The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.Author impact statementCombining biophysically realistic hippocampal neurons with coupled phase oscillators that represent theta inputs enables probing the phase-dependent effects of neurostimulation on theta-nested gamma oscillations in normal and pathological states.
publishing-routeprc
==== Body
pmcIntroduction

Neurostimulation methods have emerged as promising therapeutic modalities to restore neurological functions in a broad range of motor and cognitive disorders (Gupta et al., 2023). In the context of learning and memory, deep brain stimulation (DBS) of the entorhinal area or hippocampus has been shown to either enhance (Suthana et al., 2012; Suthana and Fried, 2014; Titiz et al., 2017; Jun et al., 2019) or disrupt (Jacobs et al., 2016; Goyal et al., 2018; Lozano et al., 2016) memory encoding. These conflicting results may originate from differences in experimental protocols and from a poor understanding of their biophysical underpinnings. Among such mechanisms, the involvement of hippocampal theta oscillations (4–12 Hz) and their interactions with higher-frequency gamma oscillations (30–120 Hz) in memory-related processes has been reported in multiple studies (Lisman et al., 2005; de Almeida et al., 2007; Lega et al., 2012; Lin et al., 2017; Malkov et al., 2022; Abbaspoor et al., 2023). Moreover, the modulation of gamma oscillations by the phase of theta oscillations in hippocampal circuits, a phenomenon termed theta-gamma phase-amplitude coupling (PAC), correlates with the efficacy of memory encoding and retrieval (Jensen and Colgin, 2007; Tort et al., 2009; Canolty and Knight, 2010; Axmacher et al., 2010; Fell and Axmacher, 2011; Lisman and Jensen, 2013; Lega et al., 2016). Experimental and computational work on the coupling between oscillatory rhythms has indicated that it originates from different neural architectures and correlates with a range of behavioral and cognitive functions, enabling the long-range synchronization of cortical areas and facilitating multi-item encoding in the context of memory (Hyafil et al., 2015).

Neurostimulation protocols that affect these rhythms, such as theta burst stimulation (Titiz et al., 2017), have also been shown to optimally induce long-term potentiation (LTP) (Larson and Munkácsy, 2015). Another potential mechanism underlying the effect of hippocampal neurostimulation might be the reset of the phase of theta oscillations in response to exogenous inputs, such as a novel sensory input or a pulse of electrical stimulation applied to the fornix or perforant path (Buño et al., 1978; Williams and Givens, 2003). Theta phase reset is known to facilitate LTP (McCartney et al., 2004) and naturally occurs during both encoding and retrieval of associative memories (Kota et al., 2020). In this context, the design of a computational model that replicates memory-related theta-gamma oscillations and theta phase reset is of uttermost importance to investigate the effects of electrical stimulation on the hippocampal formation and possibly optimize neurostimulation protocols for memory improvement.

Models of memory-related theta-gamma oscillations in the hippocampal formation have been developed across different resolution levels, ranging from abstract mean-field approaches (Traub et al., 1997; Onslow et al., 2014; Segneri et al., 2020) to biophysically realistic conductance-based models (Lundqvist et al., 2006; Herman et al., 2013; Aussel et al., 2018). Neural masses, which represent the mean activity of a neuronal population, can generate gamma oscillations through reciprocal interactions between an excitatory and inhibitory population (Traub et al., 1997; Onslow et al., 2014), or even using a self-projecting inhibitory population (Segneri et al., 2020). Under excitatory oscillatory input at theta frequencies, these models are capable of generating theta-nested gamma oscillations. Similarly, these oscillations can be observed in more complex models of the hippocampal formation composed of single-compartment excitatory and inhibitory neurons connected through conductance-based synapses following the Hodgkin-Huxley formalism (Hodgkin and Huxley, 1952), and driven by a fixed oscillatory theta input (Aussel et al., 2018). Theta-nested gamma oscillations also appear in multi-compartment models of prefrontal cortex activity during memory retrieval (Lundqvist et al., 2006; Herman et al., 2013). Finally, several models have investigated the functional link between the hippocampal theta rhythm and memory processes, showing that encoding and retrieval occur at different phases within each theta cycle through phasic changes in neuronal activity and LTP (e.g. Hasselmo et al., 2002; Cutsuridis et al., 2010).

In terms of neurostimulation, most computational work has focused on the mechanisms underlying DBS of the basal ganglia for motor disorders such as Parkinson’s disease (Rubin and Terman, 2004; Pirini et al., 2009; Mina et al., 2013; Ebert et al., 2014), peripheral nerve stimulation (Rattay et al., 2003; Kipping and Nogueira, 2022), spinal cord stimulation (Rattay et al., 2000; Capogrosso et al., 2013), or has remained generic (Basu et al., 2018). However, models investigating neurostimulation of hippocampal circuits are scarce (Hendrickson et al., 2016; Bingham et al., 2018) and do not take into account the effects on theta-gamma oscillations. For example, a detailed multicompartment model of the rat dentate gyrus was able to replicate experimentally recorded local field potentials induced by different electrode placements and pulse amplitudes during stimulation of the perforant path (Bingham et al., 2018). To model the impact of neurostimulation on neuronal oscillations, a more abstract formalism based on Kuramoto phase oscillators (Kuramoto, 1984) has been introduced in the context of Parkinson’s disease and essential tremor, which enabled the design of novel neurostimulation paradigms that enhance or disrupt neuronal synchrony in basal ganglia circuits (Tass, 2003; Ebert et al., 2014; Asllani et al., 2018; Weerasinghe et al., 2019).

To our knowledge, there is currently no model of the hippocampal formation that is able to replicate both theta-gamma PAC and theta phase reset during neurostimulation. To investigate the effects of neurostimulation on hippocampal circuits while taking into account these two mechanisms, we modified an existing biophysical model of the hippocampal formation (Aussel et al., 2018). In this original model, however, the theta rhythm was considered as an oscillatory input of fixed amplitude and phase velocity, which is inconsistent with theta phase reset. To circumvent this limitation, we combined this biophysical model with abstract Kuramoto oscillators that acted as a dynamical source of theta rhythm, thereby modeling medial septum inputs to the hippocampal formation. This new hybrid dynamical model could generate both theta-nested gamma oscillations and theta phase reset, following a particular phase response curve (PRC) inspired by experimental literature (Lengyel et al., 2005; Akam et al., 2012; Torben-Nielsen et al., 2010).

We then leveraged this model to explore the effect of single-pulse and pulse-train stimulation on theta-gamma oscillations. In the absence of theta input from the medial septum, single-pulse stimulation produced a transient effect consisting of one or several bursts of activity, depending on stimulation amplitude. The presence of multiple bursts depended on the single-cell calcium dynamics and M-type potassium adaptation current. In the presence of weak theta input, designed to mimic a pathological state that impairs theta-gamma oscillations, single-pulse stimulation could produce long-lasting or even persistent activity by switching the network to a highly synchronized state characterized by theta-nested gamma oscillations. When phase reset was not included in the model, this effect was more pronounced when the stimulation pulse was delivered at the peak of the theta rhythm. However, when strong theta reset was considered, the phase at which stimulation was delivered did not influence the outcome. In the presence of an even weaker theta input, mimicking a pathological state that completely abolishes theta-gamma oscillations, only pulse train stimulation could restore physiological theta-gamma oscillations during the stimulation period. As for the previous results, this effect was phase-dependent only when theta reset was not included.

These results provide a new framework to interpret neurostimulation interventions that interfere with hippocampal oscillations and aim at improving memory function. It can be further extended to investigate the effects of more complex neurostimulation protocols, and the impact of stimulation location and amplitude on the observed network dynamics.

Results

A computational model of the hippocampal formation with dynamical theta input

We first developed a computational model of hippocampal circuits, able to generate both theta-nested gamma oscillations and theta phase reset. To achieve this, we combined an existing conductance-based model of the hippocampal formation (Aussel et al., 2018) with a set of Kuramoto phase oscillators (Kuramoto, 1984) that were used to model the dynamical theta input originating from pacemaker neurons in the medial septum (Wang, 2002, Figure 1A). Such ensembles of oscillators are designed to exhibit a strong phase reset of their collective rhythm in response to a perturbation (Levnajić and Pikovsky, 2010). The hippocampal model contained excitatory and inhibitory neuronal populations in the entorhinal cortex (EC), dentate gyrus (DG), CA3 and CA1 fields within a coronal slice of the human hippocampus (Figure 1B–C). The output of the hippocampal formation (i.e. the activity of CA1 pyramidal neurons) was provided as an input to the Kuramoto oscillators, simulating the hippocampal-septal projections through the fornix (Williams and Givens, 2003; Nuñez and Buño, 2021; Takeuchi et al., 2021). The oscillations produced by the collective behavior of Kuramoto oscillators represented a population average of their activity: highly synchronized and desynchronized states generated respectively high and low amplitude theta oscillations (Figure 1D). The number of oscillators was chosen so that their synchronization level (i.e. their order parameter) and their frequency distribution were sufficiently close to their asymptotic behavior for a large number of oscillators (Figure 1—figure supplement 1).

Figure 1. Dynamical computational model of the hippocampal formation and medial septum oscillatory drive.

(A) Anatomical representation of the neuronal types and interconnections within and between the medial septum and the hippocampal formation (EC: entorhinal cortex, DG: dentate gyrus, CA3 and CA1 fields of the hippocampus). (B) Simplified anatomy of the hippocampal formation, modeled as a 15-mm-thick cylindrical slice, with spatially segregated excitatory and inhibitory neurons (blue: excitatory neurons, consisting of granule cells in DG and pyramidal cells in other areas; red: inhibitory basket cells). (C) Model architecture and connectivity. Each area is comprised of one excitatory and one inhibitory neuronal population. Theta drive is provided through input from the medial septum, which is modeled as a set of 250 Kuramoto oscillators and receives feedback connections from CA1. Electrical stimulation Is modeled as an intracellular current affecting both excitatory and inhibitory populations in the targeted area (shown here for CA1). (D) Illustration of Kuramoto oscillators with two different levels of synchronization. Each dot represents one oscillator, its position on the circle indicates its phase and its color its angular velocity. Higher synchronization corresponds to a clustering of the dots around a similar phase. ‘r’ indicates the order parameter, which is a measure of synchronization.

Figure 1—figure supplement 1. Influence of the number of Kuramoto oscillators on their collective dynamics.

Left: influence of the number of Kuramoto oscillators (N) on the convergence over time of the order parameter (r), which is a measure of synchronization of the set of oscillators. The above results were obtained using a synchronization ratio k/N of 25. Upon initialization, the oscillators’ natural frequencies (ωi) were normally distributed with a center frequency (f0) of 6 Hz. Initial phases were uniformly distributed around the unit circle (U∼[0,2π]). Right: histograms of the natural frequencies (ωi), together with the distribution mean (μ) and standard error of the mean (σ/N). As the number of oscillators increases, the mean frequency appears closer to the desired center frequency of 6 Hz, and the standard error of the mean is reduced. Overall, increasing the number of oscillators past 250 did not yield any substantial change in the convergence of the order parameter or the distribution of natural frequencies. This number of oscillators was therefore chosen for all subsequent analyses.

Biologically, GABAergic neurons from the medial septum project to the EC, CA3, and CA1 fields of the hippocampus (Tóth et al., 1993; Hajós et al., 2004; Manseau et al., 2008; Hangya et al., 2009; Unal et al., 2015; Müller and Remy, 2018). Although the respective roles of these different projections are not fully understood, previous computational studies have suggested that the direct projection from the medial septum to CA1 is not essential for the production of theta in CA1 microcircuits (Mysin et al., 2019). Since our modeling of the medial septum is only used to generate a dynamic theta rhythm, we opted for a simplified representation where the medial septum projects only to the EC, which in turn drives the different subfields of the hippocampus. In our model, Kuramoto oscillators are therefore connected to the EC neurons and they receive projections from CA1 neurons (see Materials and methods for more details).

From a conceptual point of view, our model is thus composed of excitatory-inhibitory (E-I) circuits connected in series, with a feedback loop going through a population of coupled phase oscillators. In the next sections, we first describe the generation of gamma oscillations by individual E-I circuits (Figure 2), and illustrate their behavior when driven by an oscillatory input such as theta oscillations (Figure 3). We then present a thorough characterization of the effects of theta input and stimulation amplitude on theta-nested gamma oscillations (Figure 4 and Figure 5). Finally, we present some results on the effects of neurostimulation protocols for restoring theta-nested gamma oscillations in pathological states (Figure 6 and Figure 7).

Figure 2. Emergence of gamma oscillations in coupled excitatory-inhibitory populations under ramping input to both populations.

(A) Two coupled populations of excitatory pyramidal neurons (NE = 1000) and inhibitory interneurons (NI = 100) are driven by a ramping current input (0 nA to 1 nA) for 5 s. As the input becomes stronger, oscillations start to emerge (shaded green area), driven by the interactions between excitatory and inhibitory populations. The green inset shows the raster plot (neuronal spikes across time) of the two populations during the green shaded period (red for inhibitory; blue for excitatory). When the input becomes sufficiently strong (shaded magenta area), the populations become highly synchronized and produce oscillations in the gamma range (at approximately 50 Hz). The spectrogram (bottom panel) shows the power of the instantaneous firing rate of the pyramidal population as a function of time and frequency. It reveals the presence of gamma oscillations that emerge around 2 s and increase in frequency until 4 s, when they settle at approximately 60 Hz. (B) Similar depiction as in panel A. with the pyramidal-interneuronal populations decoupled. The absence of coupling leads to the abolition of gamma oscillations, each cell spiking activity being driven by its own inputs and intrinsic properties.

Figure 2—figure supplement 1. Emergence of gamma oscillations in coupled excitatory-inhibitory populations under ramping input to the excitatory population.

Similar representation as in Figure 2, but with the input provided only to the excitatory population. All conclusions remain the same. In addition, the inhibitory population does not show any spiking activity in the decoupled case.

Figure 2—figure supplement 2. Cell-intrinsic spiking activity in decoupled excitatory and inhibitory populations under ramping input.

(A) Input-Frequency (I–F) curves for excitatory cells (left panel; pyramidal neurons with ICAN) and inhibitory cells (right panel; interneurons, fast-spiking) used in the model. Above a certain tonic input (around 0.35 nA for excitatory and 0.1 nA for inhibitory neurons), neurons can spike in the gamma range. (B) Raster plot showing the spiking activity of excitatory (blue,NE = 1000) and inhibitory (red,NI = 100) neurons in decoupled populations under ramping input (top trace) and in the absence of noise in the membrane potential. Despite random initial conditions across neurons, oscillations emerge in both populations due to the intrinsic properties of the cells, with a frequency that is predicted by the respective I-F curves (panel A.). (C) Similar representation as panel B. but with the addition of stochastic noise in the membrane potential of each neuron. The presence of noise disrupts the emergence of oscillations in these decoupled populations.

Figure 2—figure supplement 3. Neuronal entrainment of fast-spiking interneurons by pulsed optogenetic stimulation in the high gamma range.

All panels show the effects of pulsed optogenetic stimulation of specific cell types (A. excitatory only, B. inhibitory only, C. both excitatory and inhibitory cells) at different frequencies, similar to the experiments reported by Cardin and colleagues (Cardin et al., 2009). Pulsed stimulation was delivered continuously with a frequency ranging from 10 to 200 Hz. The relative power was calculated similarly to Cardin et al., 2009. Specifically, for each stimulation frequency, we computed first the power spectral density of the population firing rates using Welch’s method, and then the ratio between the power within a 10-Hz band centered around the stimulation frequency and the total power across all frequencies. (A) Pulsed stimulation of excitatory cells shows a higher relative power in the stimulated neurons at low frequencies (20 and 40 Hz). (B) Pulsed stimulation of inhibitory cells reveals neuronal entrainment, that is, a peak in their relative power, in the high gamma range (80-100 Hz). (C) Simultaneous pulsed stimulation of both populations reveals a strong entrainment of inhibitory cells in the high gamma range, and a more modest entrainment of excitatory cells at low frequencies.

Figure 3. Theta-nested gamma oscillations and theta phase reset in response to a single stimulation pulse.

(A) Representative example of the network behavior, spontaneously producing theta-nested gamma oscillations and characterized by a reset of the theta phase following a single stimulation pulse (vertical grey line, applied here in CA1 at a theta phase of π/2, that is in the middle of the descending slope). Top to bottom: theta rhythm originating from the medial septum and provided as an input to the EC (fθ: mean oscillation frequency), instantaneous phase of the theta rhythm, raster plots indicating the spiking activity of CA1 excitatory (blue) and inhibitory (red) neurons (μE and μI: mean firing rates within the shaded area), average population firing rates in CA1 (computed as a windowed moving average with a sliding window of 100ms with 99% overlap), spectrograms for each CA1 population (windowed short-time fast Fourier transform using a Hann sliding window: 100ms with 99% overlap). Spectrograms show gamma oscillations (around 60 Hz) modulated by the underlying theta rhythm (∼4 Hz), indicating theta-gamma PAC. Theta phase reset after stimulation is associated with a rebound of spiking activity and theta-nested gamma oscillations. (B) Power spectral densities of the CA1 firing rates. Theta peaks are found at 4 Hz for excitatory and inhibitory cells. Gamma activity is located between 40 and 80 Hz. (C) PAC as a function of theta phase and gamma frequency. The polar plot represents the amplitude of gamma oscillations (averaged across all theta cycles, see Materials and methods) at each phase of theta (theta range: 3–9 Hz, phase indicated as angular coordinate) and for different gamma frequencies (radial coordinate, binned in 10 Hz ranges), indicating that gamma oscillations between 40 and 80 Hz occur preferentially around the peak of theta. The MI gives an overall quantification of how the phase of low-frequency oscillations (3–9 Hz) modulates the amplitude of higher-frequency oscillations (40–80 Hz) (see Materials and methods and Figure 3—figure supplement 2). (D) PRC in response to a single stimulation pulse applied in CA1 at various phases of the ongoing theta rhythm and for various stimulation amplitudes (color-coded). The phase difference (left y-axis) shows the theta phase induced by the stimulation pulse (computed 2.5ms after the pulse), compared to the phase computed at the same time in a scenario without stimulation. Positive and negative phase differences indicate phase advances and delays, respectively. The grey trace shows the normalized amplitude of theta (right y axis) for different phases, used to indicate the peak and trough of the rhythm. Stimulation applied in the ascending slope of theta ([−π,0]) produced a phase advance and accelerated the rhythm towards its peak (0 radians). Conversely, stimulation during the descending slope ([0,π]) produced a phase delay that slowed down the rhythm. Higher stimulation amplitudes yielded a stronger effect.

Figure 3—figure supplement 1. Whole-network behavior producing theta-nested gamma oscillations and theta phase reset in response to a single stimulation pulse.

This figure shows the behavior of all areas of the hippocampal formation during the simulation illustrated in Figure 3A. All areas display bursts of gamma activity that occur almost synchronously at the peak of each theta cycle in both excitatory and inhibitory neuronal populations (μE and:μI mean firing rates within the shaded area). The CA1 excitatory population firing rate is supplied back to the medial septum. A single stimulation pulse (vertical line and grey triangle) delivered to both CA1 populations at a phase of π/2 (i.e. halfway through the descending slope) resets the phase of the theta rhythm according to the PRC shown in Figure 3D.

Figure 3—figure supplement 2. Quantification of PAC with and without noise.

(A) Quantifying PAC in the absence of noise produced inaccurate identification of the coupled frequency bands, due to the complete absence of oscillations at some frequencies. All analyses are based on the CA1 firing rates (top traces) during a representative simulation. Power spectral densities of these firing rates (left) indicate that some frequencies have a power of 0. PAC of the excitatory population was assessed using two graphical representations, the polar plot (middle) and comodulogram (right), and quantified using the MI. The comodulogram was calculated by computing the MI across 80% overlapping 1-Hz frequency bands in the theta range and across 90% overlapping 10 Hz frequency bands in the gamma range and subsequently plotted as a heat map. In the absence of noise, a slow theta frequency centered around 5 Hz is found to modulate a broad range of gamma frequencies between 40 and 100 Hz. The value indicated on the comodulogram indicates the average MI in the 3-9 Hz theta range and 40-80 Hz gamma range. As in Figure 3, the polar plot represents the amplitude of gamma oscillations (averaged across all theta cycles) at each phase of theta (theta range: 3–9 Hz, phase indicated as angular coordinate) and for different gamma frequencies (radial coordinate, binned in 1 Hz ranges). (B) Adding uniform noise to the firing rate (with an amplitude ranging between 15 and 25% of the maximum firing rate) improved the identification of the coupled frequency bands. In this case, the slower theta frequency centered around 5 Hz modulates a gamma band located between 45 and 75 Hz.

Figure 3—figure supplement 3. Network behavior generated by Kuramoto oscillators with non-physiological phase response functions.

Each panel is similar to Figure 3A, but with a different offset added to the phase response function of the Kuramoto oscillators (see methods, Equation 4). The center frequency was set to 6 Hz in all of these simulations. Overall, theta oscillations in these cases are less sinusoidal and show more abrupt phase changes than in the physiological case. (A) A phase offset of −π/2 leads to an overall theta oscillation of 4 Hz, with a second peak following the main theta peak. (B). A phase offset of +π/2 reduces the peak of theta, resetting the rhythm to the middle of the ascending phase. (C) A phase offset of π or −π leads to the CA1 output resetting the theta rhythm to the trough of theta.

Figure 3—figure supplement 4. A neural mass model of coupled excitatory and inhibitory neurons driven by Kuramoto oscillators generates theta-nested gamma oscillations and theta phase reset.

(A) Two coupled neural masses (one excitatory and one inhibitory) driven by Kuramoto oscillators, which represent a dynamical oscillatory drive in the theta range, were used to implement a neural mass equivalent to our conductance-based model represented in Figure 1. Neural masses were modeled using the Wilson-Cowan formalism, with parameters adapted from Onslow et al., 2014 (WEE = 4.8, WEI = WIE = 4, WII = 0). (B) The normalized population firing rates exhibit theta-nested gamma oscillations (middle and bottom panels) in response to the dynamic theta rhythm (top panel). A stimulation pulse delivered at the descending phase of the rhythm to both populations (marked by the inverted red triangle) produces a robust theta phase reset, similarly to Figure 3A.

Figure 4. CAN and adaptation (M) currents modulate neuronal responses to single-pulse stimulation in the absence of theta input.

(A) Network response to single-pulse stimulation in the absence of medial septum input. Stimulation (grey vertical line, 10 nA, applied in CA1) induced an instantaneous burst of activity lasting about 20ms in both excitatory and inhibitory CA1 neurons, followed by a secondary burst approximately 200 ms later (raster plot and firing rate traces), associated with specific CAN- and M- currents dynamics (bottom traces, illustrated for a representative CA1 excitatory neuron). Positive (IM) and negative (ICAN) currents indicate respectively a hyperpolarizing and depolarizing effect on the cell membrane potential. (B) Similar representation as in A., but in the absence of CAN channel. Stimulation induced only a single burst of activity, indicating that CAN channels are necessary to observe a rebound of activity. (C) Number of bursts in CA1 spiking activity following a single stimulation pulse at various amplitudes (x-axis), shown both in the presence and absence of CAN channels in excitatory neurons. The absence of the CAN current leads to the abolition of the second burst, irrespective of stimulation amplitude.

Figure 5. Medial septum oscillatory drive and stimulation amplitude govern the steady-state response to single-pulse stimulation.

(A-C) Network responses to single-pulse stimulation (vertical line) under medial septum input (with phase reset), shown for different amplitudes of the medial septum oscillatory drive (A-C: increasing oscillator amplitudes). (A) Low oscillatory input: stimulation induces only two bursts of spiking activity as in Figure 4. (B) Medium oscillatory input: a single stimulation pulse switches network behavior from no activity to sustained oscillations driven by the medial septum. (C) Higher oscillatory input: theta drive is capable of inducing self-sustained theta-nested gamma oscillations. In this case, stimulation is delivered at the peak of theta oscillations and does not show a pronounced effect on theta-gamma oscillations. (D) Steady-state response to single-pulse stimulation as a function of medial septum oscillatory input (x-axis) and stimulation amplitude (y-axis), characterized by three metrics: theta power (3–9 Hz), gamma power (40–80 Hz), and PAC (quantified using the MI). White dots: parameter combinations corresponding to panels A-C.

Figure 5—figure supplement 1. Whole-network behavior in response to single-pulse stimulation under low theta oscillatory drive.

This figure shows the behavior of all areas of the hippocampal formation during the simulation illustrated in Figure 5A. A single stimulation pulse delivered to both CA1 populations induces two bursts of spiking activity that propagate only from CA1 to EC.

Figure 5—figure supplement 2. Whole-network behavior in response to single-pulse stimulation under medium theta oscillatory drive.

This figure shows the behavior of all areas of the hippocampal formation during the simulation illustrated in Figure 5B. A single stimulation pulse switches network behavior from no activity to sustained oscillations driven by the medial septum.

Figure 5—figure supplement 3. Whole-network behavior in response to single-pulse stimulation under higher theta oscillatory drive.

This figure shows the behavior of all areas of the hippocampal formation during the simulation illustrated in Figure 5C. Theta oscillatory drive is enough to induce self-sustained theta-nested gamma oscillations. In this case, stimulation is delivered at the peak of theta and does not induce a pronounced effect on theta-gamma oscillations.

Figure 5—figure supplement 4. CAN currents are necessary for the production of self-sustained theta-gamma oscillations in response to single-pulse stimulation.

(A) Same as Figure 5B. (B) Similar simulation as panel A., but without the presence of CAN currents in the EC, CA3, and CA1 fields of the hippocampus. Removing CAN currents from the model abolishes self-sustained theta-nested gamma oscillations in response to a single stimulation pulse (for the parameters represented in Figure 5, point B).

Figure 5—figure supplement 5. Synchronization level of the oscillators minimally affects the production of theta-nested gamma oscillations.

(A, B) Network responses to single-pulse stimulation (vertical line) under medial septum input (with phase reset), shown for different synchronization parameter ratios of the Kuramoto oscillators (A,B: k/N of 5 and 35 respectively; stimulation amplitude: 7.0 nA). Panel B. indicates that the model can produce oscillations without external stimulation when the synchronization level is higher. (C) Steady-state response to single-pulse stimulation as a function of medial septum oscillatory input (x-axis) and synchronization parameter ratio (k/N, y-axis), characterized by three metrics: theta power (3–9 Hz), gamma power (40–80 Hz), and PAC (quantified using the MI). White dots: parameter combinations corresponding to panels A and B.

Figure 6. Single-pulse stimulation phase differentially affects network responses depending on the presence of theta phase reset.

All results shown here were obtained for parameters from Figure 5B (theta oscillation amplitude: 0.13 nA, stimulation amplitude: 7.0 nA). A single stimulation pulse was delivered at the peak (A, B) or trough (C, D) of the underlying theta rhythm, either in the presence (A, C) or absence (B, D) of theta phase reset. With phase reset, both peak and trough stimulation switch network behavior from no activity to sustained oscillations. Without phase reset, only peak stimulation can induce sustained oscillations. (E). Quantification of theta power, gamma power, and PAC (measured using the MI) in CA1 excitatory (blue) and inhibitory (red) populations in all four cases (metrics are computed in the shaded areas of panels A-D).

Figure 7. Pulse train stimulation restores theta-nested gamma oscillations depending on stimulation timing and theta phase reset.

All results shown here were obtained for parameters from Figure 5A (theta oscillation amplitude: 0.05 nA; stimulation amplitude: 7.0 nA). Representations are similar to Figure 6, with the difference that stimulation consisted of a pulse train delivered at 6 Hz for a duration of 2 s (individual pulses indicated by grey dots, first pulse by a triangle). The pulse train was delivered at the peak (A, B) or trough (C, D) of the underlying theta rhythm, either in the presence (A, C) or absence (B, D) of theta phase reset. With phase reset, both peak and trough stimulation switch network behavior from no activity to sustained oscillations. Without phase reset, only peak stimulation can induce sustained oscillations. E. Quantification of theta power, gamma power, and PAC (measured using the MI) in CA1 excitatory (blue) and inhibitory (red) populations in all four cases (metrics are computed during the pulse train, within the shaded areas of panels A-D).

Generation of gamma oscillations by E-I circuits

It is well established that a network of interconnected pyramidal neurons and interneurons can give rise to oscillations in the gamma range, a mechanism termed pyramidal-interneuronal network gamma (PING) (Traub et al., 2004; Onslow et al., 2014; Segneri et al., 2020). This mechanism has been observed in several optogenetic studies with gradually increasing light intensity (i.e. under a ramp input) affecting multiple different circuits, such as layer 2–3 pyramidal neurons of the mouse somatosensory cortex (Adesnik and Scanziani, 2010), the CA3 field of the hippocampus in rat in vitro slices (Akam et al., 2012), and in the non-human primate motor cortex (Lu et al., 2015). In all cases, gamma oscillations emerged above a certain threshold in terms of photostimulation intensity, and the frequency of these oscillations was either stable or slightly increased when increasing the intensity further. We sought to replicate these findings with our elementary E-I circuits composed of single-compartment conductance-based neurons driven by a ramping input current (Figure 2 and Figure 2—figure supplement 1). As an example, all the results in this section will be shown for an E-I circuit that has similar connectivity parameters as the CA1 field of the hippocampus in our complete model (see section ‘Hippocampal formation: inputs and connectivity’ in the Materials and methods).

For low input currents provided to both neuronal populations, only the highly excitable interneurons were activated (Figure 2A). For a sufficiently high input current (i.e. a strong input that could overcome the inhibition from the fast-spiking interneurons), the pyramidal neurons started spiking as well. As the amplitude of the input increased, the activity of both neuronal populations became synchronized in the gamma range, asymptotically reaching a frequency of about 60 Hz (Figure 2A bottom panel). Decoupling the populations led to the abolition of gamma oscillations (Figure 2B), as neuronal activity was determined solely by the intrinsic properties of each cell. Interestingly, when the ramp input was provided solely to the excitatory population, we observed that the activity of the pyramidal neurons preceded the activity of the inhibitory neurons, while still preserving the emergence of gamma oscillations (Figure 2—figure supplement 1A). As expected, decoupling the populations also abolished gamma oscillations, with the excitatory neurons spiking a frequency determined by their intrinsic properties and the inhibitory population remaining silent (Figure 2—figure supplement 1B).

To further characterize the intrinsic properties of individual inhibitory and excitatory neurons, we derived their input-frequency (I-F) curves, which represent the firing rate of individual neurons in response to a tonic input (Figure 2—figure supplement 2A). We observed that for certain input amplitudes, the firing rates of both types of neurons was within the gamma range. Interestingly, in the absence of noise, each population could generate by itself gamma oscillations that were purely driven by the input and determined by the intrinsic properties of the neurons (Figure 2—figure supplement 2B). Adding stochastic Gaussian noise in the membrane potential disrupted these artificial oscillations in decoupled populations (Figure 2—figure supplement 2C). All subsequent simulations were run with similar noise levels to prevent the emergence of artificial gamma oscillations.

Another potent way to induce gamma oscillations is to drive fast-spiking inhibitory neurons using pulsed optogenetic stimulation at gamma frequencies, a strategy that has been used both in the neocortex (Cardin et al., 2009) and hippocampal CA1 (Iaccarino et al., 2016). In particular, Cardin and colleagues systematically investigated the effect of driving either excitatory or fast-spiking inhibitory neocortical neurons at frequencies between 10 and 200 Hz (Cardin et al., 2009). They showed that fast-spiking interneurons are preferentially entrained around 40–50 Hz, while excitatory neurons respond better to lower frequencies. To verify the behavior of our model against these experimental data, we simulated pulsed optogenetic stimulation as an intracellular current provided to our reduced model of a single E-I circuit. Stimulation was applied at frequencies between 10 and 200 Hz to excitatory cells only, to inhibitory cells only, or to both at the same time (Figure 2—figure supplement 3). The population firing rates were used as a proxy for the local field potentials (LFP), and we computed the relative power in a 10 Hz band centered around the stimulation frequency, similarly to the method proposed in Cardin et al., 2009. When presented with continuous stimulation across a range of frequencies in the gamma range, interneurons showed the greatest degree of gamma power modulation (Figure 2—figure supplement 3). Furthermore, when the stimulation was delivered to the excitatory population, the relative power around the stimulation frequency dropped significantly in frequencies above 10 Hz, similar to the reported experimental data (Cardin et al., 2009). The main difference between our simulation results and these experimental data is the specific frequencies at which fast-spiking interneurons showed resonance, which was around 40 Hz in the mouse barrel cortex and around 90 Hz in our model, a fast gamma rhythm. This could be attributed to several factors, such as differences in the cellular properties between cortical and hippocampal fast-spiking interneurons, or the differences between the size of the populations and their relevant connectivity in the cortex and the hippocampus.

Theta-gamma oscillations and theta phase reset under dynamical theta input

Once we validated that our elementary E-I circuits were able to generate gamma oscillations under an input of sufficiently high amplitude, we studied the behavior of the whole model when the theta input was dynamically provided by the Kuramoto oscillators as described above (Figure 1). As in the original model (Aussel et al., 2018), the input theta rhythm drove the network to produce spiking activity in the gamma range preferentially around the peak of theta in all excitatory and inhibitory populations (Figure 3A and Figure 3—figure supplement 1). Spectrograms of the CA1 population firing rates revealed that these bursts of activity around each theta peak were characterized by oscillations around 60 Hz (Figure 3A). This was confirmed by power spectral densities, which showed a clear peak at 4 Hz (corresponding to the theta drive) and increases in gamma band activity between 40 and 80 Hz (Figure 3B). To quantify PAC between theta and gamma oscillations, we used the modulation index (MI) (Tort et al., 2008; Tort et al., 2010), which has been shown to outperform other similar measures (Hülsemann et al., 2019). However, we discovered that this metric would give erroneous results in our simulated datasets due to the absence of certain frequency components. To overcome this limitation and avoid artifacts, uniform noise was added to the firing rates prior to computing the MI (Figure 3—figure supplement 2). We first visualized the quantification of PAC using the comodulogram, which indicates the MI as a function of two frequencies corresponding to the modulating phase signal and the modulated amplitude signal. This analysis confirmed that gamma-band signals between 45 and 75 Hz were modulated by lower theta frequencies between 3 and 6 Hz, which we quantified by computing a global MI in frequency ranges encompassing theta (3–9 Hz) and gamma (40–80 Hz) (Figure 3—figure supplement 2B). To refine this analysis, we also computed a similar MI between the amplitude of various gamma frequencies and the phase of theta, which indicated that oscillations between 40 and 80 Hz occur preferentially around the peak of theta (Figure 3C and Figure 3—figure supplement 2B).

Next, we verified that our model was able to display theta phase reset during single-pulse stimulation (depolarizing pulse, 1 ms duration). This mechanism is tightly linked with the concept of PRC, which characterizes the phase delay or advancement that follows a single pulse delivered to an oscillatory system, as a function of the phase at which this input is delivered. Although there is no direct measurement of the PRC of septal neurons, such characterizations have been performed for individual pyramidal cells in the CA3 and CA1 fields of the hippocampus (Lengyel et al., 2005; Kwag and Paulsen, 2009; Akam et al., 2012). These PRCs appear biphasic and show a phase advancement (respectively delay) for stimuli delivered in the ascending (respectively descending) slope of theta. We modeled this behavior by a specific term (which we called the phase response function) in the general equation of the Kuramoto oscillators (see methods, Equation 1). Importantly, introducing a phase offset in the phase response function disrupted theta-nested gamma oscillations (Figure 3—figure supplement 3), which suggests that the septohippocampal circuitry must be critically tuned to be able to generate such oscillations. The strength of phase reset could also be adjusted by a gain that was manually tuned. In the presence of the physiological phase response function and of a sufficiently high reset gain, a single stimulation pulse delivered to all excitatory and inhibitory CA1 neurons could reset the phase of theta to a value close to its peaks (Figure 3A). We computed the PRC of our simulated data for different stimulation amplitudes and validated that our neuronal network behaved according to the phase response function set in our Kuramoto oscillators (Figure 3D). It should be noted that including this phase reset mechanism affected the generated theta rhythm even in the absence of stimulation, extending the duration of the theta peak and thereby slowing down the frequency of the generated theta rhythm.

Importantly, our approach is generalizable and can be applied to other models producing theta-nested gamma oscillations. For instance, we adapted the neural mass model by Onslow and colleagues (Onslow et al., 2014), replaced the fixed theta input by a set of Kuramoto oscillators, and demonstrated that it could also generate theta phase reset in response to single-pulse stimulation (Figure 3—figure supplement 4). These results illustrate that the general behavior of our model is not specific to the tuning of individual parameters in the conductance-based neurons, but follows general rules that are captured by the level of abstraction of the Kuramoto formalism.

Overall, we successfully developed a new model of the hippocampal formation able to exhibit both theta-nested gamma oscillations and theta phase reset in response to stimulation. We then decided to explore further the effects of various stimulation protocols on its dynamics.

Effects of theta input and stimulation amplitudes on theta-gamma oscillations

We investigated the behavior of the model across multiple states of varying septal theta input amplitude in response to single-pulse stimulation delivered to CA1. In the absence of theta drive, a single stimulation pulse elicited either zero or two bursts of spiking activity (depending on stimulation amplitude) separated by about 200 ms (Figure 4A and C). We first sought to understand the origin of this second burst, as it showed that even a single pulse could induce transient periodic activity around 5 Hz in our model, a frequency within the theta range. A previous model has shown that the presence of Calcium-Activated Non-specific cationic (CAN) currents can lead to self-sustained theta oscillations in the hippocampus (Giovannini et al., 2017). Moreover, this study showed a direct link between the increased excitation provided by the CAN current and the spike-frequency adaptation properties of the M current, directly affecting transitions from an asynchronous low-firing regime to synchronous bursting. We tested the role of the CAN current in the response to single-pulse stimulation by completely removing it from our simulations, which abolished the second burst of activity (Figure 4B). Moreover, the time interval separating the two bursts likely resulted from the interplay between the depolarizing CAN current and hyperpolarizing M current (Figure 4A).

In the presence of dynamic theta input, the effects of single-pulse stimulation depended both on theta input amplitude and stimulation amplitude, highlighting different regimes of network activity (Figure 5 and Figure 5—figure supplement 1, Figure 5—figure supplement 2, Figure 5—figure supplement 3). For low theta input, theta-nested gamma oscillations were initially absent and could not be induced by stimulation (Figure 5A). At most, the stimulation could only elicit a few bursts of spiking activity that faded away after approximately 250ms, similar to the rebound of activity seen in the absence of theta drive. For increasing theta input, the network switched to an intermediate regime: upon initialization at a state with no spiking activity, it could be kicked to a state with self-sustained theta-nested gamma oscillations by a single stimulation pulse of sufficiently high amplitude (Figure 5B). This regime existed for a range of septal theta inputs located just below the threshold to induce self-sustained theta-gamma oscillations without additional stimulation, as characterized by the post-stimulation theta power, gamma power, and theta-gamma PAC (Figure 5D). Removing CAN currents from all areas of the model abolished this behavior (Figure 5 - figure supplement 4), which is interesting given the role of this current in the multistability of EC neurons (Egorov et al., 2002; Fransén et al., 2006) and in the intrinsic ability of the hippocampus to generate theta-nested gamma oscillations (Giovannini et al., 2017). For the highest theta input, the network became able to spontaneously generate theta-nested gamma oscillations, even when initialized at a state with no spiking activity and without additional neurostimulation Figure 5C.

Neurostimulation for restoring theta-gamma oscillations in pathological states

Based on the above analyses, we considered two pathological states: one with a moderate theta input (i.e. moderately weak projections from the medial septum to the EC) that allowed the initiation of self-sustained oscillations by single stimulation pulses (Figure 5, point B), and one with a weaker theta input characterized by the complete absence of self-sustained oscillations even following transient stimulation (Figure 5, point A). In each case, we sought to assess whether single-pulse or pulse train stimulation could induce or restore theta-nested gamma oscillations and whether this effect depended on the phase at which stimulation was delivered (i.e. at the peak or trough of the theta cycle). We hypothesized that any possible phase relationship would also depend on the phase reset mechanism. To test this hypothesis, we ran a series of simulations using two different models: one without phase reset and one with strong phase reset (i.e. the reset gain was set at the value used in Figure 3).

In the case of a moderate theta input and in the presence of phase reset, delivering a pulse at either the peak or trough of theta could induce theta-nested gamma oscillations (Figure 6A and C). By contrast, in the absence of phase reset, only stimulation delivered at the peak of theta was able to induce such oscillations, after some time delay (Figure 6B and D). Quantification of these results in terms of theta power, gamma power, and theta-gamma PAC showed strongly similar responses between stimulation delivered at the peak or trough of theta with phase reset enabled, similar but weaker responses with stimulation delivered at the peak of theta with phase reset disabled, and no response in the case of trough stimulation with phase reset disabled (Figure 6E).

In the case of a weak theta input that completely abolished neuronal oscillations, we delivered stimulation pulses continuously at a frequency matching that of the underlying theta rhythm for a duration of 2 s, in order to restore physiological oscillations (Figure 7). The stimulation onset was timed to either the peak or trough of the ongoing theta cycle. The continuous delivery of pulses produced similar results as single-pulse stimulation. With phase reset, pulse train stimulation restored theta-nested gamma oscillations within the whole network, irrespective of the phase of stimulation onset (Figure 7A and C). Notably, stimulation delivered at the trough forced a reset of the phase of theta rhythm and led to the subsequent delivery of pulses at the peak of theta. Interestingly, in the absence of phase reset, peak-targeted stimulation induced theta-gamma oscillations in all network areas easily, while trough-targeted stimulation created artificial bursts that propagated in other areas with difficulty, requiring multiple pulses to achieve a fraction of the results of peak-targeted stimulation (Figure 7B and D). Comparing these simulations based on theta power, gamma power, and theta-gamma PAC within CA1 (Figure 7E) showed similar albeit less striking differences as single-pulse stimulation, possibly because CA1 was driven directly by the stimulation. Notably, gamma power in the absence of theta phase reset, was higher when utilizing pulse trains. These differences were even more pronounced in areas other than CA1, as can be seen by the gradual emergence of oscillations for trough stimulation without phase reset (Figure 7D).

Discussion

Highlights

In summary, we have developed a novel computational model to investigate the effects of electrical stimulation on a slice of the hippocampal formation, incorporating two important features related to memory: theta-nested gamma oscillations and theta phase reset. The key innovation compared to previous models (e.g. Aussel et al., 2018) is the introduction of a set of abstract Kuramoto oscillators, which represent pacemaker neurons in the medial septum and are interfaced with biophysically realistic neuronal models in the hippocampus. From a methodological point of view, this hybrid interfacing between two levels of abstraction represents an innovation in itself and could be applied to other systems or brain structures that are driven by dynamical rhythms. The main outcomes reported here relate to the importance of the theta reset mechanism when examining the effects of neurostimulation on hippocampal oscillations.

A very interesting finding concerns the behavior of the model in response to single-pulse stimulation for certain values of the theta amplitude (Figure 5). For low theta amplitudes, a single stimulation pulse was capable of switching the network behavior from a state with no spiking activity to one with prominent theta-nested gamma oscillations. Whether such an effect can be induced in vivo in the context of memory processes remains an open question. Nevertheless, delivering a single stimulation pulse bilaterally to the human hippocampus during a memory task is sufficient to impair memory encoding (Lacruz et al., 2010), suggesting that even single-pulse stimulation can indeed have wide network effects that are behaviorally relevant.

The second main finding is that the timing of individual stimulation pulses with respect to the phase of the ongoing theta rhythm matters differently depending on the presence or absence of phase reset (Figure 6 and Figure 7). Human intracranial stimulation data indicate that the receptivity of hippocampal circuits to single-pulse stimulation is modulated by the phase of theta (Lurie et al., 2022). A number of studies have also reported conflicting results in terms of memory outcomes, which could potentially be attributed to the induction of phase reset through stimulation (Suthana et al., 2012; Jacobs et al., 2016). To our knowledge, however, the degree of phase reset that follows each stimulation pulse remains unknown and should be investigated in future experimental studies. From a technological point of view, the two regimes with and without phase reset have opposite predictions concerning the need for closed-loop stimulation protocols that would trigger stimulation in real-time based on the phase of ongoing theta oscillations. Such phase-triggered stimulation would be most useful if the phase reset mechanism remains relatively limited. When this mechanism becomes too strong, regular continuous stimulation appears sufficient to restore physiological theta-gamma oscillations (Figure 7).

It should also be noted that the reset gain in our simulations (Equation 1, gain Greset) was either completely turned off or set at a high value, producing a strong effect that always reset the phase of theta to its peak. In reality, the degree of theta phase reset is dynamic, depending on the environment and associated task requirements (Rizzuto et al., 2003; Mormann et al., 2005; Jackson et al., 2008), and may be affected by neurodegenerative disorders that affect the connections between the hippocampal formation and the medial septum. Importantly, our proposed framework can simulate intermediate values of the reset gain, which should ideally be fitted to experimental data in future applications.

Finally, we modeled pathological states by reducing the maximum amplitude of the theta input (Equation 2, gain Gθ) until theta-nested gamma oscillations were impaired or even abolished (Figure 5). This choice was meant to simulate neurodegeneration in the medial septum, which is known to be affected in Alzheimer’s disease, leading to oscillatory disruptions (Nelson et al., 2014; Hampel et al., 2018; Takeuchi et al., 2021). Another possibility would be that neurodegeneration limits the ability of the septal pacemaker neurons to synchronize, thus producing a weaker collective theta rhythm without affecting the maximum amplitude of individual oscillators. Although we simulated this change in our model by reducing the synchronization parameter, the effects on hippocampal oscillations were less pronounced (Figure 5 - figure supplement 5). Linking these different modeling parameters to experimental biomarkers will be important in future work.

Limitations

Even though we took great care in developing a precise representation of the hippocampal formation, the resulting model remains a simplification that could be further enriched. In particular, we deliberately modeled only a single theta generator, while multiple intra- and extra-hippocampal generators are known to co-exist (Kocsis et al., 1999; Hummos and Nair, 2017). We decided to model septal pacemaker neurons projecting to the EC as the main source of hippocampal theta as reported in multiple experimental studies (Buzsáki, 2002; Buzsáki et al., 2003; Hangya et al., 2009; Colgin, 2013). However, experimental findings and previous models have also proposed that direct septal inputs are not essential for theta generation (Wang, 2002; Colgin, 2013; Mysin et al., 2019), but play an important role in phase synchronization of hippocampal neurons. Furthermore, the model does not account for the connections between the lateral and medial septum and the hippocampus (Takeuchi et al., 2021). These connections include the inhibitory projections from the lateral to the medial septum and the monosynaptic projections from the hippocampal CA3 field to the lateral septum. An experimental study has highlighted the importance of the lateral septum in regulating the hippocampal theta rhythm Bender et al., 2015, an area that has not been included in the model. Specifically, theta-rhythmic optogenetic stimulation of the axonal projections from the lateral septum to the hippocampus was shown to entrain theta oscillations and lead to behavioral changes during exploration in transgenic mice. To account for these discrepancies, our model could be extended by considering more realistic connectivity patterns between the medial/lateral septum and the hippocampal formation, including glutamatergic, cholinergic, and GABAergic reciprocal connections (Müller and Remy, 2018), or by considering multiple sets of oscillators each representing one theta generator.

In terms of neuronal cell types, we also made an important simplification by considering only basket cells as the main class of inhibitory interneuron in the whole hippocampal formation. However, it should be noted that many other types of interneurons exist in the hippocampus and have been modeled in various works with higher computational complexity (e.g. Bezaire et al., 2016; Chatzikalymniou et al., 2021). Among these various interneurons, oriens-lacunosum moleculare (OLM) neurons in the CA1 field have been shown to play a crucial role in synchronizing the activity of pyramidal neurons at gamma frequencies (Tort et al., 2007), and in generating theta-gamma PAC (e.g. Neymotin et al., 2011; Ponzi et al., 2023). Additionally, these cells may contribute to the formation of specific phase relationships within CA1 neuronal populations, through the integration between inputs from the medial septum, the EC, and CA3 (Mysin et al., 2019). Future work is needed to include more diverse cell types and detailed morphologies modeled through multiple compartments.

Another limitation of our model concerns synaptic transmission delays, which have been largely neglected and could affect the phase relationships between the medial septum and different hippocampal subfields. Experimental studies have indeed reported time delays in the population activities of connected anatomical structures (i.e. from EC to DG; Mizuseki et al., 2009), with pyramidal cells in downstream areas like CA3 and CA1 preferentially firing at different phases of theta (Dragoi and Buzsáki, 2006). Propagation effects could also depend on the spatial scale of the model. We also decided to represent only a thin coronal slice of the hippocampal formation, and it remains unclear how an anatomically accurate model of the whole structure would behave in terms of propagation of spontaneous and electrically-induced neuronal activity.

Importantly, we did not consider learning through synaptic plasticity, even though such mechanisms could drastically modify synaptic conduction for the whole network (Borges et al., 2017). Even more interestingly, the inclusion of spike-timing-dependent plasticity would enable the investigation of stimulation protocols aimed at promoting LTP, such as theta-burst stimulation (Larson and Munkácsy, 2015). This aspect would be of uttermost importance to make a link with memory encoding and retrieval processes (Axmacher et al., 2006; Tsanov and Manahan-Vaughan, 2009; Jutras et al., 2013) and with neurostimulation studies for memory improvement (Titiz et al., 2017; Solomon et al., 2021).

From the point of view of neurostimulation, future work is needed to extend the current model to a multi-compartment representation of neurites, since axons are known to be preferentially activated by extracellular electrical stimulation (Rattay et al., 2003). Specifically, multi-compartment cable models have been developed to investigate spike initiation and propagation by modeling the axon as a series of resistance-capacitance circuits following its trajectory (Rattay et al., 2003; Joucla and Yvert, 2012; Ashida and Nogueira, 2018). These models are particularly suited to study the effects of extracellular stimulation, which are non-intuitive and depend as a first approximation on the second spatial derivative of the electrical potential along the cell membrane (Rattay, 1986; Rattay et al., 2003; McIntyre et al., 2004; Rattay et al., 2018). Here, we have modeled electrical stimulation as an intracellular current applied equally across all neurons in the targeted area, which is extremely simplified but enables computational tractability. Future developments will focus on developing an equivalent multicompartment model with realistic axonal trajectories while making sufficient simplifications to allow for realistic computation times.

Finally, we likened conditions of low theta input to pathological states characteristic of oscillopathies such as Alzheimer’s disease, as these conditions disrupted all aspects of theta-gamma oscillations in our model: theta power, gamma power, and theta-gamma PAC (Figure 5). However, it should be noted that changes in theta or gamma power in these pathologies are often unclear, and that the most consistent alteration that has been reported in Alzheimer’s disease is a reduction of theta-gamma PAC (for review, see Kitchigina, 2018). Future work should explore the effects of cellular alterations intrinsic to the hippocampal formation and their impact on theta-gamma oscillations.

Outlook

Overall, this new model of the hippocampal formation represents a methodologically innovative basis to further explore multiple neurostimulation strategies that target hippocampal oscillations. Moreover, the limitations discussed above represent important avenues for future refinements, which will require significant work to overcome the costs in terms of computational tractability (e.g. modeling the whole hippocampus, or using multicompartment models of axonal trajectories and dendritic trees). Ultimately, such model refinements should allow the investigation of the effects of extracellular stimulation using charge-balanced biphasic pulses and multipolar electrode configurations with realistic electrode geometries. Most importantly, we believe that our current work may also serve as an inspiration for future computational models of oscillopathies (not necessarily limited to the hippocampus), which could benefit from interfacing abstract sets of synchronizing oscillators and investigating their interactions with biophysically-realistic neurons.

Materials and methods

Computational model

Overall architecture

We aimed to develop a computational model of the hippocampal formation that is able to generate theta-nested gamma oscillations and takes into account the dynamic nature of the theta input from the medial septum, that is the reset of the theta phase following strong activity in the perforant path or fornix (Buño et al., 1978; Williams and Givens, 2003). To this end, we adapted a previous biophysical model of the human hippocampal formation under fixed sinusoidal input (Aussel et al., 2018) and interconnected it with an abstract representation of the medial septum, modeled as an assembly of Kuramoto oscillators (Kuramoto, 1984; Figure 1).

More precisely, we modeled a 15-mm-thick coronal slice of the hippocampal formation which comprises the DG, the CA3 and CA1 subfields of the hippocampus, and the EC, all populated with excitatory and inhibitory Hodgkin-Huxley neurons. For the medial septum, we opted for a simplified representation where the medial septum projects only to the EC, which in turn drives the different subfields of the hippocampus (see also ‘A computational model of the hippocampal formation with dynamical theta input’, second paragraph). In our model, Kuramoto oscillators are therefore connected to the EC neurons and they receive projections from CA1 neurons (see sections below for more details).

Medial septum: Kuramoto oscillators

The medial septum contains pacemaker neurons (Varga et al., 2008; Hangya et al., 2009) that synchronize with one another and generate a global theta rhythm. In turn, these pacemaker neurons drive downstream areas that receive projections from the medial septum. Here, we modeled the entire medial septum neuronal assembly as a set of coupled phase oscillators following the Kuramoto model (Kuramoto, 1984), which generated the driving theta rhythm that was provided to the hippocampal formation through the EC (Breakspear et al., 2010).

Kuramoto oscillators have also been used to investigate the effects of neurostimulation on synchronized brain rhythms in the context of Parkinson’s disease or essential tremor (Tass, 2003; Weerasinghe et al., 2019). Here, each Kuramoto oscillator was described by its phase θi, which evolves over time according to the following equation:(1) dθidt=ωi+kN∑j=1Nsin⁡(θj−θi)+GresetX(t)Z(θi)

The term ωi denotes the natural frequency of oscillator i, and is normally distributed around the center frequency f0 and with standard deviation σ. The synchronization parameter k represents the coupling strength of the group of oscillators, and N is the number of oscillators. Higher values for the synchronization parameter indicate stronger coupling between pairs of oscillators, thus affecting how fast or slow their phases tend to synchronize. The final product GresetX(t)Z(θi) is used to describe the effects of external inputs on the phase of the oscillators. Here, the only external input originates in the projections from CA1 to the medial septum. It is described by the instantaneous firing rate X(t) of the CA1 excitatory population. Greset is an arbitrary gain that determines the strength of the CA1 input to the oscillators. Finally, the function Z(θi) describes how the effects of CA1 inputs on the oscillators depend on the phase of the ongoing theta rhythm.

The coherence and phase of the driving theta rhythm Iθ were computed using the order parameter r to extract the mean amplitude A and the mean phase ϕ of the ensemble of Kuramoto oscillators as follows:(2) {r(t)=1N∑i=1Nejθi(t)A(t)=Re⁡(r)ϕ(t)=Im⁡(r)Iθ(t)=GθA(t)cos⁡(ϕ(t))+12

where the mean amplitude A(t) and the mean phase ϕ(t) are derived by taking the real and imaginary part of the order parameter respectively, and the output theta rhythm Iθ(t) is a rectified cosine, multiplied by a gain Gθ that controls the maximum output amplitude (in nA).

To simulate the effect of the projections from CA1 to the medial septum on the phase of the oscillators, an approximation of the instantaneous firing rate of the CA1 excitatory population was used as the term X(t). To obtain the instantaneous firing rate, we convolved the CA1 population spike train with an exponential kernel using the following equations:(3) {dXdt=−XτFRXt+1=Xt+1NτFR,∀t∈tS

where tS denotes the ordered set of the spike timings, and τFR determines the exponential decay of the kernel, which was set to 10ms to compute population activity (for single neurons, a typical value is 100ms) (Gerstner et al., 2014).

Hereafter, we call the term Z(θ) the phase response function, to distinguish it from the PRC obtained from experimental data or simulations (see section below ‘Data Analysis’, ‘Phase Response Curve’). Briefly, the PRC of an oscillatory system indicates the phase delay or advancement that follows a single pulse, as a function of the phase at which this input is delivered. The phase response function Z(θ) was chosen to mimic as well as possible experimental PRCs reported in the literature (Lengyel et al., 2005; Kwag and Paulsen, 2009; Akam et al., 2012). These PRCs appear biphasic and show a phase advancement (respectively delay) for stimuli delivered in the ascending (respectively descending) slope of theta. To accurately model this behavior, we used the following equation for the phase response function, where θpeak represents the phase at which the theta rhythm reaches its maximum and the parameter ϕoffset controls the desired phase offset from the peak:(4) Z(θi)=−sin⁡(θi−(θpeak+ϕoffset))

An overview of the default parameters for the Kuramoto oscillators and the bidirectional connections between medial septum-hippocampal formation can be found in Table 1.

Table 1. Full list of the default parameter values for the Kuramoto oscillators.

Parameter	Value	
Number of oscillators (N)	250	
Center frequency (f0)	6 Hz	
Standard deviation (σ)	0.5 HZ	
Synchronization ratio (k/N)	15	
Phase reset gain (Greset)	4	
Peak phase (θpeak)	0 rad	
Firing rate time constant (τFR)	10 ms	

Hippocampal formation: Hodgkin-Huxley Neurons

The following sections describe in detail how individual neurons and synapses were modeled, and are adapted from the original work by Aussel and colleagues (Aussel et al., 2018). Neurons were modeled as conductance-based single compartments, following the Hodgkin-Huxley formalism (Hodgkin and Huxley, 1952), in line with previous work (Aussel et al., 2022; Aussel et al., 2018). The temporal evolution of the membrane potential of each neuron is described by a differential equation whose general form reads:(5) CmdVmdt=−IL−∑channelIchannel−∑j∈[E,I]Isynj+Iθ+Istim+η

IL denotes the leakage current. Ichannel are currents associated with specific ion channels, namely potassium (IK), fast sodium (INa), and low-threshold calcium (ICa) currents, the CAN current (ICAN) (Giovannini et al., 2017), and the M-type potassium channel current (IM) responsible for spike adaptation (Kosenko et al., 2012; Sun and Kapur, 2012; Kwag et al., 2014). Isyn represents the currents originating from synaptic inputs to the cell and can be either depolarizing (negative sign) or hyperpolarizing (positive sign). η is a Gaussian random noise term accounting for other external inputs and synaptic fluctuations. Theta input from the medial septum is modeled as a depolarizing current and is denoted by Iθ, while electrical stimulation is denoted by Istim. Excitatory neurons represent pyramidal cells in EC, CA3, and CA1, and granule cells in DG. They were modeled with INa, IK, and ICa, and IM currents, with the addition of ICAN for pyramidal cells. Fast-spiking interneurons in all areas were modeled with INa and IK currents. The complete description for all of the above ionic channels and their corresponding currents can be found in (Giovannini et al., 2017). Leakage currents followed the following equation:(6) IL=(gL×A)×(Vm−EL)

where gL is the maximum leakage conductance, A is the area of the single compartment corresponding to the membrane of a neuron, and EL is the reversal potential of the leakage channel.

Channel currents IK, IM, ICAN obey the following set of equations:(7) IK=gK×A×n4×(Vm−EK)IM=gM×A×p×(Vm−EM)ICAN=gCAN×A×mCAN2×(Vm−ECAN)

where gK, gM, gCAN are the maximum conductances for the respective channel, and n, p, mCAN are the respective gating variables defined by the following differential equations:(8) dndt=n∞−nτndpdt=p∞−pτpdmCANdt=mCAN,∞−mCANτmCAN

For the potassium and CAN currents, the steady-state values for their corresponding gating variables n∞ and mCAN,∞ and their corresponding time constants τK and τCAN depend on the following functions of the transition rate constants:(9) n∞=αnαn+βnmCAN,∞=αmCANαmCAN+βmCANτn=0.2αn+βnτmCAN=0.2αmCAN+βmCAN

The sodium current (INa) and calcium current (ICa) follow a similar set of equations:(10) INa=gNa×A×m3×h×(Vm−ENa)ICa=gCa×A×m2×h×(Vm−ECa)

with two gating variables m and h defined by the following differential equations:(11) dmdt=nα−nτndndt=nα−nτnm∞=αmαm+βmh∞=αhαh+βhτm=0.2αm+βmτh=0.2αh+βh

The gating variable of ICAN depends on the calcium concentration within the neuron ([Ca]i2+), given by:(12) d[Ca]i2+dt=γ(ICa)+([Ca]∞2+−[Ca]i2+)τ[Ca]2+γ(ICa)=−ku×ICa2×F×d×A

where τ[Ca]2+ = 1s represents the rate of calcium removal from the cell, [Ca]∞2+ = 0.24 mol/L is the calcium concentration if the calcium channel remains open for a duration of ΔT→∞, ku=104 is a unit conversion constant, F is the Faraday constant, and d = 1 μm is the depth at which the calcium is stored inside the cell.

Noise (η), accounting for random inputs to the network, was simulated as intracellular current acting on the membrane voltage and following the properties of a Gaussian random variable with a mean of 0 μV and a standard deviation of 1000 μV (ηE∼N(0,1000) μV) for excitatory neurons and a mean of 0 μV and standard deviation of 100 μV (ηI∼N(0,100) μV) for inhibitory neurons. The ratio of 1:10 between the noise terms was adapted from the original work and it accounts for the higher excitability of the inhibitory neurons as well as the E-I population size ratios.

The original model introduced some parameters representing the vigilance state (i.e. active wakefulness vs slow-wave sleep). However, the present model only focused on the state of active wakefulness, since this is when memory-related theta-nested gamma oscillations occur. For all simulations, the parameters were set so that the network operated in the wakefulness regime in a healthy hippocampus (Aussel et al., 2018). The full expressions for all the parameters defined above can be found in Table 2 for pyramidal cells and in Table 3 for interneurons.

Table 2. Full list of parameter values and expressions for pyramidal neurons.

Parameter	Expression	
A	29.103μm2	
Cm	1μF/cm2	
gL	0.01mS/cm2	
EL	−70mV	
gK	5mS/cm2	
EK	−100mV	
αn,K	−0.032Vm+40mV−1+e−0.02(Vm+40mV)
	
βn,K	0.5e−(Vm+45mV)40mV
	
gNa	50mS/cm2	
ENa	50mV	
αm,Na	−0.32Vm+42mVe−Vm+42mV4mV−1
	
βm,Na	0.28Vm+15mVe−Vm+15mV5mV−1
	
αh,Na	0.128e−Vm+38mV18mV
	
βh,Na	41+e−Vm+15mV5mV
	
gM	90μS/cm2	
EM	−100mV	
p∞	11+e−0.01(Vm+35mV)
	
τp	13.3e(Vm+35mV)20mV+e−(Vm+35mV)20mV
	
gCa	0.1mS/cm2	
ECa	120mV	
αm,Ca	−0.055Vm+27mVe−Vm+27mV17mV−1
	
βm,Ca	−0.94eVm+75mV17mV
	
αh,Ca	−0.000457eVm+13mV50mV
	
βh,Ca	0.0065e−Vm+15mV28mV
	
gCAN	25μS/cm2	
ECAN	−20mV	
αm,CAN	0.0002e1.4[Ca]i2+0.5mol/L	
βm,CAN	0.0002e1.4	

Table 3. Full list of parameter values and expressions for interneurons.

Parameter	Expression	
A	14.103μm2	
Cm	1μF/cm2	
gL	0.1mS/cm2	
EL	−65mV	
gK	9mS/cm2	
EK	−90mV	
αn,K	0.01Vm+34mV1−e−0.1(Vm+34mV)
	
βn,K	0.125e−Vm+44mV80mV
	
gNa	35mS/cm2	
ENa	55mV	
αm,Na	0.1Vm+35mV1−e−0.1(Vm+35mV)
	
βm,Na	4e−Vm+60mV18mV	
αh,Na	0.07e−Vm+58mV20mV	
βh,Na	11+e−0.1(Vm+28mV)
	

Synaptic models

Inter-neuronal interactions were modeled as instantaneous AMPA and GABA-A synapses using the synaptic currents IsynE and IsynI, respectively. Synaptic currents were described by the following bi-exponential differential equations:(13) IsynI,E=gI,E(Vm−EI,E)dgI,Edt=1τgI,E(−gI,E+hI,E)dhI,Edt=−hI,E1τhI,E

where EI,E are the synaptic resting potentials, and τgI,E and τhI,E are the synaptic time constants of rise and decay for inhibitory and excitatory neurons respectively. The occurrence of a pre-synaptic spike leads to an increase of the values hI or hE in the post-synaptic neuron by a fixed amount, which depends on the type of synapse and the region (due to the presence of cholinergic effects described in the initial model). Specific values for the intra-area and inter-area synaptic connections are given in Table 4 and Table 5, respectively.

Table 4. A pre-synaptic spike causes an increase in the conductances he and hi in the post-synaptic neuron.

The values for the intra-area connections are given here. Empty cells indicate no connection between the populations.

	E→E	E→I	I→E	I→I	
EC		20 pS	600 pS		
DG		180 pS	1800 pS		
CA3	20 pS	20 pS	600 pS		
CA1		60 pS	1800 pS	1800 pS	

Table 5. A pre-synaptic spike causes an increase in the conductances he and hi in the post-synaptic neuron.

The values for the inter-area connections are given here. Empty cells indicate no connection between the areas. Recurrent projections are not allowed and are marked with dashes.

Source	Target	
EC	DG	CA3	CA1	
EC	-	20 pS	20 pS	20 pS	
DG		-	180 pS		
CA3			-	20 pS	
CA1	60 pS			-	

Hippocampal formation: neuron types and numbers

Each area of the network is comprised of two populations, one excitatory and one inhibitory. Excitatory cells in the DG represent granule cells and pyramidal neurons in all other areas. Interneurons represent basket cells across all areas. The ratio between pyramidal neurons and interneurons was directly adapted from Aussel et al., 2018. The ratio between pyramidal neurons and interneurons was kept as a ratio of 10:1 for all areas except the dentate gyrus, where the ratio was 100:1. The number of neurons per subfield of the hippocampal formation is summarized in Table 6.

Table 6. Number of neurons per subfield of the hippocampal formation, divided by neuron type.

Area	NExc	NInh	
EC	10,000	1,000	
DG	10,000	100	
CA3	1,000	100	
CA1	10,000	1,000	

A two-dimensional simplified image depicting a coronal slice of the hippocampal formation (Aussel et al., 2018) was used as a basis for a two-dimensional manifold that was uniformly populated by neurons following a density-driven approach (Rougier, 2018). Pyramidal neurons were uniformly distributed within the stratum pyramidale (or within the stratum granulosum for the dentate gyrus) and interneurons were uniformly distributed within the stratum oriens. Initial neuron positions were drawn from a blue noise distribution and a Voronoi diagram was computed. To adjust the positions of the neurons over a centroidal Voronoi diagram, the Lloyd relaxation algorithm was applied for 1000 iterations. Transitioning from a two-dimensional manifold to a 3D reconstruction of the hippocampal formation was achieved through the addition of the third coordinate with values uniformly distributed between 0 and 15 mm.

Hippocampal formation: inputs and connectivity

The hippocampal formation receives most of its external inputs from the EC. Activity from the EC is projected to all hippocampal subfields, starting with the DG. DG granule cells project onto the CA3 pyramidal cells via mossy cell fibers. CA3 projects in turn to CA1 via Schaffer collaterals. These connections form the tri-synaptic pathway. Direct connections from the EC towards the CA3 and CA1 subfields through the monosynaptic pathway are also considered. Pyramidal neurons from CA1 project to pyramidal neurons and interneurons in the EC, closing the hippocampal-entorhinal loop. CA1 pyramidal neurons also project to the medial septum through the fornix. An overview of the connections between areas is presented in panel B of Figure 1.

In the current model, EC pyramidal neurons and interneurons receive oscillatory theta input from the medial septum in the form of an excitatory intracellular current as described in Equation 2 and Equation 5. Projections from CA1 towards the medial septum were modeled as a signal representing the collective firing rate of the CA1 pyramidal neurons. All connections towards and from the medial septum are summarized in panels A and B of Figure 1.

Synaptic connectivity between neurons within a region is characterized using a probability p and is distance-based, following a Gaussian-like distribution (Equation 14) with a width σ of 2500 μm (excitatory synapses) and 350 μm (inhibitory synapses). The value of Aintra defines the maximum probability of connection between two neurons separated by an infinitesimal distance (i.e. for D→0). The maximum values for Aintra are given in Table 7.(14) p=Aintraexp⁡(−D22σ2)

Table 7. Maximum probability of connection (Aintra, Equation 14) between neurons within each region.

	Py-Py	Py-Inh	Inh-Py	Inh-Inh	
EC	0	0.37	0.54	0	
DG	0	0.06	0.14	0	
CA3	0.56	0.75	0.75	0	
CA1	0	0.28	0.3	0.7	

Inter-area connectivity followed a similar Gaussian-like distribution (Equation 15), however, only excitatory projections were considered and the distance Dz was computed only across the z-coordinate. Pyramidal neurons from the source area projected to pyramidal neurons and interneurons in the target area, with connectivity probabilities drawn from the said distribution with a width σ of 1000 μm.(15) p=min(1,Ainterexp⁡(−Dz22σ2))

Tuning input gain and connection strengths: targeted firing rates

We adjusted the input gain (Gθ) and inter-area connection strengths (Ainter) to target an overall oscillatory rhythm fosc,targ at the driving frequency of the input, that is, 6 Hz, and a mean firing rate of each excitatory population fexc,targ also at 6 Hz (meaning that each excitatory cell should spike on average once per theta cycle). Because of the ratio between excitatory and inhibitory neurons in the model, this resulted in firing rates of about 60 Hz in inhibitory neurons. These targeted values were inspired by literature in behaving rodents showing that hippocampal pyramidal neurons typically fire at rates below 10 Hz, usually between 1 and 2 Hz, and that interneurons fire at rates between 20 and 80 Hz (Hirase et al., 2001). In practice, the obtained firing rates were constrained by the simplifications made in the model.

The mean population firing rate fexc was computed by counting the number of spikes in the last second of a 3 s simulation run, to avoid edge effects due to the non-physiological initial conditions, and by averaging this number over time and across neurons. The oscillatory frequency fosc was computed as the mean of the inverse of the timing between all pairs of consecutive peaks in the theta rhythm. Finally, the metric used to adjust parameters (input strength and connection strengths) was calculated as the Euclidean distance between the targeted and obtained firing rate and oscillatory rate:(16) J=(fexc−fexc,targ)2+(fosc−fosc,targ)2

Tuning input and connection strengths: detailed procedure

The input gain (Gθ, Equation 2) and inter-area connection strengths (Ainter, Equation 15) were sequentially adjusted using the following heuristics. All simulations were performed for a duration of 3 s, and the first 2 s were excluded from the analysis to avoid edge effects. The network was initialized with membrane voltages uniformly distributed in the range [-70, -60] mV. Our initial point was a fully uncoupled model in which all the connection strengths were set to 0. The tuning procedure was performed in the absence of noise.

First, we adjusted the amplitude of the theta input to a value that would produce a mean firing rate of 6 Hz in the excitatory population of the EC. To this end, the input strength was progressively increased from 0 to 0.5 nA in steps of 0.01 nA, and the value that maximized our metric was selected.

We then adjusted similarly the connection strength from the EC to the DG, increasing it from 0 to 20 in steps of 0.1.

We repeated the procedure to jointly tune the connection strengths from EC and DG to CA3, assuming an equal contribution from these two areas. These two connection strengths were set to the same value and were increased from 0 to 5 in steps of 0.1.

We followed the same procedure to jointly tune the connection strengths from CA3 and EC to CA1, assuming again an equal contribution from these two areas.

The last connection from CA1 to EC was more difficult to tune because it closes the feedback loop and induces complex dynamics when too strong. To tune it, we temporarily uncoupled CA1 from all other structures and temporarily added a fixed sinusoidal input at the target frequency of 6 Hz to both excitatory and inhibitory CA1 neurons. We first adjusted this input amplitude to generate a mean firing rate of about 6 Hz in CA1 excitatory neurons. Next, the connection strength from CA1 to EC was adjusted to also achieve a mean firing rate of 6 Hz in EC excitatory neurons.

Because EC receives both an external input and feedback from CA1, we assumed that both these contributions should be decreased for the system to behave in a physiological range in the presence of the feedback loop. We, therefore, decided to divide by two the connection strength from CA1 to EC obtained in step 5. All other connection strengths were then set again to the values tuned in the previous steps, and the temporary external input to CA1 was removed. The external theta input to the EC was reinstated and tuned again as in step 1.

Following the above procedure, the inter-area synaptic connectivity parameters were set for all subsequent simulations. The values are summarized in Table 8. Empty cells denote no effective connectivity.

Table 8. Inter-area connection strengths (Ainter).

The source is always the excitatory population of the subfield. The same values are used when targeting excitatory and inhibitory populations. Empty cells indicate no connections. EC: Entorhinal Cortex, DG: Dentate Gyrus.

Source	Target	
EC	DG	CA3	CA1	
EC	-	13.0	0.14	1.1	
DG		-	0.14		
CA3			-	1.1	
CA1	0.2			-	

Numerical implementation

The model was implemented with the Brian2 libraries for Python (Stimberg et al., 2019; see Data availability section for access to the code). Simulations were performed using a timestep of 0.1ms and a total simulation duration ranging between 3 and 10 s (depending on the experiment, see results section). The average time for simulating 1 s for the complete model locally was 10 min.

Neural mass model interfaced with Kuramoto oscillators

The neural masses represented in Figure 3—figure supplement 4, were modeled using the Wilson-Cowan formalism, with parameters adapted from Onslow et al., 2014. Specifically, the firing rates of the excitatory and inhibitory populations were determined by the following equations:(17) {τEdEdt=−E+f(gEθE+WEEE−WIEI+stimE(t))τIdIdt=−I+f(gIθI+WEIE)

with the following parameters: τE=τI=3.2ms, gE=0.7, gI=0, WEE=4.8, WEI=WIE=4, and WII=0. The sigmoid response function was not modified from the original work, and was defined as:(18) f(x)=11+e−β(x−xm)

with parameters β and xm set to 4 and 1 respectively, according to the original model. The Kuramoto oscillators were modeled according to Equation 1. For these simulations we used a set of N=100 oscillators with a center frequency f0 of 4 Hz, a synchronization ratio kN of 25, and a strong phase reset gain Greset of 90. Other parameters were kept as previously shown in Table 1. The neural masses were coupled to the Kuramoto oscillators by linking the variable X(t) in Equation 1 with the variable E(t) in Equation 17.

We increased the phase reset gain significantly compared to the Hodgkin-Huxley model, as the Onslow model utilized a sigmoid function with values between 0 and 1, whereas the instantaneous firing rate of the populations of single-compartment neurons was much higher and therefore had a stronger phase resetting effect.

Data analysis

During each simulation, we monitored and exported the following data for subsequent analysis: (i) spike timings per neuron, (ii) time series of ionic currents, (iii) septal input theta rhythm and phase, and (iv) time series of electrical stimulation. Theta phase was wrapped between [−π,π] with a phase of 0 radians corresponding to the peak of theta rhythm. All PAC analyses were performed using the TensorPAC toolbox for Python (Combrisson et al., 2020).

Firing rates

To obtain instantaneous firing rates, the corresponding spike trains were binned in 5 ms rectangular windows with a 90% overlap (i.e. consecutive windows were spaced by 0.5ms). The number of spikes within each bin was normalized by the bin size and the number of neurons in the group, yielding instantaneous population firing rates. Where reported (i.e. values μI and μE in Figure 3A), the mean population firing rates within a given time window (typically lasting several seconds) were computed by binning all spikes in that time window and normalizing by the window width and the number of neurons in the population.

Spectral analyses

Power spectral density (PSD) estimates were calculated based on Welch’s method. The periodogram was computed using 1 s windows with a 90% overlap, yielding a frequency resolution of 1 Hz. The average spectral power within a specific frequency band was calculated using Simpson’s rule within the desired frequency band μI. Spectrograms were computed using the short-time Fourier transform with a sliding Hann window of 100 ms width and 99% overlap, yielding a frequency resolution of 10 Hz.

Modulation index

The MI (Tort et al., 2008) was used to estimate the degree of PAC between theta and gamma oscillations in the model. To compute the MI, the normalized firing rate traces were band-pass filtered in the frequency ranges of interest: 3–9 Hz for theta (referred to as the ‘phase signal’) and 40–80 Hz for gamma (referred to as the ‘amplitude signal’). Then, the phase and amplitude time series were extracted from the filtered signals using the Hilbert transform. A histogram of the mean amplitude of gamma over the phase of theta was then extracted, using phase bins of 5 degrees. The MI was finally calculated as the Kullback-Leibler divergence between the mean amplitude distribution and the uniform distribution. A higher MI indicates stronger PAC. For a schematic representation regarding the computation of the MI, refer to Figure 1 in Tort et al., 2010. These computations were performed using the provided PreferredPhase and Pac methods from the TensorPAC library for Python (Combrisson et al., 2020).

Comodulograms

Comodulograms represent the amount of PAC between two ranges of frequencies, used to extract respectively a phase and an amplitude signal. Specifically, we computed the MI between 80% overlapping 1 Hz frequency bands used to compute the phase of the signal (in the theta range), and 90% overlapping 10 Hz frequency bands used to compute the amplitude of the signal (in the gamma range). The resulting distribution of MI values was subsequently plotted as heat maps. An inherent problem with the computation of the MI using simulated data was the lack of frequency components in some frequency bands. Filtering the data within 1 Hz frequency bands thus created a flat signal, resulting in high values of the MI despite the absence of modulation. To overcome this limitation, we added uniform noise to the firing rate signals prior to computing the MI, with an amplitude of approximately 20% of the maximum instantaneous firing rates. The results with and without added noise are presented in Figure 3—figure supplement 2.

Phase dependency of PAC

For a given theta frequency, PAC depends on the phase of the underlying theta oscillation. To identify this relationship and the theta phase that maximizes coupling, we used the PreferredPhase function of the TensorPAC toolbox. More precisely, we applied a Hilbert transform to extract the phase of the theta signal (firing rate band-pass filtered between 3 and 9 Hz) and the amplitude of the gamma signal band-pass filtered within narrow (10 Hz wide) frequency ranges between 20 and 100 Hz with a 90% overlap. For each narrow gamma range, we binned the amplitude with respect to the phase in a similar way as in the calculation of the MI. We obtained a vector of the binned high-frequency amplitudes with respect to the phase of the low-frequency phase, represented as a polar plot as in Figure 3C.

Phase response curves

The PRC of an oscillatory system indicates the phase delay or advancement that follows a single pulse, as a function of the phase at which this input is delivered. To characterize the PRC of our computational model, we applied a single stimulation pulse to CA1 across different phases of the theta rhythm and calculated the resulting change in the theta phase. We split a single theta cycle into intervals of width π/8 radians and applied a single stimulation pulse of a given amplitude. For each case, we ran two simulations: one with a stimulation pulse and one without. Finally, we compared the theta phase 2.5ms post-stimulation and at the same time but in the absence of stimulation. The resulting values for the phase difference Δϕ were plotted against the stimulation phase and are presented in Figure 3D for varying stimulation amplitudes.

Funding Information

This paper was supported by the following grants:

http://dx.doi.org/10.13039/501100009468 Conseil Régional Aquitaine Bordeaux Neurocampus chair of excellence to Fabien B Wagner.

http://dx.doi.org/10.13039/501100006251 Université de Bordeaux Bordeaux Neurocampus chair of excellence to Fabien B Wagner.

http://dx.doi.org/10.13039/501100000781 European Research Council ERC Starting Grant #101040391 (MEMOPROSTHETICS project) to Fabien B Wagner.

Acknowledgements

Simulations presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d’Aquitaine (see https://www.plafrim.fr).

Additional information

Competing interests

Author contributions

Additional files

MDAR checklist

Data availability

All model files and analysis scripts are available at the following repositories under the GNU General Public License v3.0 license on Zenodo (https://www.doi.org/10.5281/zenodo.8354987) and GitHub (https://github.com/NikVard/memstim-hh; copy archived at Vardalakis, 2023).

The following dataset was generated:

Vardalakis N Aussel A Rougier NP Wagner FB 2023 Dataset associated with Vardalakis et al. (2024), eLife Zenodo 10.5281/zenodo.8354987

10.7554/eLife.87356.3.sa0
eLife assessment
Toth Katalin Reviewing Editor University of Ottawa Canada

Convincing
Valuable
This study presents a computational model to explore how neurostimulation could impact hippocampal theta oscillations. The computational model combines a detailed physiologically realistic hippocampus model and an abstract theta oscillator. The study could provide valuable predictions on pathological changes in this network. The modelling is based on convincing approaches that could be improved with experimental validation in future experiments.

10.7554/eLife.87356.3.sa1
Reviewer #1 (Public Review):
Reviewer
In this article, Vardakalis et al. propose a novel model of hippocampal oscillations whereby an external input (emulating the medial septum) can drive theta rhythms. This model displays phase-amplitude coupling of gamma oscillations, as well as theta resetting, which are known features of physiological theta that have been missing in previous models. The end goal proposed by the authors is to have a framework to explore the mechanisms of neurostimulation, which have shown promising applications in pathological conditions, but for which the underlying dynamics remain largely unknown. To reach this objective, the authors implement an existing biophysical model of the hippocampus that is able to generate gamma oscillations, and receives inputs from a set of Kuramoto oscillators to emulate theta drive originating from the medial septum.

Overall, the hypotheses and results are clearly presented and supported by high quality figures. The study is presented in a didactic way, making it easy for a broad audience to understand the significance of the results. The study does present some weaknesses that could easily be addressed by the authors. First, there are some anatomical inaccuracies: line 129 and fig1C, the authors omit medial septum projections to area CA1 (in addition to the entorhinal cortex). Moreover, in addition to CA1, CA3 also provides monosynaptic feedback projections to the medial septum CA3. Finally, an indirect projection from CA1/3 excitatory neurons to the lateral septum, which in turn sends inhibitory projections to the medial septum could be included or mentioned by the authors. This could be of particular relevance to support claims related to effects of neurostimulations, whereby minutious implementation of anatomical data could be key. If not updating their model, the authors could add this point to their limitation section, where they already do a good job of mentioning some limitations of using the EC as a sole oscillatory input to CA1. The authors test conditions of low theta inputs, which they liken to pathological states (line 112). It is not clear what pathology the authors are referring to, especially since a large amount of 'oscillopathies' in the septohippocampal system are associated with decreased gamma/PAC, but not theta oscillations (e.g. Alzheimer's disease conditions). While relevant for the clinical field, there is overall a missed opportunity to explain many experimental accounts with this novel model. Although to this day, clinical use of DBS is mostly restricted to electrical (and thus cell-type agnostic) stimulation, recent studies focusing on mechanisms of neurostimulations have manipulated specific subtypes in the medial septum and observed effects on hippocampal oscillations (e.g. see Muller & Remy, 2017 for review). Focusing stimulations in CA1 is of course relevant for clinical studies but testing mechanistic hypotheses by focusing stimulation on specific cell types could be highly informative. For instance, could the author reproduce recent optogenetic studies (e.g. Bender et al. 2015 for stimulation of fornix fibers; Etter et al., 2019 & Zutshi et al. 2018 for stimulation of septal inhibitory neurons)? Cell specific manipulations should at least be discussed by the authors.

Beyond these weaknesses, this study has a strong utility for researchers wanting to explore hypotheses in the field of neurostimulations. In particular, I see value in such models for exploring more intricate, phase specific effects of continuous, as well as close loop stimulations which are on the rise in systems neuroscience.

10.7554/eLife.87356.3.sa2
Reviewer #2 (Public Review):
Reviewer
Theta-nested gamma oscillations (TNGO) play an important role in hippocampal memory and cognitive processes and are disrupted in pathology. Deep brain stimulation has been shown to affect memory encoding. To investigate the effect of pulsed CA1 neurostimulation on hippocampal TNGO the authors coupled a physiologically realistic model of the hippocampus comprising EC, DG, CA1, and CA3 subfields with an abstract theta oscillator model of the medial septum (MS). Pathology was modeled as weakened theta input from the MS to EC simulating MS neurodegeneration known to occur in Alzheimer's disease. The authors show that if the input from the MS to EC is strong (the healthy state) the model autonomously generates TNGO in all hippocampal subfields while a single neurostimulation pulse has the effect of resetting the TNGO phase. When the MS input strength is weaker the network is quiescent but the authors find that a single CA1 neurostimulation pulse can switch it into the persistent TNGO state, provided the neurostimulation pulse is applied at the peak of the EC theta. If the MS theta oscillator model is supplemented by an additional phase-reset mechanism a single CA1 neurostimulation pulse applied at the trough of EC theta also produces the same effect. If the MS input to EC is weaker still, only a short burst of TNGO is generated by a single neurostimulation pulse. The authors investigate the physiological origin of this burst and find it results from an interplay of CAN and M currents in the CA1 excitatory cells. In this case, the authors find that TNGO can only be rescued by a theta frequency train of CA1 pulses applied at the peak of the EC theta or again at either the peak or trough if the MS oscillator model is supplemented by the phase-reset mechanism.

The main strength of this model is its use of a fairly physiologically detailed model of the hippocampus. The cells are single-compartment models but do include multiple ion channels and are spatially arranged in accordance with the hippocampal structure. This allows the understanding of how ion channels (possibly modifiable by pharmacological agents) interact with system-level oscillations and neurostimulation. The model also includes all the main hippocampal subfields. The other strength is its attention to an important topic, which may be relevant for dementia treatment or prevention, which few modeling studies have addressed.

The work has several weaknesses. First, while investigations of hippocampal neurostimulation are important there are few experimental studies from which one could judge the validity of the model findings. All its findings are therefore predictions. It would be much more convincing to first show the model is able to reproduce some measured empirical neurostimulation effect before proceeding to make predictions. Second, the model is very specific. Or if its behavior is to be considered general it has not been explained why. For example, the model shows bistability between quiescence and TNGO, however what aspect of the model underlies this, be it some particular network structure or particular ion channel, for example, is not addressed. Similarly for the various phase reset behaviors that are found. We may wonder whether a different hippocampal model of TNGO, of which there are many published (for example [1-6]) would show the same effect under neurostimulation. This seems very unlikely and indeed the quiescent state itself shown by this model seems quite artificial. Some indication that particular ion channels, CAN and M are relevant is briefly provided and the work would be much improved by examining this aspect in more detail. In summary, the work would benefit from an intuitive analysis of the basic model ingredients underlying its neurostimulation response properties. Third, while the model is fairly realistic, considerable important factors are not included and in fact, there are much more detailed hippocampal models out there (for example [5,6]). In particular, it includes only excitatory cells and a single type of inhibitory cell. This is particularly important since there are many models and experimental studies where specific cell types, for example, OLM and VIP cells, are strongly implicated in TNGO. Other missing ingredients one may think might have a strong impact on model response to neurostimulation (in particular stimulation trains) include the well-known short-term plasticity between different hippocampal cell types and active dendritic properties. Fourth the MS model seems somewhat unsupported. It is modeled as a set of coupled oscillators that synchronize. However, there is also a phase reset mechanism included. This mechanism is important because it underlies several of the phase reset behaviors shown by the full model. However, it is not derived from experimental phase response curves of septal neurons of which there is no direct measurement. The work would benefit from the use of a more biologically validated MS model.

[1] Hyafil A, Giraud AL, Fontolan L, Gutkin B. Neural cross-frequency coupling: connecting architectures, mechanisms, and functions. Trends in neurosciences. 2015 Nov 1;38(11):725-40.

[2] Tort AB, Rotstein HG, Dugladze T, Gloveli T, Kopell NJ. On the formation of gamma-coherent cell assemblies by oriens lacunosum-moleculare interneurons in the hippocampus. Proceedings of the National Academy of Sciences. 2007 Aug 14;104(33):13490-5.

[3] Neymotin SA, Lazarewicz MT, Sherif M, Contreras D, Finkel LH, Lytton WW. Ketamine disrupts theta modulation of gamma in a computer model of hippocampus. Journal of Neuroscience. 2011 Aug 10;31(32):11733-43.

[4] Ponzi A, Dura-Bernal S, Migliore M. Theta-gamma phase-amplitude coupling in a hippocampal CA1 microcircuit. PLOS Computational Biology. 2023 Mar 23;19(3):e1010942.

[5] Bezaire MJ, Raikov I, Burk K, Vyas D, Soltesz I. Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent CA1 circuit. Elife. 2016 Dec 23;5:e18566.

[6] Chatzikalymniou AP, Gumus M, Skinner FK. Linking minimal and detailed models of CA1 microcircuits reveals how theta rhythms emerge and their frequencies controlled. Hippocampus. 2021 Sep;31(9):982-1002.

10.7554/eLife.87356.3.sa3
Author Response
Vardalakis Nikolaos Author University of Bordeaux Bordeaux France

Aussel Amélie Author Inria Bordeaux - Sud-Ouest Research Centre Bordeaux France

Rougier Nicolas P Author INRIA Bordeaux France

Wagner Fabien B Author University of Bordeaux Bordeaux France

The following is the authors’ response to the original reviews.

Response to reviewers

We would like to thank the reviewers for their feedback. Below we address their comments and have indicated the associated changes in our point-by-point response (blue: answers, red: changes in manuscript).

Reviewer #1:

Overall, the hypotheses and results are clearly presented and supported by high quality figures. The study is presented in a didactic way, making it easy for a broad audience to understand the significance of the results. The study does present some weaknesses that could easily be addressed by the authors.

We thank the reviewer for appreciating our work and providing useful suggestions for improvement.

1. First, there are some anatomical inaccuracies: line 129 and fig1C, the authors omit m.dial septum projections to area CA1 (in addition to the entorhinal cortex). Moreover, in addition to CA1, CA3 also provides monosynaptic feedback projections to the medial septum CA3. Finally, an indirect projection from CA1/3 excitatory neurons to the lateral septum, which in turn sends inhibitory projections to the medial septum could be included or mentioned by the authors. This could be of particular relevance to support claims related to effects of neurostimulations, whereby minutious implementation of anatomical data could be key.

If not updating their model, the authors could add this point to their limitation section, where they already do a good job of mentioning some limitations of using the EC as a sole oscillatory input to CA1.

We acknowledge that our current model strongly simplifies the interconnections between the medial septum and the hippocampal formation, but including more anatomical details is beyond the scope of this manuscript and would be a topic for future work. Nevertheless, we followed the reviewer’s advice to stress this point in our manuscript. First, we moved a paragraph that was initially in the “methods” section to the “results” section (L.141-150 of the revised manuscript):

“Biologically, GABAergic neurons from the medial septum project to the EC, CA3, and CA1 fields of the hippocampus (Toth et al., 1993; Hajós et al., 2004; Manseau et al., 2008; Hangya et al., 2009; Unal et al., 2015; Müller and Remy, 2018). Although the respective roles of these different projections are not fully understood, previous computational studies have suggested that the direct projection from the medial septum to CA1 is not essential for the production of theta in CA1 microcircuits (Mysin et al., 2019). Since our modeling of the medial septum is only used to generate a dynamic theta rhythm, we opted for a simplified representation where the medial septum projects only to the EC, which in turn drives the different fields of the hippocampus. In our model, Kuramoto oscillators are therefore connected to the EC neurons and they receive projections from CA1 neurons (see methods for more details).”

Second, we expanded the corresponding paragraph in the limitation section to discuss this point further (L.398-415 of the revised manuscript):

“We decided to model septal pacemaker neurons projecting to the EC as the main source of hippocampal theta as reported in multiple experimental studies (Buzsáki, 2002; Buzsáki et al., 2003; Hangya et al., 2009). However, experimental findings and previous models have also proposed that direct septal inputs are not essential for theta generation (Wang, 2002; Colgin et al., 2013; Mysin et al., 2019), but play an important role in phase synchronization of hippocampal neurons. Furthermore, the model does not account for the connections between the lateral and medial septum and the hippocampus (Takeuchi et al., 2021). These connections include the inhibitory projections from the lateral to the medial septum and the monosynaptic projections from the hippocampal CA3 field to the lateral septum. An experimental study has highlighted the importance of the lateral septum in regulating the hippocampal theta rhythm (Bender et al., 2015), an area that has not been included in the model. Specifically, theta-rhythmic optogenetic stimulation of the axonal projections from the lateral septum to the hippocampus was shown to entrain theta oscillations and lead to behavioral changes during exploration in transgenic mice. To account for these discrepancies, our model could be extended by considering more realistic connectivity patterns between the medial / lateral septum and the hippocampal formation, including glutamatergic, cholinergic, and GABAergic reciprocal connections (Müller and Remy, 2018), or by considering multiple sets of oscillators each representing one theta generator.”

2. The authors test conditions of low theta inputs, which they liken to pathological states (line 112). It is not clear what pathology the authors are referring to, especially since a large amount of 'oscillopathies' in the septohippocampal system are associated with decreased gamma/PAC, but not theta oscillations (e.g. Alzheimer's disease conditions).

In the manuscript, we referred to “oscillopathies” in a broad sense way as we did not want to overstate the biological implications of the model or the way we modeled pathological states. To our knowledge, several studies have yielded inconsistent results regarding the specific changes in theta or gamma power in Alzheimer’s disease, and the most convincing alteration seems to be the theta-gamma phase-amplitude coupling (PAC) (for review see e.g., Kitchigina, V. F. Alterations of Coherent Theta and Gamma Network Oscillations as an Early Biomarker of Temporal Lobe Epilepsy and Alzheimer’s Disease. Front Integr Neurosci 12, 36 (2018)), as also mentioned by the reviewer.

In this study, the most straightforward way to reduce theta-gamma PAC was to reduce the amplitude of the oscillators’ gain, which affected theta power, gamma power, and theta-gamma PAC (Figure 5 of the revised manuscript). Affecting their synchronization level (i.e., the order parameter) did not affect any of these variables (Figure 5 – Figure Supplement 4).

In order to alter theta-gamma PAC without affecting theta or gamma power, we believe that more complex changes should be performed in the model, likely at the level of individual neurons in the hippocampal formation. For example, cholinergic deprivation has been previously used in a multi-compartment model of the hippocampal CA3 to mimic Alzheimer’s disease and to draw functional implications on the slowing of theta oscillations and the storage of new information (Menschik, E. D. & Finkel, L. H. Neuromodulatory control of hippocampal function: towards a model of Alzheimer’s disease. Artif Intell Med 13, 99–121 (1998)).

This has now been added to the limitations section (L.458-465 of the revised manuscript):

“Finally, we likened conditions of low theta input to pathological states characteristic of oscillopathies such as Alzheimer’s disease, as these conditions disrupted all aspects of theta-gamma oscillations in our model: theta power, gamma power, and theta-gamma PAC (Figure 5). However, it should be noted that changes in theta or gamma power in these pathologies are often unclear, and that the most consistent alteration that has been reported in Alzheimer’s disease is a reduction of theta-gamma PAC (for review, see Kitchigina, 2018). Future work should explore the effects of cellular alterations intrinsic to the hippocampal formation and their impact on theta-gamma oscillations.”

3. While relevant for the clinical field, there is overall a missed opportunity to explain many experimental accounts with this novel model. Although to this day, clinical use of DBS is mostly restricted to electrical (and thus cell-type agnostic) stimulation, recent studies focusing on mechanisms of neurostimulations have manipulated specific subtypes in the medial septum and observed effects on hippocampal oscillations (e.g. see Muller & Remy, 2017 for review). Focusing stimulations in CA1 is of course relevant for clinical studies but testing mechanistic hypotheses by focusing stimulation on specific cell types could be highly informative. For instance, could the author reproduce recent optogenetic studies (e.g. Bender et al. 2015 for stimulation of fornix fibers; Etter et al., 2019 & Zutshi et al. 2018 for stimulation of septal inhibitory neurons)? Cell specific manipulations should at least be discussed by the authors.

We acknowledge the importance of cell-type-specific manipulation in the septo-hippocampal circuitry. However, our model was designed to study neurostimulation protocols that affect the hippocampal formation, not the medial septum, which is why only the hippocampal formation is composed of biophysically realistic (i.e., conductance-based) neuronal models. To replicate the various studies mentioned by the reviewer (which are all very relevant), we would need to implement a biophysical model of the medial septum, which would be an entirely new project.

Nevertheless, we can use the existing model to replicate optogenetic studies that induced gamma oscillations in excitatory-inhibitory circuits, using either ramped photostimulation targeting excitatory neurons (Adesnik et al., 2010; Akam et al., 2012; Lu et al., 2015), or pulsed stimulation driving inhibitory cells in the gamma range (Cardin et al., 2009; Iaccarino et al., 2016). In fact, such approaches have been demonstrated not just in the hippocampus but also in the neocortex, and represent a hallmark of local excitatory-inhibitory circuits. To account for these experimental results and replicate them, we have added 4 new figures (Figure 2 and its 3 figure supplements) and an extensive section in the results part (L.151-217 of the revised manuscript):

“From a conceptual point of view, our model is thus composed of excitatory-inhibitory (E-I) circuits connected in series, with a feedback loop going through a population of coupled phase oscillators. In the next sections, we first describe the generation of gamma oscillations by individual E-I circuits (Figure 2), and illustrate their behavior when driven by an oscillatory input such as theta oscillations (Figure 3). We then present a thorough characterization of the effects of theta input and stimulation amplitude on theta-nested gamma oscillations (Figure 4 and Figure 5). Finally, we present some results on the effects of neurostimulation protocols for restoring theta-nested gamma oscillations in pathological states (Figure 6 and Figure 7).

Generation of gamma oscillations by E-I circuits

It is well-established that a network of interconnected pyramidal neurons and interneurons can give rise to oscillations in the gamma range, a mechanism termed pyramidal-interneuronal network gamma (PING) (Traub et al., 2004; Onslow et al., 2014; Segneri et al., 2020;). This mechanism has been observed in several optogenetic studies with gradually increasing light intensity (i.e., under a ramp input) affecting multiple different circuits, such as layer 2-3 pyramidal neurons of the mouse somatosensory cortex (Adesnik et al., 2010), the CA3 field of the hippocampus in rat in vitro slices (Akam et al., 2012), and in the non-human primate motor cortex (Lu et al., 2015). In all cases, gamma oscillations emerged above a certain threshold in terms of photostimulation intensity, and the frequency of these oscillations was either stable or slightly increased when increasing the intensity further. We sought to replicate these findings with our elementary E-I circuits composed of single-compartment conductance-based neurons driven by a ramping input current (Figure 2 and Figure S2). As an example, all the results in this section will be shown for an E-I circuit that has similar connectivity parameters as the CA1 field of the hippocampus in our complete model (see section “Hippocampal formation: inputs and connectivity” in the methods).

For low input currents provided to both neuronal populations, only the highly-excitable interneurons were activated (Figure 2A). For a sufficiently high input current (i.e., a strong input that could overcome the inhibition from the fast-spiking interneurons), the pyramidal neurons started spiking as well. As the amplitude of the input increased, the activity of the both neuronal populations became synchronized in the gamma range, asymptotically reaching a frequency of about 60 Hz (Figure 2A bottom panel). Decoupling the populations led to the abolition of gamma oscillations (Figure 2B), as neuronal activity was determined solely by the intrinsic properties of each cell. Interestingly, when the ramp input was provided solely to the excitatory population, we observed that the activity of the pyramidal neurons preceded the activity of the inhibitory neurons, while still preserving the emergence of gamma oscillations (Figure S2 A). As expected, decoupling the populations also abolished gamma oscillations, with the excitatory neurons spiking a frequency determined by their intrinsic properties and the inhibitory population remaining silent (Figure S2B).

To further characterize the intrinsic properties of individual inhibitory and excitatory neurons, we derived their input-frequency (I-F) curves, which represent the firing rate of individual neurons in response to a tonic input (Figure S3A). We observed that for certain input amplitudes, the firing rates of both types of neurons was within the gamma range. Interestingly, in the absence of noise, each population could generate by itself gamma oscillations that were purely driven by the input and determined by the intrinsic properties of the neurons (Figure S3B). Adding stochastic Gaussian noise in the membrane potential disrupted these artificial oscillations in decoupled populations (Figure S3C). All subsequent simulations were run with similar noise levels to prevent the emergence of artificial gamma oscillations.

Another potent way to induce gamma oscillations is to drive fast-spiking inhibitory neurons using pulsed optogenetic stimulation at gamma frequencies, a strategy that has been used both in the neocortex (Cardin et al., 2009) and hippocampal CA1 (Iaccarino et al., 2016). In particular, Cardin and colleagues systematically investigated the effect of driving either excitatory or fast-spiking inhibitory neocortical neurons at frequencies between 10 and 200 Hz (Cardin et al., 2009). They showed that fast-spiking interneurons are preferentially entrained around 40-50 Hz, while excitatory neurons respond better to lower frequencies. To verify the behavior of our model against these experimental data, we simulated pulsed optogenetic stimulation as an intracellular current provided to our reduced model of a single E-I circuit. Stimulation was applied at frequencies between 10 and 200 Hz to excitatory cells only, to inhibitory cells only, or to both at the same time (Figure S4). The population firing rates were used as a proxy for the local field potentials (LFP), and we computed the relative power in a 10-Hz band centered around the stimulation frequency, similarly to the method proposed in (Cardin et al., 2009). When presented with continuous stimulation across a range of frequencies in the gamma range, interneurons showed the greatest degree of gamma power modulation (Figure S4). Furthermore, when the stimulation was delivered to the excitatory population, the relative power around the stimulation frequency dropped significantly in frequencies above 10 Hz, similar to the reported experimental data (Cardin et al., 2009). The main difference between our simulation results and these experimental data is the specific frequencies at which fast-spiking interneurons showed resonance, which was slow gamma around 40 Hz in the mouse barrel cortex and fast gamma around 90 Hz in our model. This could be attributed to several factors, such as differences in the cellular properties between cortical and hippocampal fast-spiking interneurons, or the differences between the size of the populations and their relevant connectivity in the cortex and the hippocampus.”

Author response image 1. Figure 2.

Emergence of gamma oscillations in coupled excitatory-inhibitory populations under ramping input to both populations. (A) Two coupled populations of excitatory pyramidal neurons (NE = 1000) and inhibitory interneurons (NI = 100) are driven by a ramping current input (0 nA to 1 nA) for 5 s. As the input becomes stronger, oscillations start to emerge (shaded green area), driven by the interactions between excitatory and inhibitory populations. The green inset shows the raster plot (neuronal spikes across time) of the two populations during the green shaded period (red for inhibitory; blue for excitatory). When the input becomes sufficiently strong (shaded magenta area), the populations become highly synchronized and produce oscillations in the gamma range (at approximately 50 Hz). The spectrogram (bottom panel) shows the power of the instantaneous firing rate of the pyramidal population as a function of time and frequency. It reveals the presence of gamma oscillations that emerge around 2s and increase in frequency until 4 s, when they settle at approximately 60 Hz. (B) Similar depiction as in panel A. with the pyramidal-interneuronal populations decoupled. The absence of coupling leads to the abolition of gamma oscillations, each cell spiking activity being driven by its own inputs and intrinsic properties.

Author response image 2. Figure S2 (Figure 2 – Figure Supplement 1).

Emergence of gamma oscillations in coupled excitatoryinhibitory populations under ramping input to the excitatory population. Similar representation as in Figure 2, but with the input provided only to the excitatory population. All conclusions remain the same. In addition, the inhibitory population does not show any spiking activity in the decoupled case.

Author response image 3. Figure S3 (Figure 2 – Figure Supplement 2).

Cell-intrinsic spiking activity in decoupled excitatory and inhibitory populations under ramping input. (A) Input-Frequency (I-F) curves for excitatory cells (left panel; pyramidal neurons with ICAN) and inhibitory cells (right panel; interneurons, fast-spiking) used in the model. Above a certain tonic input (around 0.35 nA for excitatory and 0.1 nA for inhibitory neurons), neurons can spike in the gamma range. (B) Raster plot showing the spiking activity of excitatory (blue, NE = 1000) and inhibitory (red, NI = 100) neurons in decoupled populations under ramping input (top trace) and in the absence of noise in the membrane potential. Despite random initial conditions across neurons, oscillations emerge in both populations due to the intrinsic properties of the cells, with a frequency that is predicted by the respective I-F curves (panel A.). (C) Similar representation as panel B. but with the addition of stochastic noise in the membrane potential of each neuron. The presence of noise disrupts the emergence of oscillations in these decoupled populations.

Author response image 4. Figure S3 (Figure 2 – Figure Supplement 2).

Cell-intrinsic spiking activity in decoupled excitatory and inhibitory populations under ramping input. (A) Input-Frequency (I-F) curves for excitatory cells (left panel; pyramidal neurons with ICAN) and inhibitory cells (right panel; interneurons, fast-spiking) used in the model. Above a certain tonic input (around 0.35 nA for excitatory and 0.1 nA for inhibitory neurons), neurons can spike in the gamma range. (B) Raster plot showing the spiking activity of excitatory (blue, NE = 1000) and inhibitory (red, NI = 100) neurons in decoupled populations under ramping input (top trace) and in the absence of noise in the membrane potential. Despite random initial conditions across neurons, oscillations emerge in both populations due to the intrinsic properties of the cells, with a frequency that is predicted by the respective I-F curves (panel A.). (C) Similar representation as panel B. but with the addition of stochastic noise in the membrane potential of each neuron. The presence of noise disrupts the emergence of oscillations in these decoupled populations.

Beyond these weaknesses, this study has a strong utility for researchers wanting to explore hypotheses in the field of neurostimulations. In particular, I see value in such models for exploring more intricate, phase specific effects of continuous, as well as close loop stimulations which are on the rise in systems neuroscience.

We thank the reviewer for this appreciation of our work and its future perspectives.

Recommendations For The Authors:

Line 144, the authors mention that their MI values are erroneous in absence of additive noise - could this be due to the non-sinusoidal nature of the phase signal recorded, and be fixed by upscaling model size?

We thank the reviewer for this question and suggestion. The main reason behind the errors in the computation of the MI lies in the complete absence of oscillations at specific frequencies. Filtered signals within specific bands produced a power of 0 (or extremely low values), as seen in the power spectral densities. In such cases, the phase signal was not mathematically defined, but the toolbox we used to compute it still returned a numerical result that was inaccurate (for more details on the computation of the MI see Tort et al., 2010). To mitigate this numerical artefact, we decided to add uniform noise in the computed firing rates. This strategy is illustrated on Figure S6 (Figure 3 – Figure Supplement 2), which we have copied below for reference. Alternative approaches could probably have been used, such as increasing the noise in the membrane potential so that neurons would start spiking with firing rates that show more realistic power spectra, even in the absence of external inputs.

Author response image 5. Figure S6 (Figure 3 – Figure Supplement 2).

Quantification of PAC with and without noise. (A) Quantifying PAC in the absence of noise produced inaccurate identification of the coupled frequency bands, due to the complete absence of oscillations at some frequencies. All analyses are based on the CA1 firing rates (top traces) during a representative simulation. Power spectral densities of these firing rates (left) indicate that some frequencies have a power of 0. PAC of the excitatory population was assessed using two graphical representations, the polar plot (middle) and comodulogram (right), and quantified using the MI. The comodulogram was calculated by computing the MI across 80% overlapping 1-Hz frequency bands in the theta range and across 90% overlapping 10-Hz frequency bands in the gamma range and subsequently plotted as a heat map. In the absence of noise, a slow theta frequency centered around 5 Hz is found to modulate a broad range of gamma frequencies between 40 and 100 Hz. The value indicated on the comodulogram indicates the average MI in the 3-9 Hz theta range and 40-80 Hz gamma range. As in Figure 2, the polar plot represents the amplitude of gamma oscillations (averaged across all theta cycles) at each phase of theta (theta range: 3-9 Hz, phase indicated as angular coordinate) and for different gamma frequencies (radial coordinate, binned in 1-Hz ranges). (B) Adding uniform noise to the firing rate (with an amplitude ranging between 15 and 25% of the maximum firing rate) improved the identification of the coupled frequency bands. In this case, the slower theta frequency centered around 5 Hz modulates a gamma band located between 45 and 75 Hz.

Reviewer #2:

The main strength of this model is its use of a fairly physiologically detailed model of the hippocampus. The cells are single-compartment models but do include multiple ion channels and are spatially arranged in accordance with the hippocampal structure. This allows the understanding of how ion channels (possibly modifiable by pharmacological agents) interact with system-level oscillations and neurostimulation. The model also includes all the main hippocampal subfields. The other strength is its attention to an important topic, which may be relevant for dementia treatment or prevention, which few modeling studies have addressed. The work has several weaknesses.

We thank the reviewer for appreciating our detailed description of the hippocampal formation and the focus on neurostimulation applications that aim at treating oscillopathies, especially dementia.

1. First, while investigations of hippocampal neurostimulation are important there are few experimental studies from which one could judge the validity of the model findings. All its findings are therefore predictions. It would be much more convincing to first show the model is able to reproduce some measured empirical neurostimulation effect before proceeding to make predictions.

We acknowledge that the results presented in Figures 4-7 of the revised manuscript cannot be compared to existing experimental data, and are therefore purely predictive. Future experimental work is needed to verify these predictions.

Yet, we would also like to stress that the motivation behind this project was the inadequacy of previous models of theta-nested gamma oscillations (Onslow et al., 2014; Aussel et al., 2018; Segneri et al., 2020) to account for the mechanism of theta phase reset that occurs during electrical stimulation of the fornix or perforant path (Williams and Givens, 2003). Since we could not use these previous models to study the effects of neurostimulation on theta-nested gamma oscillations, we had to modify them to account for a dynamical theta input, which is the main methodological novelty that is reported in our manuscript (Figures 1 and 3 of the revised manuscript).

Despite the scarcity of experimental studies that could confirm the full model, we sought to replicate a few experimental findings that employed optogenetic stimulation to induce gamma oscillations in individual excitatory-inhibitory circuits. Although not specific to the hippocampus, these studies have shown that gamma oscillations can be induced using either ramped photostimulation targeting excitatory neurons (Adesnik et al., 2010; Akam et al., 2012; Lu et al., 2015), or pulsed stimulation driving inhibitory cells in the gamma range (Cardin et al., 2009; Iaccarino et al., 2016). To account for these experimental results and replicate them, we have added 4 new figures (Figure 2 and its 3 figure supplements) and an extensive section in the results part (L.141-217 of the revised manuscript). The added section and related figures are indicated in our response to reviewer 1, comment 3 (p 2-7).

2.1. Second, the model is very specific. Or if its behavior is to be considered general it has not been explained why.

Although the spatial organization and cellular details of the model are indeed very specific, its general behavior, i.e., the production of theta-nested gamma oscillations and theta phase reset, are common to any excitatory-inhibitory circuit interconnected with Kuramoto oscillators. To illustrate this point, we have generalized our approach to the neural mass model developed by Onslow and colleagues (Onslow ACE, Jones MW, Bogacz R. A Canonical Circuit for Generating Phase-Amplitude Coupling. PLoS ONE. 2014 Aug; 9(8):e102591). These results are represented in a new supplementary figure (Figure3 – Figure Supplement 4), and briefly described in a new paragraph of the results section (L.262-268 of the revised manuscript):

“Importantly, our approach is generalizable and can be applied to other models producing theta-nested gamma oscillations. For instance, we adapted the neural mass model by Onslow and colleagues (Onslow et al., 2014), replaced the fixed theta input by a set of Kuramoto oscillators, and demonstrated that it could also generate theta phase reset in response to single-pulse stimulation (Figure S8). These results illustrate that the general behavior of our model is not specific to the tuning of individual parameters in the conductancebased neurons, but follows general rules that are captured by the level of abstraction of the Kuramoto formalism.”

Author response image 6. Figure S8 (Figure 3 – Figure Supplement 4).

A neural mass model of coupled excitatory and inhibitory neurons driven by Kuramoto oscillators generates theta-nested gamma oscillations and theta phase reset. (A) Two coupled neural masses (one excitatory and one inhibitory) driven by Kuramoto oscillators, which represent a dynamical oscillatory drive in the theta range, were used to implement a neural mass equivalent to our conductance-based model represented in Figure 1. Neural masses were modeled using the WilsonCowan formalism, with parameters adapted from Onslow et al. (2014) (WEE = 4.8, WEI = WIE = 4, WII = 0). (B) The normalized population firing rates exhibit theta-nested gamma oscillations (middle and bottom panels) in response to the dynamic theta rhythm (top panel). A stimulation pulse delivered at the descending phase of the rhythm to both populations (marked by the inverted red triangle) produces a robust theta phase reset, similarly to Figure 3A.

This simplified model is described in more details in the methods (L.694-710 of the revised manuscript). Additionally, the generation of gamma oscillations by individual excitatory-inhibitory circuits is now described in details in the added section “Generation of gamma oscillations by E-I circuits” (L.159-217 of the revised manuscript), which has already been discussed in our response to reviewer 1, comment 3 (p 2-7).

2.2. For example, the model shows bistability between quiescence and TNGO, however what aspect of the model underlies this, be it some particular network structure or particular ion channel, for example, is not addressed.

We thank the reviewer for mentioning this point, which we have now addressed. The “bistable” behavior that we reported occurs for values of the theta input that are just below the threshold to induce selfsustained theta-gamma oscillations (Figure 5 of the revised manuscript, point B). Moreover, the presence of the Calcium-Activated-Nonspecific (CAN) cationic channel, which is expressed by pyramidal neurons in the entorhinal cortex, CA3, and CA1 fields of the hippocampus, is necessary for this behavior to occur. Indeed, abolishing CAN channels in all areas of the model suppresses this behavior. We have now addressed this point in a new supplementary figure (Figure 5 – Figure Supplement 4) and a short description in the text (L.287-303 of the revised manuscript).

“In the presence of dynamic theta input, the effects of single-pulse stimulation depended both on theta input amplitude and stimulation amplitude, highlighting different regimes of network activity (Figure 5 and Figure S9, Figure S10, Figure S11). For low theta input, theta-nested gamma oscillations were initially absent and could not be induced by stimulation (Figure 5A). At most, the stimulation could only elicit a few bursts of spiking activity that faded away after approximately 250 ms, similar to the rebound of activity seen in the absence of theta drive. For increasing theta input, the network switched to an intermediate regime: upon initialization at a state with no spiking activity, it could be kicked to a state with self-sustained theta-nested gamma oscillations by a single stimulation pulse of sufficiently high amplitude (Figure 5B). This regime existed for a range of septal theta inputs located just below the threshold to induce self-sustained theta-gamma oscillations without additional stimulation, as characterized by the post-stimulation theta power, gamma power, and theta-gamma PAC (Figure 5D). Removing CAN currents from all areas of the model abolished this behavior (Figure S12), which is interesting given the role of this current in the multistability of EC neurons (Egorov et al., 2002; Fransen et al., 2006) and in the intrinsic ability of the hippocampus to generate thetanested gamma oscillations (Giovannini et al., 2017). For the highest theta input, the network became able to spontaneously generate theta-nested gamma oscillations, even when initialized at a state with no spiking activity and without additional neurostimulation (Figure 5C).”

Author response image 7. Figure S12 (Figure 5 – Figure Supplement 4).

CAN currents are necessary for the production of selfsustained theta-gamma oscillations in response to single-pulse stimulation. (A) Same as Figure 5B. (B) Similar simulation as panel A., but without the presence of CAN currents in the EC, CA3 and CA1 fields of the hippocampus. Removing CAN currents from the model abolishes self-sustained theta-nested gamma oscillations in response to a single stimulation pulse (for the parameters represented in Figure 5, point B).

Furthermore, we realized that the terminology “bistable” may not be justified as we could not perform a systematic bifurcation analysis, which is typically carried out in simpler neural mass models (e.g., Onslow et al., 2014; Segneri et al., 2020). Therefore, we decided to rephrase the sentences about “bistability” to keep a more general terminology. The following sentences were revised:

L.20-23: “We showed that, for theta inputs just below the threshold to induce self-sustained theta-nested gamma oscillations, a single stimulation pulse could switch the network behavior from non-oscillatory to a state producing sustained oscillations.”

L.305-309: “Based on the above analyses, we considered two pathological states: one with a moderate theta input (i.e., moderately weak projections from the medial septum to the EC) that allowed the initiation of selfsustained oscillations by single stimulation pulses (Figure 5, point B), and one with a weaker theta input characterized by the complete absence of self-sustained oscillations even following transient stimulation(Figure 5, point A).”

L.316-317: “In the case of a moderate theta input and in the presence of phase reset, delivering a pulse at either the peak or trough of theta could induce theta-nested gamma oscillations (Figure 6A and 6C).”

L.353-357: “A very interesting finding concerns the behavior of the model in response to single-pulse stimulation for certain values of the theta amplitude (Figure5). For low theta amplitudes, a single stimulation pulse was capable of switching the network behavior from a state with no spiking activity to one with prominent theta-nested gamma oscillations. Whether such an effect can be induced in vivo in the context of memory processes remains an open question.”

2.3. Similarly for the various phase reset behaviors that are found.

We would like to clarify the fact that the observed phase reset curves (reported in Figure 3D) are a direct consequence of the choice of an appropriate phase response function for the Kuramoto oscillators representing the medial septum. This choice is inspired by experimentally measured phase response curves from CA3 neurons. These aspects are described briefly in the introduction and in more details in the methods, as indicated below:

L.101: “This new hybrid dynamical model could generate both theta-nested gamma oscillations and theta phase reset, following a particular phase response curve (PRC) inspired by experimental literature (Lengyel et al., 2005; Akam et al., 2012; Torben-Nielsen et al., 2010).”

L.528-537: “Hereafter, we call the term Z(θ) the phase response function, to distinguish it from the PRC obtained from experimental data or simulations (see section below "Data Analysis", "Phase Response Curve"). Briefly, the PRC of an oscillatory system indicates the phase delay or advancement that follows a single pulse, as a function of the phase at which this input is delivered. The phase response function Z(θ) was chosen to mimic as well as possible experimental PRCs reported in the literature (Lengyel et al., 2005; Kwag and Paulsen, 2009; Akam et al., 2012). These PRCs appear biphasic and show a phase advancement (respectively delay) for stimuli delivered in the ascending (respectively descending) slope of theta. To accurately model this behavior, we used the following equation for the phase response function, where ϕpeak represents the phase at which the theta rhythm reaches its maximum and the parameter ϕoffset controls the desired phase offset from the peak:

Author response image 8. On the figure below, we illustrate the phase response curves of CA3 neurons measured by Lengyel et al., 2005 (A), and compare it with our simulated phase response curves (B).

Note that the conventions for phase advance and phase delay are reversed between the two panels.

Finally, we would like to acknowledge that the model “is not derived from experimental phase response curves of septal neurons of which there is no direct measurement”, as mentioned by the reviewer in their comment 4 below. Despite the lack of experimental data specific to medial septum neurons, we argue that this phase response function is the only one that mathematically supports the generation of self-sustained theta-nested gamma oscillations in our current model. This statement is illustrated by Figure S7 (Figure 3 – Figure Supplement 3) and is mentioned in the results (L.249-261 of the revised manuscript):

We modeled this behavior by a specific term (which we called the phase response function) in the general equation of the Kuramoto oscillators (see methods, Equation 1). Importantly, introducing a phase offset in the phase response function disrupted theta-nested gamma oscillations (Figure S7), which suggests that the septohippocampal circuitry must be critically tuned to be able to generate such oscillations. The strength of phase reset could also be adjusted by a gain that was manually tuned. In the presence of the physiological phase response function and of a sufficiently high reset gain, a single stimulation pulse delivered to all excitatory and inhibitory CA1 neurons could reset the phase of theta to a value close to its peaks (Figure 3A). We computed the PRC of our simulated data for different stimulation amplitudes and validated that our neuronal network behaved according to the phase response function set in our Kuramoto oscillators (Figure 3D). It should be noted that including this phase reset mechanism affected the generated theta rhythm even in the absence of stimulation, extending the duration of the theta peak and thereby slowing down the frequency of the generated theta rhythm.

Author response image 9. Figure S7 (Figure 3 – Figure Supplement 3).

Network behavior generated by Kuramoto oscillators with nonphysiological phase response functions. Each panel is similar to Figure 3A, but with a different offset added to the phase response function of the Kuramoto oscillators (see methods, Equation 4). The center frequency was set to 6 Hz in all of these simulations. Overall, theta oscillations in these cases are less sinusoidal and show more abrupt phase changes than in the physiological case. (A) A phase offset of −π/2 leads to an overall theta oscillation of 4 Hz, with a second peak following the main theta peak. (B) A phase offset of +π/2 reduces the peak of theta, resetting the rhythm to the middle of the ascending phase. (C) A phase offset of π or−π leads to the CA1 output resetting the theta rhythm to the trough of theta.

2.4. We may wonder whether a different hippocampal model of TNGO, of which there are many published (for example [1-6]) would show the same effect under neurostimulation. This seems very unlikely […]

[1] Hyafil A, Giraud AL, Fontolan L, Gutkin B. Neural cross-frequency coupling: connecting architectures, mechanisms, and functions. Trends in neurosciences. 2015 Nov 1;38(11):725-40.

[2] Tort AB, Rotstein HG, Dugladze T, Gloveli T, Kopell NJ. On the formation of gamma-coherent cell assemblies by oriens lacunosum-moleculare interneurons in the hippocampus. Proceedings of the National Academy of Sciences. 2007 Aug 14;104(33):13490-5.

[3] Neymotin SA, Lazarewicz MT, Sherif M, Contreras D, Finkel LH, Lytton WW. Ketamine disrupts theta modulation of gamma in a computer model of hippocampus. Journal of Neuroscience. 2011 Aug 10;31(32):11733-43.

[4] Ponzi A, Dura-Bernal S, Migliore M. Theta-gamma phase-amplitude coupling in a hippocampal CA1 microcircuit. PLOS Computational Biology. 2023 Mar 23;19(3):e1010942.

[5] Bezaire MJ, Raikov I, Burk K, Vyas D, Soltesz I. Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent CA1 circuit. Elife. 2016 Dec 23;5:e18566.

[6] Chatzikalymniou AP, Gumus M, Skinner FK. Linking minimal and detailed models of CA1 microcircuits reveals how theta rhythms emerge and their frequencies controlled. Hippocampus. 2021 Sep;31(9):982-1002.

The highlighted publications, while very important in their findings regarding theta-gamma phase-amplitude coupling, focused on specific subfields of the hippocampus. In our work, we aimed to develop a model that includes the different anatomical divisions of the hippocampal formation, while still exhibiting theta-nested gamma oscillations, which is why we decided to expand the model by Aussel et al. (2018). Exploring the behavior of all these different hippocampal models under neurostimulation is beyond the scope of the current manuscript.

Nevertheless, we have added a new figure (Figure 3 – Figure Supplement 4) showing an adaptation of our modeling approach to a generic neural mass model of theta-nested gamma oscillations (Onslow et al., 2014), which illustrates the generalizability of our findings and is described in details in our response to comment 2.1. Moreover, we have further addressed the comments of the reviewers regarding bistability and phase response curves in our responses to comments 2.2 and 2.3.

Furthermore, we have added references to all 6 of these publications in the revised version of the manuscript:

L.43-50: Moreover, the modulation of gamma oscillations by the phase of theta oscillations in hippocampal circuits, a phenomenon termed theta-gamma phase-amplitude coupling (PAC), correlates with the efficacy of memory encoding and retrieval (Jensen and Colgin, 2007; Tort et al., 2009; Canolty and Knight, 2010; Axmacher et al., 2010; Fell and Axmacher, 2011; Lisman and Jensen, 2013; Lega et al., 2016). Experimental and computational work on the coupling between oscillatory rhythms has indicated that it originates from different neural architectures and correlates with a range of behavioral and cognitive functions, enabling the long-range synchronization of cortical areas and facilitating multi-item encoding in the context of memory (Hyafil et al., 2015)."

L.415-426: “In terms of neuronal cell types, we also made an important simplification by considering only basket cells as the main class of inhibitory interneuron in the whole hippocampal formation. However, it should be noted that many other types of interneurons exist in the hippocampus and have been modeled in various works with higher computational complexity (e.g., Bezaire et al., 2016; Chatzikalymniou et al., 2021). Among these various interneurons, oriens-lacunosum moleculare (OLM) neurons in the CA1 field have been shown to play a crucial role in synchronizing the activity of pyramidal neurons at gamma frequencies (Tort et al., 2007), and in generating theta-gamma PAC (e.g., Neymotin et al., 2011; Ponzi et al., 2023). Additionally, these cells may contribute to the formation of specific phase relationships within CA1 neuronal populations, through the integration between inputs from the medial septum, the EC, and CA3 (Mysin et al., 2019). Future work is needed to include more diverse cell types and detailed morphologies modeled through multiple compartments.”

2.5. […] and indeed the quiescent state itself shown by this model seems quite artificial.

We would like to clarify the fact that the “quiescent state” mentioned by the reviewer is a simply a state where the theta input is too low to induce theta-nested gamma oscillations. In this regime, neurons are active only due to the noise term in the membrane potential, which was adjusted based on Figure S3 (Figure 2 – Figure Supplement 2, shown below), at the minimal level needed to disrupt artificial synchronization in decoupled populations. For an input of 0 nA, we acknowledge that this network is indeed fully quiescent (i.e., does not show any spiking activity). However, as soon as the input increases, spontaneous spiking activity starts to appear with an average firing rate that depends on the input amplitude and is characterized by the input-frequency curves (panel A.). Please note that adding more noise could eliminate the observed quiescence in the absence of any input, but that it would not affect qualitatively the reported results.

Author response image 10. Figure S3 (Figure 2 – Supplement 2).

Cell-intrinsic spiking activity in decoupled excitatory and inhibitory populations under ramping input. (A) Input-Frequency (I-F) curves for excitatory cells (left panel; pyramidal neurons with ICAN) and inhibitory cells (right panel; interneurons, fast-spiking) used in the model. Above a certain tonic input (around 0.35 nA for excitatory and 0.1 nA for inhibitory neurons), neurons can spike in the gamma range. (B) Raster plot showing the spiking activity of excitatory (blue, NE = 1000) and inhibitory (red, NI = 100) neurons in decoupled populations under ramping input (top trace) and in the absence of noise in the membrane potential. Despite random initial conditions across neurons, oscillations emerge in both populations due to the intrinsic properties of the cells, with a frequency that is predicted by the respective IF curves (panel A.). (C) Similar representation as panel B. but with the addition of stochastic noise in the membrane potential of each neuron. The presence of noise disrupts the emergence of oscillations in these decoupled populations.

2.6. Some indication that particular ion channels, CAN and M are relevant is briefly provided and the work would be much improved by examining this aspect in more detail.

We thank the reviewer for acknowledging the importance of these ion channels. We have now added a new supplementary figure (Figure 5 – Figure Supplement 4), which is described in more details in our response to comment 2.2 and illustrates the role of the CAN current in the generation of theta-nested gamma oscillations following a single stimulation pulse. Moreover, we would like to stress that the impact of CAN currents in the ability of the hippocampus to generate theta-nested gamma oscillations intrinsically, i.e., in the absence of persistent external input, has already been investigated in details by a previous computational study cited in our manuscript (Giovannini F, Knauer B, Yoshida M, Buhry L. The CAN-In network: A biologically inspired model for self-sustained theta oscillations and memory maintenance in the hippocampus. Hippocampus. 2017 Apr;809 27(4):450–463).

2.7. In summary, the work would benefit from an intuitive analysis of the basic model ingredients underlying its neurostimulation response properties.

We thank the reviewer for this suggestion. By addressing the reviewer’s previous comments (reviewer 2, comments 2.1 and 2.2), which overlap partly with the first reviewer (reviewer 1, comment 3), we believe we have improved the manuscript and have provided key information related to the way the model responds to neurostimulation.

3.1. Third, while the model is fairly realistic, considerable important factors are not included and in fact, there are much more detailed hippocampal models out there (for example [5,6]). In particular, it includes only excitatory cells and a single type of inhibitory cell. This is particularly important since there are many models and experimental studies where specific cell types, for example, OLM and VIP cells, are strongly implicated in TNGO.

[5] Bezaire MJ, Raikov I, Burk K, Vyas D, Soltesz I. Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent CA1 circuit. Elife. 2016 Dec 23;5:e18566.

[6] Chatzikalymniou AP, Gumus M, Skinner FK. Linking minimal and detailed models of CA1 microcircuits reveals how theta rhythms emerge and their frequencies controlled. Hippocampus. 2021 Sep;31(9):982-1002.

We thank the reviewer for pointing out these interesting avenues for future studies. As indicated in previous responses (reviewer 1, comment 1; reviewer 2, comment 2.4), we have added several paragraphs to discuss these limitations, the rationale behind our simplifications, and potential improvements. In particular, we have added the following paragraphs to discuss our simplifications in terms of connectivity and cell types:

Anatomical connectivity:

L.141-150: “Biologically, GABAergic neurons from the medial septum project to the EC, CA3, and CA1 fields of the hippocampus (Toth et al., 1993; Hajós et al., 2004; Manseau et al., 2008; Hangya et al., 2009; Unal et al., 2015; Müller and Remy, 2018). Although the respective roles of these different projections are not fully understood, previous computational studies have suggested that the direct projection from the medial septum to CA1 is not essential for the production of theta in CA1 microcircuits (Mysin et al., 2019). Since our modeling of the medial septum is only used to generate a dynamic theta rhythm, we opted for a simplified representation where the medial septum projects only to the EC, which in turn drives the different subfields of the hippocampus. In our model, Kuramoto oscillators are therefore connected to the EC neurons and they receive projections from CA1 neurons (see methods for more details).”

Cell types:

L.415-426: “In terms of neuronal cell types, we also made an important simplification by considering only basket cells as the main class of inhibitory interneuron in the whole hippocampal formation. However, it should be noted that many other types of interneurons exist in the hippocampus and have been modeled in various works with higher computational complexity (e.g., Bezaire et al., 2016; Chatzikalymniou et al., 2021). Among these various interneurons, oriens-lacunosum moleculare (OLM) neurons in the CA1 field have been shown to play a crucial role in synchronizing the activity of pyramidal neurons at gamma frequencies (Tort et al., 2007), and in generating theta-gamma PAC (e.g., Neymotin et al., 2011; Ponzi et al., 2023). Additionally, these cells may contribute to the formation of specific phase relationships within CA1 neuronal populations, through the integration between inputs from the medial septum, the EC, and CA3 (Mysin et al., 2019). Future work is needed to include more diverse cell types and detailed morphologies modeled through multiple compartments.”

3.2. Other missing ingredients one may think might have a strong impact on model response to neurostimulation (in particular stimulation trains) include the well-known short-term plasticity between different hippocampal cell types and active dendritic properties.

We agree with the reviewer that plasticity mechanisms are important to include in future work, which we had already mentioned in the limitations section of the manuscript:

L.436-443: “Importantly, we did not consider learning through synaptic plasticity, even though such mechanisms could drastically modify synaptic conduction for the whole network (Borges et al., 2017). Even more interestingly, the inclusion of spike-timing-dependent plasticity would enable the investigation of stimulation protocols aimed at promoting LTP, such as theta-burst stimulation (Larson et al., 2015). This aspect would be of uttermost importance to make a link with memory encoding and retrieval processes (Axmacher et al., 2006; Tsanov et al., 2009; Jutras et al., 2013) and with neurostimulation studies for memory improvement (Titiz et al., 2017; Solomon et al., 2021).”

4. Fourth the MS model seems somewhat unsupported. It is modeled as a set of coupled oscillators that synchronize. However, there is also a phase reset mechanism included. This mechanism is important because it underlies several of the phase reset behaviors shown by the full model. However, it is not derived from experimental phase response curves of septal neurons of which there is no direct measurement. The work would benefit from the use of a more biologically validated MS model.

We would like to confirm that the phase reset mechanism is indeed at the core of using Kuramoto oscillators to model a particular system. For more details about our choice of a phase response function and the obtained results in terms of phase response curves, we refer the reader to our response to comment 2.3.

Generally speaking, we chose to use Kuramoto oscillators as it is the simplest model that can provide an oscillatory input to another system while including a phase reset mechanism. This set of oscillators was used to replace the fixed sinusoidal wave that represented theta inputs in previous models (Onslow et al., 2014; Aussel et al., 2018; Segneri et al., 2020). Kuramoto oscillators are a well-established model of synchronization in various fields of physics. They have also been used in neuroscience to model the phase reset of collective rhythms (Levnajić et al. 2010), and the effects of DBS on the basal ganglia network in Parkinson’s disease (Tass et al. 2003, Ebert et al. 2014, Weerasinghe et al. 2019).

More detailed models of the medial septum exist in the literature (e.g., Wang et al. 2002, Hajós et al. 2004) and model the GABAergic effects of the septal projections onto the hippocampal formation. However, it is not trivial to infer the connectivity parameters and the degree of innervation between the hippocampus and the medial septum. Furthermore, the claims made in our study do not necessarily depend on the nature of the projections between the two areas. Therefore, we decided to represent the medial septum in a conceptual way and focus mostly on the effects of these projections rather than replicating them in detail.

Aussel, Amélie, Laure Buhry, Louise Tyvaert, and Radu Ranta. “A Detailed Anatomical and Mathematical Model of the Hippocampal Formation for the Generation of Sharp-Wave Ripples and Theta-Nested Gamma Oscillations.” Journal of Computational Neuroscience 45, no. 3 (December 2018): 207–21. https://doi.org/10.1007/s10827-018-0704-x.

Ebert, Martin, Christian Hauptmann, and Peter A. Tass. “Coordinated Reset Stimulation in a Large-Scale Model of the STN-GPe Circuit.” Frontiers in Computational Neuroscience 8 (2014): 154. https://doi.org/10.3389/fncom.2014.00154.

Hajós, M., W.E. Hoffmann, G. Orbán, T. Kiss, and P. Érdi. “Modulation of Septo-Hippocampal θ Activity by GABAA Receptors: An Experimental and Computational Approach.” Neuroscience 126, no. 3 (January 2004): 599–610. https://doi.org/10.1016/j.neuroscience.2004.03.043.

Levnajić, Zoran, and Arkady Pikovsky. “Phase Resetting of Collective Rhythm in Ensembles of Oscillators.” Physical Review E 82, no. 5 (November 3, 2010): 056202. https://doi.org/10.1103/PhysRevE.82.056202.

Onslow, Angela C. E., Matthew W. Jones, and Rafal Bogacz. “A Canonical Circuit for Generating PhaseAmplitude Coupling.” Edited by Adriano B. L. Tort. PLoS ONE 9, no. 8 (August 19, 2014): e102591. https://doi.org/10.1371/journal.pone.0102591.

Segneri, Marco, Hongjie Bi, Simona Olmi, and Alessandro Torcini. “Theta-Nested Gamma Oscillations in Next Generation Neural Mass Models.” Frontiers in Computational Neuroscience 14 (2020). https://doi.org/10.3389/fncom.2020.00047.Tass, Peter A. “A Model of Desynchronizing Deep Brain Stimulation with a Demand-Controlled Coordinated Reset of Neural Subpopulations.” Biological Cybernetics 89, no. 2 (August 1, 2003): 81–88. https://doi.org/10.1007/s00422-003-0425-7.

Wang, Xiao-Jing. “Pacemaker Neurons for the Theta Rhythm and Their Synchronization in the Septohippocampal Reciprocal Loop.” Journal of Neurophysiology 87, no. 2 (February 1, 2002): 889–900. https://doi.org/10.1152/jn.00135.2001.

Weerasinghe, Gihan, Benoit Duchet, Hayriye Cagnan, Peter Brown, Christian Bick, and Rafal Bogacz. “Predicting the Effects of Deep Brain Stimulation Using a Reduced Coupled Oscillator Model.” PLoS Computational Biology 15, no. 8 (August 8, 2019): e1006575. https://doi.org/10.1371/journal.pcbi.1006575.

No competing interests declared.

Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing.

Software, Formal analysis, Supervision, Investigation, Methodology, Writing – original draft, Writing – review and editing.

Conceptualization, Formal analysis, Supervision, Investigation, Methodology, Writing – original draft, Writing – review and editing.

Conceptualization, Formal analysis, Supervision, Funding acquisition, Investigation, Methodology, Writing – original draft, Writing – review and editing.
==== Refs
References

Abbaspoor S Hussin AT Hoffman KL 2023 Theta- and gamma-band oscillatory uncoupling in the macaque hippocampus eLife 12 e86548 10.7554/eLife.86548 37139864
Adesnik H Scanziani M 2010 Lateral competition for cortical space by layer-specific horizontal circuits Nature 464 1155 1160 10.1038/nature08935 20414303
Akam T Oren I Mantoan L Ferenczi E Kullmann DM 2012 Oscillatory dynamics in the hippocampus support dentate gyrus–CA3 coupling Nature Neuroscience 15 763 768 10.1038/nn.3081 22466505
Ashida G Nogueira W 2018 Spike-Conducting Integrate-and-Fire Model eNeuro 5 ENEURO.0112-18.2018 10.1523/ENEURO.0112-18.2018 30225348
Asllani M Expert P Carletti T 2018 A minimally invasive neurostimulation method for controlling abnormal synchronisation in the neuronal activity PLOS Computational Biology 14 e1006296 10.1371/journal.pcbi.1006296 30024878
Aussel A Buhry L Tyvaert L Ranta R 2018 A detailed anatomical and mathematical model of the hippocampal formation for the generation of sharp-wave ripples and theta-nested gamma oscillations Journal of Computational Neuroscience 45 207 221 10.1007/s10827-018-0704-x 30382451
Aussel A Ranta R Aron O Colnat-Coulbois S Maillard L Buhry L 2022 Cell to network computational model of the epileptic human hippocampus suggests specific roles of network and channel dysfunctions in the ictal and interictal oscillations Journal of Computational Neuroscience 50 519 535 10.1007/s10827-022-00829-5 35971033
Axmacher N Mormann F Fernández G Elger CE Fell J 2006 Memory formation by neuronal synchronization Brain Research Reviews 52 170 182 10.1016/j.brainresrev.2006.01.007 16545463
Axmacher N Henseler MM Jensen O Weinreich I Elger CE Fell J 2010 Cross-frequency coupling supports multi-item working memory in the human hippocampus PNAS 107 3228 3233 10.1073/pnas.0911531107 20133762
Basu I Crocker B Farnes K Robertson MM Paulk AC Vallejo DI Dougherty DD Cash SS Eskandar EN Kramer MM Widge AS 2018 A neural mass model to predict electrical stimulation evoked responses in human and non-human primate brain Journal of Neural Engineering 15 066012 10.1088/1741-2552/aae136 30211694
Bender F Gorbati M Cadavieco MC Denisova N Gao X Holman C Korotkova T Ponomarenko A 2015 Theta oscillations regulate the speed of locomotion via a hippocampus to lateral septum pathway Nature Communications 6 8521 10.1038/ncomms9521 26455912
Bezaire MJ Raikov I Burk K Vyas D Soltesz I 2016 Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent CA1 circuit eLife 5 e18566 10.7554/eLife.18566 28009257
Bingham CS Loizos K Yu GJ Gilbert A Bouteiller JMC Song D Lazzi G Berger TW 2018 Model-Based analysis of electrode placement and pulse amplitude for hippocampal stimulation IEEE Transactions on Bio-Medical Engineering 65 2278 2289 10.1109/TBME.2018.2791860 29993519
Borges RR Borges FS Lameu EL Batista AM Iarosz KC Caldas IL Antonopoulos CG Baptista MS 2017 Spike timing-dependent plasticity induces non-trivial topology in the brain Neural Networks 88 58 64 10.1016/j.neunet.2017.01.010 28189840
Breakspear M Heitmann S Daffertshofer A 2010 Generative models of cortical oscillations: neurobiological implications of the kuramoto model Frontiers in Human Neuroscience 4 190 10.3389/fnhum.2010.00190 21151358
Buño W Garcia-Sanchez JL Garcia-Austt E 1978 Reset of hippocampal rhythmical activities by afferent stimulation Brain Research Bulletin 3 21 28 10.1016/0361-9230(78)90057-6 630419
Buzsáki G 2002 Theta oscillations in the hippocampus Neuron 33 325 340 10.1016/s0896-6273(02)00586-x 11832222
Buzsáki G Buhl DL Harris KD Csicsvari J Czéh B Morozov A 2003 Hippocampal network patterns of activity in the mouse Neuroscience 116 201 211 10.1016/s0306-4522(02)00669-3 12535953
Canolty RT Knight RT 2010 The functional role of cross-frequency coupling Trends in Cognitive Sciences 14 506 515 10.1016/j.tics.2010.09.001 20932795
Capogrosso M Wenger N Raspopovic S Musienko P Beauparlant J Bassi Luciani L Courtine G Micera S 2013 A computational model for epidural electrical stimulation of spinal sensorimotor circuits The Journal of Neuroscience 33 19326 19340 10.1523/JNEUROSCI.1688-13.2013 24305828
Cardin JA Carlén M Meletis K Knoblich U Zhang F Deisseroth K Tsai LH Moore CI 2009 Driving fast-spiking cells induces gamma rhythm and controls sensory responses Nature 459 663 667 10.1038/nature08002 19396156
Chatzikalymniou AP Gumus M Skinner FK 2021 Linking minimal and detailed models of CA1 microcircuits reveals how theta rhythms emerge and their frequencies controlled Hippocampus 31 982 1002 10.1002/hipo.23364 34086375
Colgin LL 2013 Mechanisms and functions of theta rhythms Annual Review of Neuroscience 36 295 312 10.1146/annurev-neuro-062012-170330 23724998
Combrisson E Nest T Brovelli A Ince RAA Soto JLP Guillot A Jerbi K 2020 Tensorpac: An open-source Python toolbox for tensor-based phase-amplitude coupling measurement in electrophysiological brain signals PLOS Computational Biology 16 e1008302 10.1371/journal.pcbi.1008302 33119593
Cutsuridis V Cobb S Graham BP 2010 Encoding and retrieval in a model of the hippocampal CA1 microcircuit Hippocampus 20 423 446 10.1002/hipo.20661 19489002
de Almeida L Idiart M Lisman JE 2007 Memory retrieval time and memory capacity of the CA3 network: role of gamma frequency oscillations Learning & Memory 14 795 806 10.1101/lm.730207 18007022
Dragoi G Buzsáki G 2006 Temporal encoding of place sequences by hippocampal cell assemblies Neuron 50 145 157 10.1016/j.neuron.2006.02.023 16600862
Ebert M Hauptmann C Tass PA 2014 Coordinated reset stimulation in a large-scale model of the STN-GPe circuit Frontiers in Computational Neuroscience 8 154 10.3389/fncom.2014.00154 25505882
Egorov AV Hamam BN Fransén E Hasselmo ME Alonso AA 2002 Graded persistent activity in entorhinal cortex neurons Nature 420 173 178 10.1038/nature01171 12432392
Fell J Axmacher N 2011 The role of phase synchronization in memory processes Nature Reviews. Neuroscience 12 105 118 10.1038/nrn2979 21248789
Fransén E Tahvildari B Egorov AV Hasselmo ME Alonso AA 2006 Mechanism of graded persistent cellular activity of entorhinal cortex layer V neurons Neuron 49 735 746 10.1016/j.neuron.2006.01.036 16504948
Gerstner W Kistler WM Naud R Paninski L 2014 Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition USA Cambridge University Press 10.1017/CBO9781107447615
Giovannini F Knauer B Yoshida M Buhry L 2017 The CAN-In network: A biologically inspired model for self-sustained theta oscillations and memory maintenance in the hippocampus Hippocampus 27 450 463 10.1002/hipo.22704 28052448
Goyal A Miller J Watrous AJ Lee SA Coffey T Sperling MR Sharan A Worrell G Berry B Lega B Jobst BC Davis KA Inman C Sheth SA Wanda PA Ezzyat Y Das SR Stein J Gorniak R Jacobs J 2018 Electrical stimulation in hippocampus and entorhinal cortex impairs spatial and temporal memory The Journal of Neuroscience 38 4471 4481 10.1523/JNEUROSCI.3049-17.2018
Gupta A Vardalakis N Wagner FB 2023 Neuroprosthetics: from sensorimotor to cognitive disorders Communications Biology 6 14 10.1038/s42003-022-04390-w 36609559
Hajós M Hoffmann WE Orbán G Kiss T Erdi P 2004 Modulation of septo-hippocampal Theta activity by GABAA receptors: an experimental and computational approach Neuroscience 126 599 610 10.1016/j.neuroscience.2004.03.043 15183510
Hampel H Mesulam MM Cuello AC Farlow MR Giacobini E Grossberg GT Khachaturian AS Vergallo A Cavedo E Snyder PJ Khachaturian ZS 2018 The cholinergic system in the pathophysiology and treatment of Alzheimer’s disease Brain 141 1917 1933 10.1093/brain/awy132 29850777
Hangya B Borhegyi Z Szilágyi N Freund TF Varga V 2009 GABAergic neurons of the medial septum lead the hippocampal network during theta activity The Journal of Neuroscience 29 8094 8102 10.1523/JNEUROSCI.5665-08.2009 19553449
Hasselmo ME Bodelón C Wyble BP 2002 A proposed function for hippocampal theta rhythm: separate phases of encoding and retrieval enhance reversal of prior learning Neural Computation 14 793 817 10.1162/089976602317318965 11936962
Hendrickson PJ Yu GJ Song D Berger TW 2016 A million-plus neuron model of the hippocampal dentate gyrus: critical role for topography in determining spatiotemporal network dynamics IEEE Transactions on Bio-Medical Engineering 63 199 209 10.1109/TBME.2015.2445771 26087482
Herman PA Lundqvist M Lansner A 2013 Nested theta to gamma oscillations and precise spatiotemporal firing during memory retrieval in a simulated attractor network Brain Research 1536 68 87 10.1016/j.brainres.2013.08.002 23939226
Hirase H Leinekugel X Czurkó A Csicsvari J Buzsáki G 2001 Firing rates of hippocampal neurons are preserved during subsequent sleep episodes and modified by novel awake experience PNAS 98 9386 9390 10.1073/pnas.161274398 11470910
Hodgkin AL Huxley AF 1952 A quantitative description of membrane current and its application to conduction and excitation in nerve The Journal of Physiology 117 500 544 10.1113/jphysiol.1952.sp004764 12991237
Hülsemann MJ Naumann E Rasch B 2019 Quantification of phase-amplitude coupling in neuronal oscillations: comparison of phase-locking value, mean vector length, modulation index, and generalized-linear-modeling-cross-frequency-coupling Frontiers in Neuroscience 13 573 10.3389/fnins.2019.00573 31275096
Hummos A Nair SS 2017 An integrative model of the intrinsic hippocampal theta rhythm PLOS ONE 12 e0182648 10.1371/journal.pone.0182648 28787026
Hyafil A Giraud AL Fontolan L Gutkin B 2015 Neural cross-frequency coupling: connecting architectures, mechanisms, and functions Trends in Neurosciences 38 725 740 10.1016/j.tins.2015.09.001 26549886
Iaccarino HF Singer AC Martorell AJ Rudenko A Gao F Gillingham TZ Mathys H Seo J Kritskiy O Abdurrob F Adaikkan C Canter RG Rueda R Brown EN Boyden ES Tsai LH 2016 Gamma frequency entrainment attenuates amyloid load and modifies microglia Nature 540 230 235 10.1038/nature20587 27929004
Jackson J Dickson CT Bland BH 2008 Median raphe stimulation disrupts hippocampal theta via rapid inhibition and state-dependent phase reset of theta-related neural circuitry Journal of Neurophysiology 99 3009 3026 10.1152/jn.00065.2008 18436639
Jacobs J Miller J Lee SA Coffey T Watrous AJ Sperling MR Sharan A Worrell G Berry B Lega B Jobst BC Davis K Gross RE Sheth SA Ezzyat Y Das SR Stein J Gorniak R Kahana MJ Rizzuto DS 2016 Direct electrical stimulation of the human entorhinal region and hippocampus impairs memory Neuron 92 983 990 10.1016/j.neuron.2016.10.062 27930911
Jensen O Colgin LL 2007 Cross-frequency coupling between neuronal oscillations Trends in Cognitive Sciences 11 267 269 10.1016/j.tics.2007.05.003 17548233
Joucla S Yvert B 2012 Modeling extracellular electrical neural stimulation: From basic understanding to MEA-based applications Journal of Physiology-Paris 106 146 158 10.1016/j.jphysparis.2011.10.003 22036892
Jun S Kim JS Chung CK 2019 Direct stimulation of human hippocampus during verbal associative encoding enhances subsequent memory recollection Frontiers in Human Neuroscience 13 23 10.3389/fnhum.2019.00023 30804768
Jutras MJ Fries P Buffalo EA 2013 Oscillatory activity in the monkey hippocampus during visual exploration and memory formation PNAS 110 13144 13149 10.1073/pnas.1302351110 23878251
Kipping D Nogueira W 2022 A computational model of a single auditory nerve fiber for electric-acoustic stimulation Journal of the Association for Research in Otolaryngology 23 835 858 10.1007/s10162-022-00870-2 36333573
Kitchigina VF 2018 Alterations of coherent theta and gamma network oscillations as an early biomarker of temporal lobe epilepsy and alzheimer’s disease Frontiers in Integrative Neuroscience 12 36 10.3389/fnint.2018.00036 30210311
Kocsis B Bragin A Buzsáki G 1999 Interdependence of multiple theta generators in the hippocampus: a partial coherence analysis The Journal of Neuroscience 19 6200 6212 10.1523/JNEUROSCI.19-14-06200.1999 10407056
Kosenko A Kang S Smith IM Greene DL Langeberg LK Scott JD Hoshi N 2012 Coordinated signal integration at the M-type potassium channel upon muscarinic stimulation The EMBO Journal 31 3147 3156 10.1038/emboj.2012.156 22643219
Kota S Rugg MD Lega BC 2020 Hippocampal theta oscillations support successful associative memory formation The Journal of Neuroscience 40 9507 9518 10.1523/JNEUROSCI.0767-20.2020 33158958
Kuramoto Y 1984 Chemical Oscillations, Waves, and Turbulence Springer 10.1007/978-3-642-69689-3
Kwag J Paulsen O 2009 The timing of external input controls the sign of plasticity at local synapses Nature Neuroscience 12 1219 1221 10.1038/nn.2388 19734896
Kwag J Jang HJ Kim M Lee S 2014 M-type potassium conductance controls the emergence of neural phase codes: a combined experimental and neuron modelling study Journal of the Royal Society, Interface 11 99 10.1098/rsif.2014.0604 25100320
Lacruz ME Valentín A Seoane JJG Morris RG Selway RP Alarcón G 2010 Single pulse electrical stimulation of the hippocampus is sufficient to impair human episodic memory Neuroscience 170 623 632 10.1016/j.neuroscience.2010.06.042 20643192
Larson J Munkácsy E 2015 Theta-burst LTP Brain Research 1621 38 50 10.1016/j.brainres.2014.10.034 25452022
Lega BC Jacobs J Kahana M 2012 Human hippocampal theta oscillations and the formation of episodic memories Hippocampus 22 748 761 10.1002/hipo.20937 21538660
Lega B Burke J Jacobs J Kahana MJ 2016 Slow-theta-to-gamma phase–amplitude coupling in human hippocampus supports the formation of new episodic memories Cerebral Cortex 26 268 278 10.1093/cercor/bhu232 25316340
Lengyel M Kwag J Paulsen O Dayan P 2005 Matching storage and recall: hippocampal spike timing-dependent plasticity and phase response curves Nature Neuroscience 8 1677 1683 10.1038/nn1561 16261136
Levnajić Z Pikovsky A 2010 Phase resetting of collective rhythm in ensembles of oscillators arXiv http://arxiv.org/abs/1007.4097
Lin J-J Rugg MD Das S Stein J Rizzuto DS Kahana MJ Lega BC 2017 Theta band power increases in the posterior hippocampus predict successful episodic memory encoding in humans Hippocampus 27 1040 1053 10.1002/hipo.22751 28608960
Lisman JE Talamini LM Raffone A 2005 Recall of memory sequences by interaction of the dentate and CA3: A revised model of the phase precession Neural Networks 18 1191 1201 10.1016/j.neunet.2005.08.008 16233972
Lisman JE Jensen O 2013 The theta-gamma neural code Neuron 77 1002 1016 10.1016/j.neuron.2013.03.007 23522038
Lozano AM Fosdick L Chakravarty MM Leoutsakos JM Munro C Oh E Drake KE Lyman CH Rosenberg PB Anderson WS Tang-Wai DF Pendergrass JC Salloway S Asaad WF Ponce FA Burke A Sabbagh M Wolk DA Baltuch G Okun MS Foote KD McAndrews MP Giacobbe P Targum SD Lyketsos CG Smith GS 2016 A phase ii study of fornix deep brain stimulation in mild alzheimer’s disease Journal of Alzheimer’s Disease 54 777 787 10.3233/JAD-160017
Lu Y Truccolo W Wagner FB Vargas-Irwin CE Ozden I Zimmermann JB May T Agha NS Wang J Nurmikko AV 2015 Optogenetically induced spatiotemporal gamma oscillations and neuronal spiking activity in primate motor cortex Journal of Neurophysiology 113 3574 3587 10.1152/jn.00792.2014 25761956
Lundqvist M Rehn M Djurfeldt M Lansner A 2006 Attractor dynamics in a modular network model of neocortex Network 17 253 276 10.1080/09548980600774619 17162614
Lurie SM Kragel JE Schuele SU Voss JL 2022 Human hippocampal responses to network intracranial stimulation vary with theta phase eLife 11 e78395 10.7554/eLife.78395
Malkov A Shevkova L Latyshkova A Kitchigina V 2022 Theta and gamma hippocampal-neocortical oscillations during the episodic-like memory test: Impairment in epileptogenic rats Experimental Neurology 354 114110 10.1016/j.expneurol.2022.114110 35551900
Manseau F Goutagny R Danik M Williams S 2008 The hippocamposeptal pathway generates rhythmic firing of gabaergic neurons in the medial septum and diagonal bands: an investigation using a complete septohippocampal preparation in vitro The Journal of Neuroscience 28 4096 4107 10.1523/JNEUROSCI.0247-08.2008 18400909
McCartney H Johnson AD Weil ZM Givens B 2004 Theta reset produces optimal conditions for long-term potentiation Hippocampus 14 684 687 10.1002/hipo.20019 15318327
McIntyre CC Grill WM Sherman DL Thakor NV 2004 Cellular effects of deep brain stimulation: model-based analysis of activation and inhibition Journal of Neurophysiology 91 1457 1469 10.1152/jn.00989.2003
Mina F Benquet P Pasnicu A Biraben A Wendling F 2013 Modulation of epileptic activity by deep brain stimulation: a model-based study of frequency-dependent effects Frontiers in Computational Neuroscience 7 7 10.3389/fncom.2013.00094
Mizuseki K Sirota A Pastalkova E Buzsáki G 2009 Theta oscillations provide temporal windows for local circuit computation in the entorhinal-hippocampal loop Neuron 64 267 280 10.1016/j.neuron.2009.08.037 19874793
Mormann F Fell J Axmacher N Weber B Lehnertz K Elger CE Fernández G 2005 Phase/amplitude reset and theta-gamma interaction in the human medial temporal lobe during a continuous word recognition memory task Hippocampus 15 890 900 10.1002/hipo.20117
Müller C Remy S 2018 Septo-hippocampal interaction Cell and Tissue Research 373 565 575 10.1007/s00441-017-2745-2 29250747
Mysin IE Kitchigina VF Kazanovich YB 2019 Phase relations of theta oscillations in a computer model of the hippocampal CA1 field: Key role of Schaffer collaterals Neural Networks 116 119 138 10.1016/j.neunet.2019.04.004
Nelson AR Kolasa K McMahon LL 2014 Noradrenergic sympathetic sprouting and cholinergic reinnervation maintains non-amyloidogenic processing of AβPP Journal of Alzheimer’s Disease 38 867 879 10.3233/JAD-130608
Neymotin SA Lazarewicz MT Sherif M Contreras D Finkel LH Lytton WW 2011 Ketamine disrupts theta modulation of gamma in a computer model of hippocampus The Journal of Neuroscience 31 11733 11743 10.1523/JNEUROSCI.0501-11.2011 21832203
Nuñez A Buño W 2021 The theta rhythm of the hippocampus: from neuronal and circuit mechanisms to behavior Frontiers in Cellular Neuroscience 15 649262 10.3389/fncel.2021.649262 33746716
Onslow ACE Jones MW Bogacz R 2014 A canonical circuit for generating phase-amplitude coupling PLOS ONE 9 e102591 10.1371/journal.pone.0102591
Pirini M Rocchi L Sensi M Chiari L 2009 A computational modelling approach to investigate different targets in deep brain stimulation for Parkinson’s disease Journal of Computational Neuroscience 26 91 107 10.1007/s10827-008-0100-z 18553128
Ponzi A Dura-Bernal S Migliore M 2023 Theta-gamma phase amplitude coupling in a hippocampal CA1 microcircuit PLOS Computational Biology 19 e1010942 10.1371/journal.pcbi.1010942 36952558
Rattay F 1986 Analysis of models for external stimulation of axons IEEE Transactions on Bio-Medical Engineering 33 974 977 10.1109/TBME.1986.325670
Rattay F Minassian K Dimitrijevic MR 2000 Epidural electrical stimulation of posterior structures of the human lumbosacral cord: 2. quantitative analysis by computer modeling Spinal Cord 38 473 489 10.1038/sj.sc.3101039 10962608
Rattay F Resatz S Lutter P Minassian K Jilge B Dimitrijevic MR 2003 Mechanisms of electrical stimulation with neural prostheses: mechanisms in electrical stimulation Neuromodulation: Technology at the Neural Interface 6 42 56 10.1046/j.1525-1403.2003.03006.x
Rattay F Bassereh H Stiennon I 2018 Compartment models for the electrical stimulation of retinal bipolar cells PLOS ONE 13 e0209123 10.1371/journal.pone.0209123 30557410
Rizzuto DS Madsen JR Bromfield EB Schulze-Bonhage A Seelig D Aschenbrenner-Scheibe R Kahana MJ 2003 Reset of human neocortical oscillations during a working memory task PNAS 100 7931 7936 10.1073/pnas.0732061100 12792019
Rougier NP 2018 A density-driven method for the placement of biological cells over two-dimensional manifolds Frontiers in Neuroinformatics 12 12 10.3389/fninf.2018.00012 29615887
Rubin JE Terman D 2004 High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model Journal of Computational Neuroscience 16 211 235 10.1023/B:JCNS.0000025686.47117.67 15114047
Segneri M Bi H Olmi S Torcini A 2020 Theta-nested gamma oscillations in next generation neural mass models Frontiers in Computational Neuroscience 14 47 10.3389/fncom.2020.00047 32547379
Solomon EA Sperling MR Sharan AD Wanda PA Levy DF Lyalenko A Pedisich I Rizzuto DS Kahana MJ 2021 Theta-burst stimulation entrains frequency-specific oscillatory responses Brain Stimulation 14 1271 1284 10.1016/j.brs.2021.08.014 34428553
Stimberg M Brette R Goodman DF 2019 Brian 2, an intuitive and efficient neural simulator eLife 8 e47314 10.7554/eLife.47314 31429824
Sun J Kapur J 2012 M‐type potassium channels modulate Schaffer collateral–CA1 glutamatergic synaptic transmission The Journal of Physiology 590 3953 3964 10.1113/jphysiol.2012.235820 22674722
Suthana N Haneef Z Stern J Mukamel R Behnke E Knowlton B Fried I 2012 Memory enhancement and deep-brain stimulation of the entorhinal area The New England Journal of Medicine 366 502 510 10.1056/NEJMoa1107212 22316444
Suthana N Fried I 2014 Deep brain stimulation for enhancement of learning and memory NeuroImage 85 Pt 3 996 1002 10.1016/j.neuroimage.2013.07.066 23921099
Takeuchi Y Nagy AJ Barcsai L Li Q Ohsawa M Mizuseki K Berényi A 2021 The medial septum as a potential target for treating brain disorders associated with oscillopathies Frontiers in Neural Circuits 15 701080 10.3389/fncir.2021.701080 34305537
Tass PA 2003 A model of desynchronizing deep brain stimulation with A demand-controlled coordinated reset of neural subpopulations Biological Cybernetics 89 81 88 10.1007/s00422-003-0425-7 12905037
Titiz AS Hill MRH Mankin EA M Aghajan Z Eliashiv D Tchemodanov N Maoz U Stern J Tran ME Schuette P Behnke E Suthana NA Fried I 2017 Theta-burst microstimulation in the human entorhinal area improves memory specificity eLife 6 e29515 10.7554/eLife.29515 29063831
Torben-Nielsen B Uusisaari M Stiefel KM 2010 A comparison of methods to determine neuronal phase-response curves Frontiers in Neuroinformatics 4 6 10.3389/fninf.2010.00006 20431724
Tort ABL Rotstein HG Dugladze T Gloveli T Kopell NJ 2007 On the formation of gamma-coherent cell assemblies by oriens lacunosum-moleculare interneurons in the hippocampus PNAS 104 13490 13495 10.1073/pnas.0705708104 17679692
Tort ABL Kramer MA Thorn C Gibson DJ Kubota Y Graybiel AM Kopell NJ 2008 Dynamic cross-frequency couplings of local field potential oscillations in rat striatum and hippocampus during performance of a T-maze task PNAS 105 20517 20522 10.1073/pnas.0810524105 19074268
Tort ABL Komorowski RW Manns JR Kopell NJ Eichenbaum H 2009 Theta–gamma coupling increases during the learning of item–context associations PNAS 106 20942 20947 10.1073/pnas.0911331106 19934062
Tort ABL Komorowski R Eichenbaum H Kopell N 2010 Measuring phase-amplitude coupling between neuronal oscillations of different frequencies Journal of Neurophysiology 104 1195 1210 10.1152/jn.00106.2010 20463205
Tóth K Borhegyi Z Freund TF 1993 Postsynaptic targets of GABAergic hippocampal neurons in the medial septum-diagonal band of broca complex The Journal of Neuroscience 13 3712 3724 10.1523/JNEUROSCI.13-09-03712.1993 7690065
Traub RD Jefferys JGR Whittington MA 1997 Simulation of gamma rhythms in networks of interneurons and pyramidal cells Journal of Computational Neuroscience 4 141 150 10.1023/a:1008839312043 9154520
Traub RD Bibbig A LeBeau FEN Buhl EH Whittington MA 2004 Cellular mechanisms of neuronal population oscillations in the hippocampus in vitro Annual Review of Neuroscience 27 247 278 10.1146/annurev.neuro.27.070203.144303
Tsanov M Manahan-Vaughan D 2009 Long-term plasticity is proportional to theta-activity PLOS ONE 4 e5850 10.1371/journal.pone.0005850
Unal G Joshi A Viney TJ Kis V Somogyi P 2015 Synaptic targets of medial septal projections in the hippocampus and extrahippocampal cortices of the mouse The Journal of Neuroscience 35 15812 15826 10.1523/JNEUROSCI.2639-15.2015 26631464
Vardalakis N 2023 Memstim-HH swh:1:rev:563f808f6c4f40630f5b8876cc0b440cdf4159e8 Software Heritage https://archive.softwareheritage.org/swh:1:dir:98abfedbd225e7c63d89101cb0605b7c4a785aa4;origin=https://github.com/NikVard/memstim-hh;visit=swh:1:snp:c2e27ce82d59e51c3aa0d7aaece24f9f6f9a2889;anchor=swh:1:rev:563f808f6c4f40630f5b8876cc0b440cdf4159e8
Varga V Hangya B Kránitz K Ludányi A Zemankovics R Katona I Shigemoto R Freund TF Borhegyi Z 2008 The presence of pacemaker HCN channels identifies theta rhythmic GABAergic neurons in the medial septum The Journal of Physiology 586 3893 3915 10.1113/jphysiol.2008.155242 18565991
Wang XJ 2002 Pacemaker neurons for the theta rhythm and their synchronization in the septohippocampal reciprocal loop Journal of Neurophysiology 87 889 900 10.1152/jn.00135.2001 11826054
Weerasinghe G Duchet B Cagnan H Brown P Bick C Bogacz R 2019 Predicting the effects of deep brain stimulation using a reduced coupled oscillator model PLOS Computational Biology 15 e1006575 10.1371/journal.pcbi.1006575 31393880
Williams JM Givens B 2003 Stimulation-induced reset of hippocampal theta in the freely performing rat Hippocampus 13 109 116 10.1002/hipo.10082 12625462
