Introduction

Pavlovian fear conditioning is widely used as a model of associative learning across multiple species (Phelps and LeDoux, 2005), and has been used more generally to study plasticity and memory formation. The major questions, still under investigation, relate to how plasticity is instantiated in defined circuits and how it is regulated by the circuit components (Rumpel et al., 2005; Sah et al., 2008; Johansen et al., 2014; Bocchio et al., 2017; Grewe et al., 2017). The fear conditioning paradigm consists of a neutral stimulus (conditioned stimulus, CS) presented one or more times together with an aversive stimulus (unconditioned stimulus, US), which induces a fear response.

In the basolateral amygdala (BLA), the main site of fear learning in the mammalian brain (Fanselow and LeDoux, 1999; Tovote et al., 2015; Krabbe et al., 2018), local field potential recordings (LFP) show prominent low theta (∼3-6 Hz), high theta (∼6-12 Hz), and gamma (>30 Hz) rhythms (Seidenbecher et al., 2003; Courtin et al., 2014b; Stujenske et al., 2014; Davis et al., 2017). Recent rodent studies show increased low theta (Davis et al., 2017) and gamma (Courtin et al., 2014b) in local field potential (LFP) recordings after successful fear conditioning, whereas modulation of high theta is associated with fear extinction (Davis et al., 2017), a paradigm that aims to suppress the association between CS and fear (Bouton, 2004); the modulation of the power of these rhythms suggests they may be associated with BLA plasticity.

The origin of the rhythms and their potential roles in instantiating the plasticity needed for successful learning is still under investigation (Bocchio et al., 2017). Fear responses are mediated by a network of several brain areas including, but not limited to, the BLA, the prefrontal cortex (mPFC), and the hippocampus (Seidenbecher et al., 2003; Stujenske et al., 2014; Tovote et al., 2015). In experimental studies, successful fear conditioning leads fear-encoding neurons to respond to the CS. In our work, it results in strengthening the connection between neurons encoding CS and neurons encoding fear. We create a biophysically detailed model of the BLA circuit suggesting that interneurons in the BLA can be critically involved in producing the experimentally measured rhythms and are mechanistically important in the instantiation of plasticity during learning. Indeed, we find that the low theta, fast theta, and gamma rhythms in the BLA originating from the BLA interneurons promote the formation of a dedicated fear circuit shaped through rhythmic regulation of depression-dominated spike-timing-dependent plasticity (STDP). (See Discussion about non-depression-dominated plasticity related to fear conditioning.) In this model, if any of the classes of interneurons are removed from the circuit, the rhythms are changed and there is a failure of the plasticity needed for successful learning. We show that fine timing between excitatory projection neurons responding to the CS and US is necessary but not sufficient to produce associative plasticity when the latter can be affected by all the spikes, not just nearest neighbors. The other critical element for plasticity is the interaction of interneurons that creates pauses in excitatory cell activity.

The model reproduces the increase in the low theta after training that was found in the experimental data. This increase was not seen for network instantiations that did not learn; thus, the simulations suggest that the increase in low theta is a biomarker of successful fear conditioning. Furthermore, this low theta signal emanates from the newly formed synapse linking CS- and fear-encoding neurons. These results are also highly relevant to the recent study (Perrenoud and Cardin, 2023; Veit et al., 2023) in which the same types of interneurons used in this paper are shown to allow global selectivity in coherence in gamma oscillations (see Discussion).

Results

Rhythms in the BLA can be produced by interneurons

Brain rhythms are thought to be encoded and propagated largely by interneurons (Whittington et al., 2000). Identified interneurons in BLA include VIP (vasoactive intestinal peptide-expressing), SOM (somatostatin-expressing), and PV (parvalbumin-expressing) (Muller et al., 2006, 2007; Rainnie et al., 2006; Bienvenu et al., 2012; Krabbe et al., 2018, 2019), which can be further subdivided according to their electrophysiological dynamics (Sosulina et al., 2010; Spampanato et al., 2011). In the model, we show that some types of VIP, SOM, and PV can each contribute to the generation of a key rhythm involved in the BLA due to their specific intrinsic currents.

For VIP interneurons, we consider a subtype that responds to a depolarizing step current with bursting or stuttering behavior (Sosulina et al., 2010; Spampanato et al., 2011). This type of behavior can be elicited by a D-type potassium current, a current thought to be found in a similar electrophysiological subtype of VIP interneurons in the cortex (Porter et al., 1998). In our model, VIP interneurons endowed with a D-current respond to depolarizing currents with long bursts, and to hyperpolarizing currents with no action potentials, thus reproducing the electrophysiological properties of type I BLA interneurons in (Sosulina et al., 2010) (Fig. 1A, top): our VIP neurons exhibit gamma (∼52 Hz) bursting activity at low theta frequencies (∼3-6 Hz) in the baseline condition, the condition without any external input from the fear conditioning paradigm (Fig. 1B, top).

Figure 1.

SOM interneurons have been well studied in the hippocampus, especially the O-LM cells (Maccaferri and McBain, 1996; Gillies et al., 2002; Saraga et al., 2003; Rotstein et al., 2005). Although the intrinsic currents of BLA SOM cells are unknown, SOM cells in the hippocampus have a hyperpolarization-activated current, i.e., H-current, and a persistent sodium current, i.e., NaP-current. In our model, with the introduction of these currents with specific conductances (see Materials and Methods for details), the SOM cells mimic the electrophysiologic behavior of type III BLA SOM cells in (Sosulina et al., 2010), showing regular spikes with early spike-frequency adaptation in response to a depolarizing current, and pronounced inward rectification (downward deflection) and outward rectification (upward deflection) upon the initiation and release of a hyperpolarizing current (Fig. 1A, middle). In our baseline model, SOM cells have a natural frequency rhythm in the high theta range (6-12 Hz) (Fig. 1B, middle).

Our model PV interneurons are fast-spiking interneurons (FSIs) with standard action potentials produced by Hodgkin-Huxley-type sodium and potassium conductances. They are silent at baseline condition and show similar behaviors to type IV interneurons (Sosulina et al., 2010) in response to depolarizing and hyperpolarizing currents (Fig. 1A, bottom). At baseline, our PV interneurons are silent (Fig. 1B, bottom). However, when reciprocal connections are present and the excitatory projection neuron receives enough excitation to fire, the PV interneuron forms a PING rhythm (pyramidal-interneuron network gamma) (Whittington et al., 2000) (Fig. 1C); this has been suggested as a possible mechanism for the basis of gamma rhythm generation in the BLA (Feng et al., 2019). The reciprocal interactions between the projection neuron and PV form a gamma rhythm that depends sensitively on the external input to the excitatory projection neuron and the PV’s decay time constant of inhibition (8.3 ms).

Interneurons interact to modulate fear neuron output

Our BLA network consists of interneurons, detailed in the previous section, and excitatory projection neurons (Fig. 2A). In response to the noxious stimulus US, this network produces all the rhythms originating from the interneurons, as shown in Fig. 1. The simplest network we consider is made of one neuron for each cell type. We introduce some heterogeneity in the last two sections of the Results.

Figure 2.

Both the fear-encoding neuron (F), an excitatory projection neuron, and the VIP interneuron are activated by the noxious stimulus US (Krabbe et al., 2019). As shown in Fig. 2A (top, right), VIP disinhibits F by inhibiting both SOM and PV, as suggested in (Krabbe et al., 2019). Since VIP is active during both rest and upon the injection of US, it then modulates F at low theta cycles via SOM and PV. We do not include connections from PV to SOM and VIP, nor connections from SOM to PV and VIP, since those connections have been shown to be significantly weaker than the ones included (Krabbe et al., 2019). Fig. 2A (bottom) shows a typical dynamic of the network before and after the US input onset. In the baseline condition, the VIP interneuron has short gamma bursts nested in low theta rhythm. With US onset, VIP increases its burst duration and the frequency of low theta rhythm. These longer bursts make the SOM cell silent for long periods of each low theta cycle, providing F with windows of disinhibition and contributing to the abrupt increase in activity right after the US onset. Finally, in Fig. 2A, PV lacks any external input and fires only when excited by F. Thanks to their reciprocal interactions, PV forms a PING rhythm with F, as depicted in Fig.1C.

Interneuron rhythms provide the fine timing needed for depression-dominated STDP to make the association between CS and fear

We now introduce another excitatory projection neuron (ECS), as shown in Fig. 2B (top). ECS, unlike F, responds to the neutral stimulus CS, as does PV. By the end of fear conditioning, CS consistently activates the neuron F, thus eliciting the network fear response. This happens because of the formation and strengthening of the synapse from ECS to F by means of synaptic plasticity. We now show how this network, with appropriate connection strengths among neurons, can make the timing of the interneurons confer pre-post timing to ECS and F, which is conducive to spike-timing-dependent plasticity potentiation suggested to be critical for associative aversive learning (Rogan et al., 1997; Nabavi et al., 2014); in particular, we need feedback inhibition (from PV to F) to be stronger than lateral inhibition (from PV to ECS) to promote ECS firing before F. The Hebbian plasticity rule that we use is characterized by a longer time constant of depression than potentiation (and equal maximal amplitudes) and considers the whole history of ECS and F spiking activity (see Materials and Methods and Fig. S1A for more details).

Fig. 2B (bottom) shows an example of the network dynamics with CS present for 2 seconds and US injected after the first second of simulation. ECS is active during the whole 2-second interval. CS also affects PV, which is active most of the time; F is active only in the second half of the dynamics when US is also present. All the rhythms generated by the interneurons are apparent in response to CS and simultaneous CS-US inputs, and they are generated by the same mechanisms as in the fear-only network (Fig. 2A). In contrast to the fear-only network, the projection neurons F and ECS are both modulated by the VIP low theta rhythm. This is because the PV increased activity due to CS tends to silence (with the help of SOM) ECS and F during the silent VIP phase at low theta. During the active VIP phase, however, both ECS and F are active, and the simulations show that ECS fires most of the time slightly before F (see Fig. 2B,C and Fig. S1B); this fine timing needed for potentiation is established by the PING rhythm (see Fig. 1C). By contrast, in the first second of the simulation in Fig. 2B (bottom), SOM and PV prevent plasticity by silencing F, which in the absence of US receives only a weak applied current: thus, without concurrent US and CS, STDP cannot occur between ECS and F. In the next sections, we will explore the role of each interneuron and its associated rhythm in shaping the network dynamics and allowing the association between CS and fear to be instantiated.

With the depression-dominated plasticity rule, all interneuron types are needed to provide potentiation during fear learning

We now show that, in the example used above, only the network endowed with all the interneurons and their associated rhythms leads to overall potentiation of the conductance from ECS to F in the timeframe used to induce the fear learning in experimental work. (see Discussion for other plasticity rules.) In general, experimental work finds successful learning after one or very few presentations of CS and US (lasting 1.5 or 2 seconds) interspersed with CS-only intervals lasting 30-40 seconds (e.g., see (Davis et al., 2017; Krabbe et al., 2019)). The 40-second interval we consider consists of US and CS active during the entire period: since an initial bout of US is known to produce a long-lasting fear response beyond the offset of the US (Hole and Lorens, 1975), we treat the US as present (e.g., for 1.5 or 2 seconds) followed by US activity induced by memory for the remainder of the 40 seconds.

Fig. 3A shows the evolution of the average ECS to F AMPA maximal conductance during fear conditioning across 40 network realizations of the full network, as well as from networks lacking VIP, PV, SOM. The average ECS to F AMPA conductance robustly potentiates only in the network containing all three interneurons (Fig. 3A; see Fig. S3 for exceptions outside the normal physiological range). Furthermore, as shown in Fig. 3B, all the full network realizations are “learners”. We define learners as those realizations whose AMPA conductance from ECS to F is higher than 0.12 mS/cm³ at the end of the 40-second interval, which results in the systematic activation of the fear neuron F following each ECS spike when only CS is presented. This is consistent with the high rate of successful learning in rodent experiments after one pairing of CS and US, despite inter-individual differences among animals (Schafe et al., 2000). Network realizations differ from one another in the initial state of each neuron involved, and all receive independent Gaussian noise.

Figure 3.

We show in Fig. 3C results in which we relax the assumption that US is active over the entire forty seconds by allowing US to be active only during the initial 15 seconds of the fear conditioning paradigm. We find that learning may still occur under this condition. The reason is that, after 15 seconds of CS and US pairing, ECS to F may have potentiated enough that ECS can drive F some of the time (although not all the time). This allows further potentiation to occur in the presence of CS alone. We find that, if the conductance from ECS to F is higher than a threshold value of 0.037 mS/cm³ after the 15 seconds of CS and US presentation, then the network will become a learner after a further thirty seconds of CS alone; those network realizations that did not reach the threshold value are defined here as non-learners. Even with this more realistic assumption, the large majority of network realizations were learners (19 out of 20 networks). This is in agreement with the experimental fear conditioning literature showing that most of the subjects learn the association between CS and fear after only one trial (Schafe et al., 2000).

Physiology of interneuron types is critical to their role in appropriate plasticity

The PV cell is necessary to induce the correct pre-post timing between ECS and F needed for long-term potentiation of the ECS to F conductance. In our model, PV has reciprocal connections with F and provides lateral inhibition to ECS. Since the lateral inhibition is weaker than the feedback inhibition, PV tends to bias ECS to fire before F. This creates the fine timing needed for the depression-dominated rule to instantiate plasticity. If we used the classical Hebbian plasticity rule (Bi and Poo, 2001) with gamma frequency inputs, this fine timing would not be needed and ECS to F would potentiate over most of the gamma cycle, and thus we would expect random timing to lead to potentiation (Fig. S2). In this case, no interneurons are needed (See Discussion).

In this network configuration, the pre-post timing for ECS and F is repeated robustly over time due to coordinated gamma oscillations (PING, as shown in Fig. 4A, Fig. 1C, and Fig. 2B,C) arising through the reciprocal interactions between F and PV (Feng et al., 2019). PING can arise only when PV is in a sufficiently low excitation regime such that F can control PV activity (Börgers et al., 2005), as in Fig. 4A. However, despite PING establishing the correct spike-timing for potentiation, setting PV in a low excitation regime in which F-PV PING occurs continuously does not lead to potentiation (Fig. 4A, Fig. S1C): the depression-dominated rule leads to depression rather than potentiation unless the PING is periodically interrupted. During the pauses, made possible only in the full network by the presence of VIP and SOM, the history-dependent build-up of depression decays back to baseline, allowing potentiation to occur on the next ECS/F active phase. (The detailed mechanism of how this happens is in the Supplementary Information, including Fig. S1). Thus, a network without the other interneuron types cannot lead to potentiation.

Figure 4.

In our simulations (except for the one in Fig. 4A), we choose to set PV cells at a high excitation level due to a strong input from CS in order to account for: (i) the experimental literature showing a strong activation of PV after the onset of CS (Wolff et al., 2014); (ii) many other inputs that the PV cell receive besides CS. We show in Fig. S3A that increasing the PV excitation by adding other inputs to PV cells leads to similar results; all interneurons are necessary for ECS to F potentiation under these circumstances. If CS activation of PV is weakened, we find potentiation can occur in the absence of VIP and SOM (Fig. S3B). However, we consider this outside the physiological range.

In the context of high PV excitation, we find that the VIP cell plays two critical roles in enabling potentiation to occur via PING. First, it reduces the excitation of the PV cell during its active low-theta phase, enabling the PV cell to participate in PING (Fig. 4B, Fig. S1D). Second, it provides the periodic interruptions in ECS and F firing necessary for potentiation during its silent low-theta phase. Without these pauses, depression dominates (see SI section “ECS and F activity patterns determine overall potentiation or depression”). This interruption requires the participation of another interneuron, the SOM cell (Figs. 2B, 4B, S1): the pauses in inhibition from the VIP periodically interrupt ECS and F firing by releasing PV and SOM from inhibition and thus indirectly silencing ECS and F.

Network with multiple heterogeneous neurons can establish the association between CS and fear

To test the robustness of our single-cell results to heterogeneity, we expand our BLA network to include three cells of each interneuron subtype and ten of each excitatory projection neuron (Fig. 5A). Each neuron has independent noise and cellular parameters (see Materials and Methods for details). We find that the network very robustly produces potentiation between the ECS and F receiving CS and US, respectively, during fear training, with additional properties matching the experimental results.

Figure 5.

Fig. 5B shows an example of the network dynamics during fear conditioning with simultaneous CS and US inputs. As previously presented for the single neuron network (Figs. 2-3), interneurons are crucial in conferring the correct pre-post spike timing to ECS and F. We assume all the VIP interneurons receive the same US; hence, the VIP neurons tend to approximately synchronize at low theta in response to US input, allowing a window for potentiation of ECS to F conductance (Fig. 5B), as we have seen in the single cell model. The potentiation from ECS to F is specifically for the excitatory projection cells receiving CS and US, respectively (Fig. 5C). The ECS cells not receiving CS are inhibited by ongoing PV activity during the disinhibition window (Fig. 5B); they are constructed to be firing at 11 Hz in the absence of any connections from other cells. The lack of activity in those cells during fear conditioning implies that there is no plasticity from those ECS cells to F. This larger network corroborates the results obtained for the single neuron network: only the realizations of the full network learn the association between CS and fear (Fig. 5D, left) and all those network realizations become learners in less than forty seconds (Fig. 5D, right). Similarly, there is a striking failure of plasticity if any interneuron is removed from the network; even partial plasticity does not arise.

Increased low-theta frequency is a biomarker of fear learning

In this section, we explore how plasticity in the fear circuit affects the network dynamics, comparing after fear conditioning to before. We first show that fear conditioning leads to an increase in LFP low theta frequency power compared to pre-conditioned level (Fig. 6A,B); there is no change in the high theta power. These results reproduce the experimental findings in (Davis et al., 2017). We find that the newly potentiated AMPA current from ECS to F is the main source of increase in the low theta power (Fig. 6C). The increase in power is in the low theta range because ECS and F are allowed to spike only during the active phase of the low theta spiking VIP neurons.

Figure 6.

Although the experimental results in (Davis et al., 2017) show an increase in low theta after fear learning, they excluded non-learners from their analysis, and thus is unclear from their results if low theta can be used as a biomarker of fear learning. To address this question, we looked at the power spectra after fear conditioning in learning versus non-learning networks. To have an adequate number of non-learners, we ran 60 network simulations for 10 seconds and chose those whose conductance from ECS to F remained lower than 0.037 mS/cm³ (i.e., non-learners). Notably, the low theta power increase is completely absent after fear conditioning in those network realizations that display no signs of learning. This suggests low theta power change is not just an epiphenomenon but rather a biomarker of successful fear conditioning.

Finally, we find a trend towards increased power in the gamma range between 35 and 60 Hz compared to pre-conditioned level, which supports experimental results in (Courtin et al., 2014b).

Plasticity between fear neuron and VIP slows down overall potentiation

In-vivo studies documented plastic connections involving GABAergic interneurons in addition to those among excitatory projection neurons (Szinyei et al., 2007; Lucas et al., 2016). VIP interneurons play a crucial role in enabling the association between CS and fear in the fear circuit; hence, we studied the effect that potentiation in the conductance from the fear neuron F to the VIP interneurons, acquired during fear conditioning, has on the potentiation from ECS to F. We find that such a plastic connection in both the single-neuron network and the heterogeneous network can lead to preventing some realizations from learning the association between CS and fear in the first forty seconds; this is because strengthening of the F to VIP connections slow down the potentiation from ECS to F (Fig. S4-S5). Further details on the F to VIP plasticity rule and the role of the F to VIP conductance in slowing down the association between CS and fear can be found in section “Plasticity between fear neuron and VIP slows down overall potentiation” in the Supplementary Information.

We also find that plasticity between the fear neuron and VIP interneurons helps increase even further the low theta power after fear conditioning compared to before (Fig. S6). Further details can be found in section “Plastic F to VIP connections further increase low-theta frequency power after fear conditioning” in the Supplementary Information.

Discussion

Overview

Our study suggests that amygdalar rhythms play a crucial role in plasticity during fear conditioning. Prominent rhythms found in the BLA during fear conditioning include low theta (3-6 Hz), high theta (6-12 Hz), and gamma (>30 Hz) (Seidenbecher et al., 2003; Courtin et al., 2014b; Stujenske et al., 2014; Davis et al., 2017). Experimental work in rodents shows that the BLA undergoes more oscillatory firing at low theta frequency and gamma following fear conditioning, while high theta frequency remains unchanged compared to before fear conditioning (Courtin et al., 2014b; Davis et al., 2017).

To examine the origin of these rhythmic changes and their functional role in fear conditioning, we implement a biophysically detailed model of the BLA. We show that vasoactive intestinal peptide-expressing (VIP) interneurons, parvalbumin-expressing (PV) interneurons, and somatostatin-expressing (SOM) interneurons in the BLA may be centrally involved in producing the experimentally measured rhythms based on the biophysical properties of these interneurons. More specifically: the gamma oscillation is associated with PV cell interaction with excitatory cells to produce PING; VIP cells produce low theta due to their intrinsic D-current; SOM cells help produce the high theta rhythm due to persistent sodium current (NaP-current) and H-currents. Moreover, we show that the rhythmic dynamics produced by VIP and PV cells both play a crucial role in the instantiation of plasticity during associative learning by promoting the formation of a dedicated fear circuit shaped through spike-timing-dependence and the removal of any of the interneuron types from the circuit leads to failure of the plasticity needed for associative fear learning. We note that the presence of SOM cells is crucial for plasticity in our model since they help to produce the necessary pauses in the excitatory projection cell activity. However, the high theta rhythm they produce is not crucial to the plasticity: in our model, high theta or higher frequency rhythms in SOM cells are all conducive to associative fear learning. This opens the possibility that the high theta rhythm in the BLA mostly originates in the prefrontal cortex and/or the hippocampus (Stujenske et al., 2014, 2022). Finally, we replicate the experimental increase in the low theta after successful training (Davis et al., 2017) and determine that this is a biomarker of fear learning, i.e., the increase in low-theta power does not appear in the absence of learning.

Synaptic plasticity in our model

Synaptic plasticity is the mechanism underlying the association between neurons that respond to the neutral stimulus CS (ECS) and those that respond to fear (F), which instantiates the acquisition and expression of fear behavior. One form of experimentally observed long-term synaptic plasticity is spike-timing-dependent plasticity (STDP), which defines the amount of potentiation and depression for each pair of pre- and postsynaptic neuron spikes as a function of their relative timing (Bi and Poo, 2001; Caporale and Dan, 2008). There are many STDP rules in the literature (Abbott and Nelson, 2000; Feldman, 2012) but none we are aware of specifically for the amygdala. The depression-dominated Hebbian rule, which we use in this study, provides each interneuron subtype with a specific role in allowing the instantiation of fear learning. With this rule, depression is overall stronger than potentiation, unless there is fine timing among pyramidal cells in the gamma and low theta cycles.

With the depression-dominated rule, each interneuron plays a role in setting up this fine timing in the gamma and low theta cycles needed for appropriate plasticity. PV interneurons are important in creating a gamma rhythm with F when the latter is activated by the unconditioned stimulus (US). This gamma rhythm plays a central role in the fine timing between ECS and F. As mentioned earlier, the critical requirements for plasticity are a pause of VIP firing within each low-theta cycle and fine timing of the ECS (pre) and F (post) neurons in the active phase of the cycle. The pause in VIP cell firing allows for activation of SOM cells, which helps to inhibit the ECS and F cells. Without the pause, the ECS and F cells continue to fire. It is a consequence of the depression-dominated rule that this would lead to depression; this happens because the depression does not relax to zero between gamma-frequency spikes, and hence continues to build up. With a Hebbian plasticity rule characterized by a lower amplitude for depression than for potentiation, potentiation would occur at most of the phases of gamma and thus fine timing would not be needed (Fig. S2). The depression-dominated rule allows for regulation of the plasticity by modulation of all the kinds of cells that are known to be involved in fear learning as well as providing reasons for the known involvement of rhythms. It remains to be understood how lack of such regulation might produce deficits in fear learning.

VIP cells may be important in establishing the kind of depression-dominated rule we are using in this model for the synapse between ECS and F. VIP cells are known to corelease the peptide VIP along with GABA (Bayraktar et al., 1997). The amount of the release is related to the amount of high frequency firing, i.e., the duty cycle of the gamma burst in each low theta cycle (Agoston et al., 1988; Agoston and Lisziewicz, 1989). The peptide VIP can act on second messenger pathways to inhibit potentiation; the pathways are complex and not fully understood, but in hippocampus they involve GABA transmission, NMDA activation, and CaMKII (Caulino-Rocha et al., 2022). The relevant VIP receptor VPAC1 is known to exist in the amygdala (Joo et al., 2005; Boucher et al., 2021). Thus, by inhibiting potentiation, VIP may change a Hebbian plasticity rule to a depression-dominated rule.

Involvement of other brain structures

Studies using fear conditioning as a model of associative learning reveal that learning and expression of fear are not limited to the amygdala but involve a distributed network including the amygdala, the medial prefrontal cortex, and the hippocampus (Seidenbecher et al., 2003; Bocchio and Capogna, 2014; Courtin et al., 2014a; Stujenske et al., 2014; Tovote et al., 2015; Karalis et al., 2016; Chen et al., 2021). In our model, the fear-encoding neuron F fires for at least 13 seconds due to the US (first ∼2 seconds), followed by the fear memory that could excite F in a similar way to US. This is suggested by the observation that during fear conditioning, animals are hyperactive for about 10 seconds after US presentation and then freeze while they are maintained in the box in which they received the shock (Hole and Lorens, 1975). The medial prefrontal cortex and the hippocampus may provide the substrates for the continued firing of F needed for ECS to F potentiation after the 2-second US stimulation.

It has also been suggested that ventral tegmental area has a role in fear expression (Lesas et al., 2023). Furthermore, it has been reported that the prelimbic cortex (PL) modulates the BLA SOM cells during fear retrieval, and the latter cells are crucial to discriminate non-threatening cues when desynchronized by the PL inputs (Stujenske et al., 2022). Other brain structures, e.g., prefrontal cortex and hippocampus, have been documented to play a crucial role also in fear extinction, the paradigm following fear conditioning aimed at decrementing the conditioned fearful response through repeated presentations of the CS alone. As reported by several studies, fear extinction suppresses the fear memory through the acquisition of a distinct memory, instead of through the erasure of the fear memory itself (Harris et al., 2000; Bouton, 2002; Trouche et al., 2013; Thompson et al., 2018). Davis et al., 2017 found a high theta rhythm following fear extinction that was associated with the suppression of threat in rodents. In the future, our BLA network will also include structures in the prefrontal cortex and hippocampus, thus allowing us to further study the role of rhythms in fear conditioning and extinction. We hypothesize that a new population of PV interneurons plays a crucial role in mediating competition between fearful memories, associated with a low theta rhythm, and safety memories, associated with a high theta rhythm; supporting experimental evidence is in (Lucas et al., 2016; Davis et al., 2017; Chen et al., 2022).

Comparison with other models

Our modeling approach aims at investigating the genesis and the functional role of rhythms in the amygdala fear circuit before, during, and after fear learning. Several computational models related to fear conditioning have been proposed in the last years, including connectionist models, firing rate models, and biophysically detailed ones (for a review see (Nair et al., 2016)). However, to the best of our knowledge, none focuses on rhythms generated by different kinds of interneurons. Among the biophysically detailed models, some tend to focus on the connectivity among a large number of similar neurons rather than the interactions among different kinds of neurons (Feng et al., 2016; Kim et al., 2016). Our focus, using a smaller number of neurons, highlights the functional importance of the physiological differences among the different kinds of neurons in producing fine timing needed for plasticity. Feng et al., 2019 also suggests that the gamma in BLA is a PING; our paper looks at how this rhythm evolves during fear conditioning while also investigating interactions with other BLA rhythms.

The dynamics of the VIP cell play a central role in the plasticity we investigate. This is in contrast to the cortical model in (Veit et al., 2023) for which VIP is essential for locally controlling gain and globally controlling coherence in gamma oscillations. In the model by Veit et al., 2023, the global control requires differences in long-range connectivity that are known to exist and are inserted in the model by hypothesis. Our paper shows how more detailed biophysics produces rhythms among the interneurons used in (Veit et al., 2023) and how these rhythms can produce the plasticity needed to construct those differences in long-range connectivity. Thus, though (Veit et al., 2023) shows that rhythms are not needed for some kinds of control once connectivity is established, our paper suggests that the same set of interneurons, with more detailed physiology, can support the establishment of appropriate connectivity as well as the control described in (Veit et al., 2023). We note that Veit et al., 2023 deal with cortical networks, while our model describes BLA networks; however, it is known that these networks are structurally related (Sah et al., 2003; Tovote et al., 2015; Polepalli et al., 2020).

The study in (Grewe et al., 2017) suggests that fear conditioning induces both potentiation and depression among coactive excitatory neurons; coactivity was determined by calcium signaling and thus did not allow measurements of fine timing between spikes. In our model, depression happens in several circumstances: (1) in the absence of pauses of ECS and F firing and (2) in the absence of the fine timing provided by the PV neurons. This could happen among some cells responding to weaker sensory inputs that do not lead to pre-post timing with fear neurons. This timing could be modified by the “tri-conditional rule”, as suggested in (Grewe et al., 2017).

Limitations and caveats

The use of small number of neurons raises the issue of model scalability. One way in which large models can be different from much smaller ones is in heterogeneity of the neurons of any given type. By using a network with a few neurons, we have begun the study of effects of heterogeneity. In general, the use of small numbers, which effectively assumes that each cell represents a synchronous subset of a larger population, replaces the use of gap-junctions that are known to exist in the cortex among VIP cells (Francavilla et al., 2018) as well as among SOM cells and among PV cells (Tremblay et al., 2016).

We do not explicitly model the biophysics of NMDA receptors. Rather, we model the effect of such receptors using the spike-timing-dependent plasticity resulting from such biophysics, as is commonly done when modeling STDP (Song et al., 2000). Also, our neurons are single-compartment, so do not build in the spatial structure known to exist on the dendrites (Blair et al., 2001; Bennett et al., 2019).

Our model assumes that initial stages of fear learning can be accomplished entirely within the amygdala, though it is known that other structures in the brain are important for modulating networks related to fear. Much of the work involving the prefrontal cortex and the hippocampus relates to fear extinction, which is not addressed in this paper.

Summary and significance

We have shown how networks of amygdala neurons, including multiple types of interneurons, can work together to produce plasticity needed for fear learning. The coordination necessary to produce the plasticity requires the involvement of multiple rhythms. Thus, our paper both accounts for the experimental evidence showing such amygdala rhythms exist and points to their central role in the mechanisms of plasticity involved in associative learning. These mechanisms may be common to other types of associative learning, as similar interneuron subtypes and connectivity are ubiquitous in the cortex (Sah et al., 2003; Tovote et al., 2015; Polepalli et al., 2020).

Acknowledgements

Research reported here was supported by the National Institute of Neurological Disorders and Stroke of the National Institute of Health under Award Number U01NS122082 to NK and DBA.

Author contributions

Conceptualization, AC, DBA, MM, and NK; Methodology, AC, MM, and NK; Software, AC, and MM; Formal analysis, AC, and MM; Investigation, AC, MM, and MM; Data curation, AC; Writing – original draft, AC, MM, and NK; Visualization, AC, and DBA; Supervision, MM and NK; Funding acquisition: DBA and NK.

Declaration of interests

The authors declare no competing interest.

Materials and Methods

Neuron model

Our network is made of interacting single-compartment neurons modeled using conductance-based models with Hodgkin-Huxley-type dynamics. The temporal voltage change of each neuron is described by:

where, c_m is the membrane capacitance, and I_membrane are the intrinsic membrane currents, which include a fast sodium current (I_Na), a fast potassium current (I_K), and a leak current (I_L) for all neuron types. VIP interneurons additionally have a D-current, and SOM interneurons additionally have P and H-currents (Rotstein et al., 2005; Tort et al., 2007). All these currents are discussed in more detail below, where we describe each neuron individually. The synaptic currents (I_synaptic) take into account the input from the other neurons in the network and depend on the network connectivity and specific type of synaptic input, as discussed below. Finally, the background drive I_app is a constant term that determines the background excitation of a neuron, and I_noise corresponds to a Gaussian noise input with mean zero, standard deviation 1, and a specific amplitude for each neuronal cell type (specified below).

Membrane currents

The membrane currents I_Na, I_K, and I_L are modeled using Hodgkin-Huxley-type conductances formulated as:

Each membrane current has a constant maximal conductance and a reversal potential E_channel (for channel = Na, K, or L). The activation (m and n) and inactivation (h) gating variables evolve in time according to:

where x = m, n, h. The steady-state function (x_∞) and the time constant of decay (τ_x), which are taken from previous models (Mainen and Sejnowski, 1996; Olufsen et al., 2003), are formulated as rate functions for each opening (α_x) and closing (β_x) of the ionic channel through:

The specific functions and constants for each cell type in the network are given below.

Vasoactive intestinal peptide interneurons (VIP)

The membrane currents (I_membrane) of the VIP interneurons consist of a fast sodium current (I_Na) (described as in the first formulation of I_Na in Eq. (1)), a fast potassium current (I_K), a leak current (I_L), as in Eq. (1), and a potassium D-current (I_D). The formulations of these currents were derived from a previous model of cortical interneurons (Golomb et al., 2007) and subsequently used to model striatal fast spiking interneurons (Sciamanna and Wilson, 2011; Chartove et al., 2020), which are reported below.

The maximal sodium conductance is and the sodium reversal potential is E_Na = 50 mV. The steady state functions for the sodium current activation (m) and inactivation (h) variables and h time constant (τ_h) are described by:

The maximal conductance for the fast potassium channel is and the potassium reversal potential is E_K = −90 mV. The fast potassium channel has no inactivation gates and two activation gates described as follows:

The leak current (I_L) has no gating variables. The maximal leak conductance is and the leak channel reversal potential is E_L = −70 mV.

The fast-activating, slowly inactivating potassium D-current I_D is formulated as described in (Golomb et al., 2007):

with maximal conductance . The steady state functions for the activation (a) and inactivation (b) variables are described as follows:

while the time constant of the decay is τ_α = 2ms for the activation gate and τ_β = 150ms for the inactivation gate.

In the absence of US, the applied current I_app is set to 4 μA/cm². When US is present, I_app = 5 μA/cm². The Gaussian noise (I_noise) has mean 0, standard deviation 1, and an amplitude of , where δt = 0.05 ms corresponds to the time step of integration in our simulations.

Somatostatin-positive interneurons (SOM)

The membrane currents (I_membrane) of the SOM interneurons consist of a fast sodium current (I_Na) (described as in the second formulation of I_Na in Eq. (1)), a fast potassium current (I_K), and a leak current (I_L) as in Eq.(1), along with an H-current (I_H) and NaP-current (I_P). The formulations of these currents were taken from previous models of the oriens lacunosum-moleculare (SOM-positive O-LM) cells in the hippocampus (Rotstein et al., 2005; Tort et al., 2007) and are reported below.

The maximal sodium conductance is and the sodium reversal potential is E_Na = 55 mV. The rate functions for the sodium current activation (m) and inactivation (h) variables in Eqs. (2-3) are formulated as follows:

The maximal potassium conductance is and the potassium reversal potential is E_K = −90 mV. The rate functions for the potassium current activation variable (n) are formulated as follows:

The leak current (I_L) has no gating variables. The maximal leak conductance is and the leak channel reversal potential is E_L = −65 mV.

The slow hyperpolarization-activated mixed cation current I_H is formulated as described in (Rotstein et al., 2005):

with maximal conductance and E_H = −20 mV. The steady state functions for the h^f and h^s variables and their time constant of decay are described as follows:

The persistent sodium current I_P is formulated as described in (Rotstein et al., 2005, 2006):

The maximal persistent sodium conductance is and the sodium reversal potential is, as stated above, E_Na = 55 mV. The steady state function for the persistent sodium current I_P (p_∞) and the time constant (τ_p) are described by:

Throughout all simulations, the applied current I_app is set to 0.1 μA/cm². The Gaussian noise (I_noise) has mean 0, standard deviation 1, and an amplitude of , where δt = 0.05 ms corresponds to the time step of integration in our simulations. We note that the persistent sodium current can be replaced by an A-current to produce a high theta rhythm (Gloveli et al., 2005).

Parvalbumin-positive interneurons (PV)

The membrane currents (I_membrane) of the PV interneurons consist of only a fast sodium current (I_Na) (described as in the second formulation of I_Na in Eq. (1)), a fast potassium current (I_K), and a leak current (I_L), as in Eq. (1).

The maximal sodium conductance is and the sodium reversal potential is E_Na = 50 mV. The rate functions for the sodium current activation (m) and inactivation (h) variables are formulated as follows:

The maximal potassium conductance is and the potassium reversal potential is E_K = −100 mV. The rate functions for the potassium current activation (n) variables are formulated as follows:

The leak current (I_L) has no gating variables. The maximal leak conductance is and the leak channel reversal potential is E_L = −67 mV. Throughout all the simulations, I_app = 0 μA/cm². The Gaussian noise (I_noise) has mean 0, standard deviation 1, and an amplitude of , where δt = 0.05 ms corresponds to the time step of integration in our simulations.

Excitatory projection neurons (ECS and F)

The membrane currents (I_membrane) of ECS and F consist of a fast sodium current (I_Na) (described as in the first formulation of I_Na in Eq. (1)), a fast potassium current (I_K), and a leak current (I_L) as in Eq. (1).

The maximal sodium conductance is and the sodium reversal potential is E_Na = 50 mV. The rate functions for the sodium current activation (m) and inactivation (h) variables are formulated as follows:

The maximal potassium conductance is and the potassium reversal potential is E_K = −100 mV. The rate functions for the potassium current activation (n) variables are formulated as follows:

The leak current (I₀) has no gating variables. The maximal leak conductance is and the leak channel reversal potential is E_L = −67 mV. The formulations of these currents were taken from the description of excitatory/inhibitory neurons presented in (Zhou et al., 2018).

When neither US nor CS are injected, the applied current I_app,F is set to 0.35 μA/cm² and I_app,ECS is set to 0.45 μA/cm². By contrast, I_app,F is set to 0.5 μA/cm² when US is injected. For both ECS and F, the Gaussian noise (I_noise) has mean 0, standard deviation 1, and an amplitude of , where δt = 0.05 ms corresponds to the time step of integration in our simulations.

Conditioned and unconditioned stimuli

The conditioned (CS) and unconditioned (US) stimuli affect specific cell types according to the fear conditioning phase. CS consists of a Poisson spike train (λ = 800) that excites an auxiliary excitatory neuron (described by the same equations used for ECS in the previous section). The auxiliary excitatory neuron excites both ECS and PV and makes them fire, in isolation, at ∼50 Hz. The maximal AMPA conductance from the auxiliary excitatory neuron to PV and ECS is (see the next paragraph for a description of the AMPA synapses). Many AMPA conductance values from the auxiliary excitatory neuron to PV (spanning from 0.01 mS/cm² to 0.2 mS/cm²) have been considered to make Fig. S3B Similarly, US (λ = 800, independent of CS) affects an auxiliary excitatory neuron that makes F fire in isolation fires at ∼50 Hz. The maximal AMPA conductance from the auxiliary excitatory neuron to F is (see the next paragraph for a description of the AMPA synapses). Finally, US influences VIP activity by increasing its I_app set to 5 μA/cm².

Network connectivity and synaptic currents

We modeled the network connectivity as presented in Fig. 2B, derived from (Krabbe et al., 2019). We have a total of 9 types of projections between neurons: 6 inhibitory (VIP → PV, VIP → SOM, PV → F, PV → ECS, SOM → F, SOM → ECS), 3 excitatory (ECS → F, F → PV, F → VIP).

All inhibitory synapses are described as GABAa currents (I_GABAa) using a Hodgkin-Huxley-type conductance, as formulated in (Olufsen et al., 2003):

The maximal GABAa conductance VIP → PV is , VIP → SOM is , PV → F is , PV → ECS is , SOM → F is , and SOM → ECS is , where N_VIP, N_PV, N_SOM are the number of VIP, PV, and SOM cells, respectively, in the network. The GABAa current reversal potential (E_i) is set to −80 mV. The variable s_i represents the gating variable for inhibitory GABAa synaptic transmission, where i stands for inhibitory synapse. The contribution of an inhibitory synapse to a specific postsynaptic inhibitory neuron j in the network takes the following form:

The contribution to a specific postsynaptic excitatory neuron m in the network, reads as follows:

where k indexes the presynaptic inhibitory neurons. The variable in Eq. (4) describes the kinetics of the gating variables from the inhibitory presynaptic neuron k^th to the inhibitory postsynaptic neuron j. This variable evolves in time according to:

Similarly, the kinetics of the gating activation variable in Eq. (5) from the k^th interneuron to the postsynaptic excitatory neuron m is formulated as:

The GABAa decay time constant is a constant that depends on the type of presynaptic interneuron. The rate functions for the open state of the GABAa receptor (g_GABAa(V_k)) has a specific form based on the presynaptic cell type k. More specifically,

All excitatory synapses are described as AMPA currents (I_AMPA) using a Hodgkin-Huxley-type conductance, as formulated in (Olufsen et al., 2003):

At the beginning of the fear conditioning paradigm there is no connection from ECS to F, i.e., the maximal AMPA conductance ECS → F is . Since ECS to F is a plastic connection (see paragraph related to synaptic plasticity), it evolves over time up to a maximum of . The maximal AMPA conductance F → PV is , and F → VIP is . In the case of the plastic F to VIP cell connections (see Supplementary Information), the F → VIP conductances evolve over time up to a maximum of .The AMPA current reversal potential (E_e) is set to 0 mV.

The variable s_e represents the gating variable for excitatory AMPA synaptic transmission, where e stands for excitatory synapse. For a specific postsynaptic excitatory neuron m in the network:

For a specific postsynaptic inhibitory neuron j in the network:

where k indexes the presynaptic excitatory neurons.

The variable describes the kinetics of the gating variables from the excitatory presynaptic neuron k^th to the excitatory postsynaptic neuron m. This variable evolves in time according to:

Similarly, the kinetics of the synaptic activation variable of the k^th excitatory neuron to the inhibitory neuron j is denoted by and is formulated as:

The time-constant of decay for the AMPA synapse is τ_e = 2 ms. The rate functions for the open state of the AMPA receptor (g_AMPA(V_k)) follows the mathematical formulation:

Synaptic plasticity

Fear conditioning is a paradigm able to create associative learning between the neutral (CS) and the aversive (US) stimuli. Synaptic plasticity is thought to be at the basis of associative learning. In our work, synaptic plasticity takes the form of spike-timing-dependent plasticity (Song et al., 2000; Lee et al., 2009), where synaptic modifications are enforced at the synapse from ECS to F at each presynaptic (ECS) and postsynaptic (F) neuron spikes. Synaptic modification is generated in the model through two auxiliary functions: P, used for potentiation, and M, used for depression, as in standard STDP models. The major difference between our use of STDP and some others is that, for each pre-synaptic spike, we take into account all the post-synaptic spikes with which it could potentially interact; similarly for each post-synaptic spike, we take into account all presynaptic spikes. P and M are initialized at zero and are updated at each presynaptic and postsynaptic neuron spike, respectively. Between spikes, they exponentially decay to zero (see Fig. S1). The update is described as follows. At each postsynaptic neuron (F) spike, M(t) is decremented by an amount A₋ = 0.005, i.e., M(t) = M(t) − A_N. At each presynaptic neuron (ECS) spike, P(t) is incremented by an amount A₊ = 0.005, i.e., P(t) = P(t) + A₊ (Fig.S1A). This Hebbian plasticity rule is depression-dominant when A₊ = A₋ and τ₊ < τ₋, as formulated in this work. The exponential decay is described by the following equations:

with τ₋ = 28 ms and τ₊ = 14 ms. Every time the synapse receives an ECS action potential at time t, its maximal conductance is weakened according to . If is set to zero. If F fires an action potential at time t, the synapse maximal conductance is strengthened according to . If this strengthening makes , then is set to (see Fig. S1A, right, for an example). See Supplementary Information for a visualization of M and P and consequences of the depression-dominated rule.

In the Supplementary Information we introduce a plastic connection from the fear neuron F to the VIP interneurons. For that specific synapse, A₊ = 0.00065 and A₋ = 0.0003, while the decay time constants are as above. In this case, the rule is not depression-dominant.

Model simulations

Our network models were programmed in C++ and compiled using g++ compiler (Apple clang version 14.0.0) on macOS Monterey version 12.5.1. The differential equations were integrated using a fourth-order Runge Kutta algorithm. The integration time step was 0.05 ms. Model output is graphed and analyzed using MATLAB, Version R2022a. Simulation codes are made freely available (see Resource Availability section).

Local field potentials and spectral analysis

Modeling LFP

One measure of neuronal population activity in the BLA is the LFP. It has been suggested that inhibitory synaptic currents represent a major contribution to LFP in the hippocampus and the cortex because of the neurons arranged in a laminar way (Bazelot et al., 2010; Teleńczuk et al., 2017, 2020). By contrast, LFP in brain regions with no laminar structure (as the BLA), may be generated by localized excitatory inputs and intrinsic membrane currents (Reimann et al., 2013; Tanaka and Nakamura, 2019). For this reason, we model the LFP as the sum of all synaptic AMPA currents in the pyramidal cells, along with the D-current in the VIP interneurons, and NaP-current and H-current in SOM interneurons.

Spectral analysis

Stationarity of the network before and after fear conditioning is ensured after 2000 ms. Thus, to ensure elimination of transients due to initial conditions, we discard the first 2000 ms of LFP signals. LFP’s power spectra are calculated using the Thomson’s multitaper power spectral density estimate (MATLAB function pmtm) (Bokil et al., 2007) for frequencies ranging from 0.1 – 70 Hz. Analysis codes are made freely available (see Resource Availability section).

Resource availability

Data and code availability

The computer simulation code created for this paper, as well as the code to make the figures presented, are available online at https://github.com/annacatt/Basolateral_amygdala_oscillations_enable_fear_learning_in_a_biophysical_model). The simulated data that support the findings of this study are available at https://datadryad.org/stash/share/FCSYmUOhwforkleWFwlKW2hhmu2XXA7hsvsYbWnQGCU (the following DOI will be active after publication: https://doi.org/10.5061/dryad.pvmcvdnr2)