Introduction

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., 2021; Abbaspoor et al., 2022). 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). 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 inconsistentwith theta phase reset. To circumventthis 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 resetwas 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). The 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 S1).

Figure 1.

As in the original model (Aussel et al., 2018), the input theta rhythm drove the network to produce spiking activity preferentially around the peak of theta in all excitatory and inhibitory populations (Figure 2A and Figure S2). 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 1A). 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 2B). To quantify PAC between theta and gamma oscillations, we used the modulation index (MI) (Tort et al., 2008, 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 S3). 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 S2B). 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 2C and Figure S2B).

Figure 2.

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 S4), 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 2A). 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 2D). 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.

Overall, we successfully developed a new model of the hippocampal formation able to exhibit 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 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 3A, 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 thatthe 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 3B). Moreover, the time interval separating the two bursts likely resulted from the interplay between the depolarizing CAN current and hyperpolarizing M current (Figure 3A).

Figure 3.

In the presence of theta input, the effects of single-pulse stimulation depended both on theta drive amplitude and stimulation amplitude, highlighting different regimes of network activity (Figure 4 and Figure S5, Figure S6, Figure S7). For low theta input, theta-nested gamma oscillations were initially absent and could not be induced by stimulation (Figure 4A, B). At most, the stimulation could only elicit a few bursts of spiking activity that faded away after approximately 500 ms (Figure 4B). For increasing theta input, the network switched to a bistable 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 4C). This regime existed for a range of theta inputs and stimulation amplitudes, as characterized by the post-stimulation theta power, gamma power, and theta-gamma PAC (Figure 4D). For the highest theta input, the network became monostable, spontaneously generating theta-nested gamma oscillations even when initialized at a state with no spiking activity (Figure 4D, last stimulation amplitude). Additionally, we investigated the effects of various synchronization levels of the Kuramoto oscillators on the production of theta-gamma oscillations and their response to single-pulse stimulation. The results showed that the production of theta-nested gamma oscillations depended mostly on the theta input amplitude and less on the synchronization level (Figure S8). Therefore, we decided to use the amplitude of the theta input as the main variable controlling whether the neuronal activity was within a physiological or pathological state.

Figure 4.

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 associated with a bistable regime that allowed self-sustained oscillations (Figure 4, point C), and one with a weak theta input characterized by the complete absence of self-sustained oscillations (Figure 4, point B). 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 2).

In the case of a bistable network (moderate theta input) with phase reset, delivering a pulse at either the peak or trough of theta could induce theta-nested gamma oscillations (Figure 5A, C). By contrast, in the absence of phase reset, only stimulation delivered at the peak of theta was able to induce such oscillations (Figure 5B, 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 5E).

Figure 5.

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 6). 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 6A, C). There was a slight difference between the responses caused by the first pulse, as trough stimulation missed one cycle before producing network-wide oscillations. 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, only peak-targeted stimulation induced theta-gamma oscillations in all network areas (Figure 5B, D). Comparing these simulations based on theta power, gamma power, and theta-gamma PAC within CA1 (Figure 6E) showed similar differences as single-pulse stimulation, although less striking because CA1 was driven directly by the stimulation. However, these differences were even more pronounced in areas other than CA1, as can be seen by the complete absence of oscillations for trough stimulation without phase reset (Figure 6D).

Figure 6.

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 consists in the bistable nature of the model for certain values of the theta amplitude (Figure 4). In the bistable regime, a single stimulation pulse could switch the network behavior from a state with no spiking activity to one with theta-nested gamma oscillations. Whether such bistability exists 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 pulses with respect to the phase of the ongoing theta rhythm matters differently depending on the presence or absence of phase reset (Figure 5 and Figure 6). 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). 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 6).

It should also be noted that the reset gain in our simulations (Equation 1, gain G_reset) 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 4). 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 more mitigated (Figure S8). 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). 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, previous models have also proposed that direct septal inputs are not essential for theta generation (Wang, 2002; Mysin et al., 2019), but play an important role in phase synchronization of hippocampal neurons. To account for these discrepancies, our model could be extended by considering more realistic connectivity patterns between the medial 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.

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.

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, i.e. 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 underfixed 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. 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 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:

The term ω_i denotes the natural frequency of oscillator i, and is normally distributed around the center frequency f₀ 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 G_resetX(t)Z(θ) 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. G_reset is an arbitrary gain that determines the strength of the CA1 input to the oscillators. Finally, the function Z(θ) describes howthe 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:

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:

where t^s denotes the ordered set of the spike timings, and τ_FR determines the exponential decay of the kernel, which was set to 10 ms to compute population activity (for single neurons, a typical value is 100 ms) (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 ϕ_offse, controls the desired phase offset from the peak:

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.

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, 2018). The temporal evolution of the membrane potential of each neuron is described by a differential equation whose general form reads:

I_L denotes the leakage current. I_channel are currents associated with specific ion channels, namely potassium (I_K), fast sodium (I_Na), and low-threshold calcium (I_Ca) currents, the CAN current (I_CAN) (Giovannini et al., 2017), and the M-type potassium channel current (I_M) responsible for spike adaptation (Kosenko et al., 2012; Sun and Kapur, 2012; Kwag et al., 2014). I_syn 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 I_stim. Excitatory neurons represent pyramidal cells in EC, CA3, and CA1, and granule cells in DG. They were modeled with I_Na, I_κ, and I_Ca, and I_M currents, with the addition of I_CAN for pyramidal cells. Fast-spiking interneurons in all areas were modeled with I_Na and I_K 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:

where g_L is the maximum leakage conductance, A is the area of the single compartment corresponding to the membrane of a neuron, and E_L is the reversal potential of the leakage channel. Channel currents I_K, I_M, I_CAN obey the following set of equations:

where g_K, g_M, g_CAN are the maximum conductances for the respective channel, and n, p, m_CAN are the respective gating variables defined by the following differential equations:

For the potassium and CAN currents, the steady-state values for their corresponding gating variables n_∞ and m_CAN,∞ and their corresponding time constants τ_K and τ_CAN depend on the following functions of the transition rate constants:

The sodium current (I_Na) and calcium current (I_Ca) follow a similar set of equations:

with two gating variables m and h defined by the following differential equations:

The gating variable of I_CAN depends on the calcium concentration within the neuron , given by:

where τ_[Ca]²⁺ = 1s represents the rate of calcium removal from the cell, is the calcium concentration if the calcium channel remains open for a duration of ΔT → ∞, k_u = 10⁴ 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 10 μV (η_E ~ N(0, 10) μV) for excitatory neurons and a mean of 0 μV and standard deviation of 1 μV (η_I ~ N(0, 1) μ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.

Table 3.

Synaptic Models

Inter-neuronal interactions were modeled as instantaneous AMPAand GABA-A synapses usingthe synaptic currents I_synE and I_synI respectively. Synaptic currents were described by the following biexponential differential equations:

where E_I,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 h_I or h_E 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.

Table 5.

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.

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 A_intra defines the maximum probability of connection between two neurons separated by an infinitesimal distance (i.e. for D → 0). The maximum values for A_intra are given in Table 7.

Table 7

Inter-area connectivity followed a similar Gaussian-like distribution (Equation 15), however, only excitatory projections were considered and the distance IX 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.

Tuning input gain and connection strengths: targeted firing rates

We adjusted the input gain (G_θ) and inter-area connection strengths (A_inter) to target an overall oscillatory rhythm f_osc,targ at the driving frequency of the input, i.e. 6 Hz, and a mean firing rate of each excitatory population f_exc,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 f_exc 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 f_osc 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:

Tuning input and connection strengths: detailed procedure

The input gain (G_θ, Equation 2) and inter-area connection strengths (A_inter, 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.

1. 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.

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

3. 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.

4. 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.

5. 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.

6. 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.

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.1 ms and a total simulation duration ranging between 3-10 seconds (depending on the experiment, see results section). The average time for simulating 1 second for the complete model locally was 10 minutes.

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) 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.5 ms). 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 2A), 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 [f_l, f_h]. 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 S3.

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 2C.

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.5 ms 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 2D for varying stimulation amplitudes.

Acknowledgements

Simulations presented in this paperwere 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

Funding

This work was supported by the Bordeaux Neurocampus junior chair program co-funded by the Région Nouvelle-Aquitaine and the University of Bordeaux Initiative of Excellence (IdEx) (F.B.W., Neurocampus junior chair), and the European Union’s Horizon Europe research and innovation program under European Research Council (ERC) Starting Grant agreement N° 101040391 (F.B.W., MEMOPROSTHETICS).

Author contributions

N.V., Conceptualization, Methodology, Software, Formal analysis, Investigation, Visualization (lead), Writing - original draft, Writing - review and editing; A.A., Methodology, Formal analysis, Investigation, Software, Visualization (supporting), Writing - original draft, Writing - review and editing, Supervision; N.P.R., Conceptualization, Methodology, Formal analysis, Investigation, Visualization (supporting), Writing - original draft, Writing - review and editing, Supervision; F.B.W., Conceptualization, Methodology, Formal analysis, Investigation, Visualization (supporting), Writing - original draft, Writing - review and editing, Supervision, Funding acquisition.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability

All model files and analysis scripts are available at the following repositories under the GNU General Public License v3.0 license:

ZENODO https://www.doi.org/10.5281/zenodo.7702041

GITHUB https://github.com/NikVard/memstim-hh

Appendix

Supplementary Figures

Figure S1. (Figure 1 – Figure Supplement 1)

Figure S2. (Figure 2 – Figure Supplement 1)

Figure S3. (Figure 2 – Figure Supplement 2)

Figure S4. (Figure 2 – Figure Supplement 3)

Figure S5. (Figure 4 – Figure Supplement 1)

Figure S6. (Figure 4 – Figure Supplement 2)

Figure S7. (Figure 4 – Figure Supplement 3)

Figure S8. (Figure 4 – Figure Supplement 4)