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 open questions 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. The presentation of the US and CS may or may not overlap in time (e.g., see (Laxmi et al., 2003; Stujenske et al., 2014; Krabbe et al., 2019)).

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 experimental studies (Antonoudiou et al., 2022; Bratsch-Prince et al., 2024) suggest that BLA can intrinsically generate theta oscillations (3-12 Hz). Furthermore, other studies show increased low theta (Davis et al., 2017) and gamma (Courtin et al., 2014b) in 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. However, the potential roles of rhythms in instantiating plasticity needed for successful learning are still under investigation (Bocchio et al., 2017). We also note that there is not uniformity on the exact frequencies associated with low and high theta, e.g., ((Lorétan et al., 2004) used 2-6 Hz for low theta). Here, we use 2-6 Hz for the theta range and 6-14 Hz for the high theta range.

In this paper, we aim to show 1) How a variety of BLA interneurons (PV, SOM and VIP) lead to the creation of these rhythms and 2) How the interaction of the interneurons and the rhythms leads to the appropriate timing of the cells responding to the US and those responding to the CS to promote fear association through spike-timing-dependent plasticity (STDP). Since STDP requires overlap of the effects of the CS and US, and some conditioning paradigms do not have overlapping US and CS, we include as a hypothesis that the effects of the CS and US overlap even if the CS and US stimuli do not. In the Discussion, we suggest how neuromodulation by ACh and/or dopamine can provide such overlap. We create a biophysically detailed model of the BLA circuit involving all three types of interneurons and show how each may participate in producing the experimentally observed rhythms and interacting to produce the necessary timing for fear learning. In particular, we find that the low theta, high 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 STDP. 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 the entire history of spikes, not just the most recent pre- and post-synaptic spikes. 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 (Davis et al., 2017). 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 is a consequence of network dynamics that emerge due to the newly formed synapse linking CS-encoding and fear-encoding neurons. Finally, we note that the ideas in the model may apply very generally to associative learning in the cortex, which contains similar subcircuits of pyramidal cells and interneurons: PV, SOM and VIP cells.

Results

Rhythms in the BLA can be produced by interneurons

Brain rhythms are thought to be generated 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 VIP interneurons in (Sosulina et al., 2010) (Fig. 1A, top): in the baseline condition, the condition without any external input from the fear conditioning paradigm (Fig. 1B, top), our VIP neurons exhibit short bursts of gamma activity (∼38 Hz) at low theta frequencies (∼2-6 Hz) (peaking at ∼3.5 Hz) (see Fig. S1A).

Figure 1.

For SOM interneurons, we focus on 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. SOM interneurons have been mostly 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), which are known to have a hyperpolarization-activated current, i.e., H-current, and a persistent sodium current, i.e., NaP-current. An H-current has also been observed in the SOM cells in the BLA (Ünal et al., 2020). In our model, with the introduction of both H- and NaP-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) (Fig. 1A, middle). In our baseline model, SOM cells have a natural frequency of ∼12 Hz (Fig. 1B, middle; Fig. S1B), which is at the upper limit of the experimental high theta range; this motivates our choice to extend the high theta range up to 14 Hz in order to include the peak.

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 (Fig. 1B, bottom) and show similar behaviors to type IV interneurons (Sosulina et al., 2010) in response to depolarizing and hyperpolarizing currents (Fig. 1A, bottom). However, when reciprocal connections are present between the PV interneuron and the excitatory projection neuron, these neurons form a PING rhythm (pyramidal-interneuron network gamma) (Whittington et al., 2000), if the excitatory projection neuron receives enough excitation to fire (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 frequency of PING 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). 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). The US input on the pyramidal cell and VIP interneuron is modeled as a Poisson spike train at ∼ 50 Hz and an applied current, respectively. In the rest of the paper, we will use the words “US” as shorthand for “the effects of US”. As shown in Fig. 2A (top, right), VIP disinhibits F by inhibiting both SOM and PV, as suggested in (Krabbe et al., 2019). 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). The simplest network we consider is made of one neuron for each cell type. We introduce additional neurons for each cell type with some heterogeneity in the last two sections of the Results.

Figure 2.

Fig. 2A shows a typical dynamic of the network before and after the US input onset; The network produces all the rhythms originating from the interneurons alone or through their interactions with the excitatory projection neuron (shown in Fig. 1). Specifically, since VIP is active at low theta during both rest and upon the injection of US, it then modulates F at low theta cycles via SOM and PV. In the baseline condition, the VIP interneuron has short gamma bursts nested in a low theta rhythm. With US onset, VIP increases its burst duration and the frequency of the 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 the activity of F right after the US onset. Finally, in Fig. 2A, PV lacks any external input and fires only when excited by F. Due 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. Similarly to the US, in the rest of the paper, we will use the words “CS” as shorthand for “the effects of CS”. In our simulations, CS is modeled as a Poisson spike train at ∼ 50 Hz, independent of the US input. Thus, we hypothesize that the time structure of the inputs sometimes used for the training (e.g., a series of auditory pips) is not central to the formation of the plasticity in the network. Our CS input describes either the context or the cue in contextual and cued fear conditioning, respectively. For the context, the input may come from the hippocampus or other non-sensory regions, but this does not affect its role as input in the model. 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. S2A for more details).

Fig. 2B 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 simultaneous effects of the 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. 2 B,C and Fig. S2B); 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 (during CS-only), SOM and PV prevent plasticity by silencing F, which in the absence of US receives only a weak applied current.

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 has both ECS and F, as well as VIP and PV interneurons, active during the entire period: an initial bout of US is known to produce a long-lasting fear response beyond the offset of the US (Hole and Lorens, 1975) and to induce the release of neuromodulators. The latter, in particular acetylcholine and dopamine that are known to be released upon US presentation (Harmer and Phillips, 1999; Suzuki et al., 2002; Rajebhosale et al., 2024a), may induce more sustained activity in the ECS, F, VIP, and PV neurons during and after the presentation of US, thus ensuring a concomitant activation of those neurons necessary for STDP to take place (see “Assumptions and predictions of the model” in the Discussion).

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). 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 most of the ECS spikes 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. We note that once the association between CS and fear is acquired, subsequent presentations of CS and US do not weaken or erase it: the interneurons ensure the correct timing and pauses in ECS and F activity, which are conducive for potentiation. Furthermore, since (Krabbe et al., 2019) reported that a fraction of PV interneurons are affected by US, we have also run the simulations for single neuron network with the PV interneuron affected by US instead of CS. In this case as well, all the network realizations are learners (see Fig. S3).

Figure 3.

We show in Fig. 3C results in which F and VIP show an increased activity due to the US for only 15 seconds of the fear conditioning paradigm, thus relaxing the assumption of these neurons being active for the entire forty seconds. We find that learning may still occur under this condition. The reason is that, after 15 seconds of enhanced network activity due to CS and US, 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 15 seconds, 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 restrictive 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).

Mechanisms by which interneurons contribute to potentiation during depression-dominated 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 between ECS and F to lead to potentiation (Fig. S4). In this case, no interneurons are needed (See Discussion “Synaptic plasticity in our model” for the potential necessity of the depression-dominated rule).

In the PV and projection cell 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) 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, although such a low excitation regime establishes the correct fine timing for potentiation, it is not sufficient to lead to potentiation (Fig. 4A, Fig. S2C): 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. S2). Thus, a network without the other interneuron types cannot lead to potentiation. Though a low excitation level for a PV cell is necessary to produce a PING, a higher excitation level is necessary to produce a pause in the ECS and F. This higher excitation level is consistent with the experimental literature showing a strong activation of PV after the onset of CS (Wolff et al., 2014). The higher excitation happens when the VIP cell is silent, whereas a low excitation level is achieved when the VIP cell fires and partially inhibits the PV cell (Fig. 4B, Fig. S2D). The interruption in the ECS and F activity requires the participation of another interneuron, the SOM cell (Figs. 2B, S2): 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. Without these pauses, depression dominates (see SI section “ECS and F activity patterns determine overall potentiation or depression”).

Figure 4.

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 affected by CS and US, respectively, during fear training.

Figure 5.

Fig. 5B shows an example of the network dynamics during fear conditioning. 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, allowing a coordinated disinhibition 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 affected by 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 the active F. Those cells are included for the calculation of the LFP (see below in “Increased low-theta frequency is a biomarker of fear learning”). 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 type is removed from the network; even partial plasticity does not arise.

Increased low-theta frequency is a biomarker of fear learning

We find that fear conditioning leads to an increase in low theta frequency power of the network spiking activity compared to the pre-conditioned level (Fig. 6 A,B); there is no change in the high theta power. We also find that the LFP, modeled as the linear sum of all the AMPA, GABA, NaP-, D-, and H-currents in the network, similarly reveals a low theta power increase when considering the peak of the low theta power, and no significant variation in the high theta power again when considering the peak of the high theta power (Fig. 6 C,D,E). These results are consistent with the experimental findings in (Davis et al., 2017). Specifically, the newly potentiated AMPA synapse from ECS to F ensures F is active after fear conditioning, thus generating strong currents in the PV cells to which it has strong connections (Fig. 6F). It is the AMPA currents to the PV interneurons that are directly responsible for the low theta increase; it is the newly potentiated ECS to F synapse that paces the AMPA currents in the PV interneurons to go at low theta. Thus, the low theta increase is due to added excitation provided by the new learned pathway.

Figure 6.

We find that the AMPA currents are the major contributor to the low theta increase. Although both the AMPA and GABA currents contribute to the power increase in the low theta frequency range (Fig. 6G), the AMPA currents show a dramatic power increase relative to the baseline (the average power ratio of AMPA and GABA post-vs pre-conditioning across 20 network realizations is 3*10³ and 4.6, respectively). Thus, the AMPA currents are the major contributors to the low theta power increase. As a further constrain, the additional unresponsive ECS and F cells in the network were included to ensure we had not biased the LFP towards excitation. Finally, 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. We have also explored another proxy for the LFP (see Supplementary Information and Fig. S5).

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, it 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 20 of them 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 that the low theta power change is a biomarker of successful fear conditioning: it occurs when there is learning and does not occur when there is no learning.

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. The BLA SOM cells do not necessarily have to be the only source of the high theta observed in the BLA during fear learning; the high theta detected in the LFP of the BLA also originates from 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). All forms of STDP require that there be an overlap in the firing of the pre- and post-synaptic cells. In some fear learning paradigms, the US and the CS do not overlap. We address this below under “Assumptions and predictions of the model”, showing how the effects of US and CS on the spiking of the relevant neurons can overlap even in the absence of overlap of US and CS.

There are many STDP rules in the literature (Abbott and Nelson, 2000; Feldman, 2012) but none we are aware of specifically for the amygdala. With the depression-dominated rule used in this study, depression is overall stronger than potentiation, unless there is fine timing among pyramidal cells in the gamma and low theta cycles. With this 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 in the Result section “Mechanisms by which interneurons contribute to potentiation during depression-dominated plasticity”, 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. S4). 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.

With the classical Hebbian plasticity rule, we suggest that learning can occur without the involvement of the VIP and SOM cells. Although fear learning can occur without the depression-dominated rule, we hypothesize that the latter is necessary for other aspects of fear learning and regulation. Generalization of learning can be pathological, and we hypothesize that the modulation created by the involvement of VIP and SOM interneurons is normally used to prevent such overgeneralization. However, in some circumstances, it may be desirable to account for many possible threats, and then a classical Hebbian plasticity rule could be useful. We note that the involvement or not of the VIP-SOM circuit has been implicated when there are multiple strategies for solving a task (Piet et al., 2024). In our situation, the nature of the task (including reward structure) may determine whether the learning rule is depression-dominated and therefore whether the VIP-SOM circuit plays an important role.

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. (See (Cunha-Reis and Caulino-Rocha, 2020) for details about VIP and plasticity). Since high-theta is increased during fear extinction, but not learning, and our model suggests an amygdalar source of high-theta is the SOM cells, we hypothesize that SST peptide may play some role in aiding the formation of networks involved in fear extinction. Also, the CCK peptide has been proposed to promote the switch between fear and safety states after fear extinction (Krabbe et al., 2018).

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 excitatory projection neurons and VIP and PV interneurons show sustained activity during and after the US presentation, thus allowing potentiation through STDP to take place. The medial prefrontal cortex and/or the hippocampus may provide the substrates for the continued firing of the BLA neurons after the 2-second US stimulation. We also discuss below that this network sustained activity may originate from neuromodulator release induced by US (see section “Assumptions and predictions of the model” in the Discussion).

Other brain structures may be involved in other aspects of fear responsiveness, such as fear extinction and prevention of generalization. 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). Also, brain structures such as the prefrontal cortex and hippocampus have been documented to play a crucial role 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. Our model can be extended to include structures in the prefrontal cortex and the hippocampus to further investigate the role of rhythms in the context of discrimination of non-threatening cues and extinction. We hypothesize that a different 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).

Where the rhythms originate, and by what mechanisms

A recent experimental paper (Antonoudiou et al., 2022) suggests that the BLA can intrinsically generate theta oscillations (3-12 Hz) detectable by LFP recordings when inhibition is totally removed due to gabazine application. They draw this conclusion in mice by removing the hippocampus, which can volume conduct to BLA, and noticing that other nearby brain structures did not display any oscillatory activity. In our model, we note that when inhibition is removed, both AMPA and intrinsic currents contribute to the network dynamics and the LFP. Thus, interneurons with their specific intrinsic currents (i.e., D-current in the VIP interneurons, and NaP- and H-currents in SOM interneurons) can indeed affect the model LFP and support the generation of theta and gamma rhythms (Fig. 6G).

Another slice study, (Bratsch-Prince et al., 2024), shows that BLA is intrinsically capable of producing a low theta rhythm with ACh stimulation and without needing external glutamate input. ACh is produced in vivo by the basal forebrain in response to US (Rajebhosale et al., 2024b). Although we did not explicitly include the BF and ACh modulation of BLA in our model, we implicitly include the effect of ACh in BLA by increasing the activity of the VIP cells, which then produce the low theta rhythm. Indeed, low theta in the BLA is known to depend on the muscarinic activation of CCK interneurons, a group of interneurons that overlaps with the class of VIP neurons in our model (Mascagni and McDonald, 2003; Krabbe et al., 2018).

Although the BLA can produce these rhythms, this does not rule out that other brain structures also produce the same rhythms through different mechanisms, and these can be transmitted to the BLA. Specifically, it is known that the olfactory bulb produces and transmits the respiratory-related low theta (4 Hz) oscillations to the dorsomedial prefrontal cortex, where it organizes neural activity (Bagur et al., 2021). Thus, the respiratory-related low theta may be captured by BLA LFP because of volume conduction or through BLA extensive communications with the prefrontal cortex. Furthermore, high theta oscillations are known to be produced by the hippocampus during various brain functions and behavioral states, including during spatial exploration (Vanderwolf, 1969) and memory formation/retrieval (Raghavachari et al., 2001), which are both involved in fear conditioning. Similarly to the low theta rhythm, the hippocampal high theta can manifest in the BLA. It remains to understand how these other rhythms may interact with the ones described in our paper. However, we emphasize that there is also evidence (as discussed above) that these rhythms arise within the BLA.

Assumptions and predictions of the model

The interneuron descriptions in the model were constrained by the electrophysiological properties reported in response to hyperpolarizing currents (Sosulina et al., 2010). Specifically, we modeled the three subtypes of VIP, SOM, and PV interneurons displaying bursting behavior, regular spiking with early spike-frequency adaptation, and regular spiking without spike-frequency adaptation, respectively. Focusing on VIP interneurons, we were able to model the bursting behavior by including the D-type potassium current. This current is thought to exist in the VIP interneurons in the cortex (Porter et al., 1998), but whether this current is also found in the BLA is still unknown. We modeled the SOM interneurons with NaP- and H-currents, similar to the OLM cells in the hippocampus. Due to these currents, the VIP and SOM cells are able to produce low- and high-theta oscillations, respectively. The presence of these currents, and the neurons’ ability to exhibit oscillations in the theta range during fear conditioning and at baseline in BLA, which are assumptions of our model, should be tested experimentally. Our model predicts that blockade of D-current in VIP interneurons (or silencing VIP interneurons) will both diminish low theta and prevent fear learning. Finally, the model assumes the absence of significantly strong connections from the excitatory projection cells ECS to PV interneurons, unlike the ones from F to PV. Including those synapses would alter the PING rhythm created by the interactions between F and PV, which is crucial for fine timing between ECS and F needed for LTP.

Our model, which is a first effort towards a biophysically detailed description of the BLA rhythms and their functions, does not include the neuron morphology, many other cell types, conductances, and connections that are known to exist in the BLA; models such as ours are often called “minimal models” and constitute most biologically detailed models. For example, although there is considerable variability in the activity patterns of both VIP cells and SOM cells (Sosulina et al., 2010; Guthman et al., 2020; Ünal et al., 2020; Vereczki et al., 2021), our focus was specifically on those subtypes that generate critical rhythms within the BLA. Such minimal models are used to maximize the insight that can be gained by omitting details whose influence on the answers to the questions addressed in the model are believed not to be qualitatively important. We note that the absence of these omitted features constitutes hypotheses of the model: we hypothesize that the absence of these features does not materially affect the conclusions of the model about the questions we are investigating. Of course, such hypotheses can be refuted by further work showing the importance of some omitted features for these questions and may be critical for other questions. Our results hold when there is some degree of heterogeneity of cells of the same type, showing that homogeneity is not a necessary condition.

Our study suggests that all the interneurons modeled are necessary for associative learning provided that the STDP rule is depression-dominated. This prediction could be tested experimentally by selectively silencing each interneuron subtype in the BLA: if the associative learning is hampered by silencing any of the interneuron subtypes, this supports our study. Finally, the model prediction could be tested indirectly by acquiring more information about the plasticity rule involved in the BLA during associative learning. We found that all the interneurons are necessary to establish fear learning only in the case of a depression-dominated rule. This rule ensures that fine timing and pauses are always required for potentiation: interneurons provide both fine timing and pauses to pyramidal cells, making them crucial components of the fear circuit.

Finally, our model requires the effect of the CS and US inputs on the BLA neuron activity to overlap in time in order to instantiate fear learning through STDP. Such a hypothesis, that learning uses spike-timing-dependent plasticity, is common in the modeling literature (Bi and Poo, 2001; Caporale and Dan, 2008; Markram et al., 2011). Current paradigms of fear conditioning include examples in which the CS and US stimuli do not overlap (Krabbe et al., 2019). Such a condition might seem to rule out the mechanisms in our paper. Nevertheless, the argument below suggests that the effects of the CS and US can cause an overlap in neuronal spiking of ECS, F, VIP, and SOM, even when CS and US inputs do not overlap.

Experimental recordings cannot speak to the rate of spiking of BLA neurons during US due to recording interference from the shock. However, evidence suggests that ECS activity should increase during the US due to the release of acetylcholine (ACh) from neurons in the basal forebrain (BF) (Rajebhosale et al., 2024b). Pyramidal cells of the BLA robustly express M1 muscarinic ACh receptors (McDonald and Mott, 2021). Thus, ACh from BF should elicit a depolarization in pyramidal cells. Indeed, the pairing of ACh with even low levels of spiking of BLA neurons results in a membrane depolarization that can last 7 – 10 s (Unal et al., 2015). Other modulators, including dopamine, may also play a role in producing the sustained activity. Activation of US leads to increased dopamine release in the BLA (Harmer and Phillips, 1999; Suzuki et al., 2002). D1 receptors are known to increase the membrane excitability of BLA projection neurons by lowering their spiking threshold (Kröner et al., 2005). Thus, neuromodulator release should induce higher spiking rates and more sustained activity in the ECS and F neurons during and after the presentation of US, thus ensuring a concomitant activation of ECS and fear (F) neurons necessary for STDP to take place. Thus, the activation of the US can lead to continued and higher firing rates of ECS and F. The effect of dopamine can last up to 20 minutes (Kröner et al., 2005). For CS-positive neurons, the ACh modulation coming from the firing of US may lead to a temporary extension of firing that is then amplified and continued by dopaminergic effects.

Hence, we suggest that a solution to the problem apparently posed by the non-overlap US and CS in some paradigms of auditory fear conditioning (Krabbe et al., 2019) may be solved by considering the roles of ACh and dopamine in the BLA. The model we have may be considered a “minimal” model that puts in by hand the overlap in activity due to the neuromodulation without explicitly modeling it. We have used the simplest way to model overlap without assumptions about timing specificity in the overlap. We note that, even though ECS and F neurons have the ability to fire continuously when ACh and dopamine are involved, the participation of the interneurons enforces periodic silence needed for the depression-dominated STDP.

Is STDP needed in fear conditioning?

The study in (Grewe et al., 2017) questions the validity of the Hebbian model in establishing associative learning during fear conditioning. There are several critiques we discuss here. The first critique is that Hebbian plasticity does not explain the experimental finding showing that both upregulation and downregulation of stimulus-evoked responses are present between coactive neurons. The upregulation is provided by our model, so the issue is the downregulation, which is not addressed by our model. However, our model highlights that coactivity alone does not create potentiation; the fine timing of the pre- and postsynaptic spikes determines whether there is potentiation or depression. Here, we find that PING networks are instrumental in setting up the fine timing for potentiation. We suggest that networks not connected to produce the PING may undergo depression when coactive.

The second critique raised by (Grewe et al., 2017) is that Hebbian plasticity alone does not explain why most of the cells exhibiting enhanced responses to the CS did not react to the US before fear conditioning. They suggest that neuromodulators may provide a third condition (besides the activity of the pre- and postsynaptic neurons) that changes the plasticity rule. Our model also does not explicitly address this experimental finding since it requires F to be initially activated by US in order for the fear association to be established. We agree that the fear cells described in (Grewe et al., 2017) may be depolarized by the US without reaching the spiking threshold; however, with neuromodulation provided during the fear training, the same input can lead to spiking, enabling the conditions for Hebbian plasticity. Our discussions above about how neuromodulators affect excitability are relevant to this point. We do not exclude that other forms of plasticity may play a role during fear conditioning in cells not initially activated by the US, but this is not the topic of our modeling study.

The third critique raised by (Grewe et al., 2017) is that Hebbian plasticity cannot explain why the majority of cells that were US- and CS-responsive before training have a reduced CS-evoked response afterward. The reduced response happens over multiple exposures of CS without US; this can involve processes similar to those present in fear extinction, which require plasticity in further networks, especially involving the infralimbic cortex (Milad and Quirk, 2002; Burgos-Robles et al., 2007). An extension of our model could investigate such mechanisms. In the fourth critique, (Grewe et al., 2017) suggests that the Hebbian plasticity rule cannot easily account for the reduction of the responses of many CS+-responsive cells, but not of the CS−-responsive cells. We suggest that the circuits involving paradigms similar to fear extinction do not involve the CS-cells.

Overall, we agree with (Grewe et al., 2017) that neuromodulators play a crucial role in fear conditioning, especially in prolonging the US- and CS-encoding activity as discussed in (see section “Assumptions and predictions of the model” in the Discussion), or even participating in changing the details of the plasticity rule. A possible follow-up of our work involves investigating how fear ensembles form and modify through fear conditioning and later stages. This follow-up work may involve using a tri-conditional rule, as suggested in (Grewe et al., 2017), in which the potential role of neuromodulators is taken into account in the plasticity rule in addition to the pre- and postsynaptic neuron activity. Another direction is to investigate a possible relationship between neuromodulation and a depression-dominated Hebbian rule.

Comparison with other models

Many computational models that study fear conditioning have been proposed in the last years; the list includes biophysically detailed models (e.g., (Li et al., 2009; Kim et al., 2013), firing rate models (e.g., (Krasne et al., 2011; Vlachos et al., 2011; Ball et al., 2012)), and connectionist models (e.g., (Edeline and Weinberger, 1992; Armony et al., 1997; Moustafa et al., 2013) (for a review see (Nair et al., 2016)). Both firing rate models and connectionist models use an abstract description of the interacting neurons or regions. The omission of biophysical details prevents such models from addressing questions concerning the roles of dynamics and biophysical details in fear conditioning, which is the aim of our model. There are also biophysically detailed models (Li et al., 2009; Kim et al., 2013, 2016; Feng et al., 2019), which differ from ours in both the physiology included in the model and the description of how plastic changes take place. One main difference in the physiology is that we differentiated among types of interneurons, since the fine timing produced for the latter was key to our use of rhythms to produce spike-timing-dependent plasticity. The origin of the gamma rhythm (but not the other rhythms) was investigated in (Feng et al., 2019), but none of these papers connected the rhythms to plasticity.

The most interesting difference between our work and that in (Li et al., 2009; Kim et al., 2013, 2016) is the modeling of plasticity. We use spike-timing-dependent plasticity rules. The models in (Li et al., 2009; Kim et al., 2013, 2016) were more mechanistic about how the plasticity takes place, starting with the known involvement of calcium with plasticity. Using a hypothesis about back propagation of spikes, the set of papers together develop a theory that is consistent with STDP and other instantiations of plasticity (Shouval et al., 2002a, 2002b). For the purposes of our paper, this level of detail, though very interesting, was not necessary for our conclusions. By contrast, in order for the rhythms and the interneurons to have the dynamic roles they play in the model, we needed to restrict our STDP rule to ones that are depression-dominated. Our reading of (Shouval et al., 2002b) suggests to us that such subrules are possible outcomes of the general theory. Thus, there is no contradiction between the models, just a difference in focus; our focus was on the importance of the much-documented rhythms (Seidenbecher et al., 2003; Courtin et al., 2014a; Stujenske et al., 2014; Davis et al., 2017) in providing the correct spike timing. We showed in the Supplementary Information (“Classical Hebbian plasticity rule, unlike the depression-dominated one, shows potentiation even with no strict pre and postsynaptic spike timing”) that if the STDP rule was not depression dominated, the rhythms need not be necessary. We hypothesize that the necessity of strict timing enforced by the depression-dominated rule may foster the most appropriate association with fear at the expense of less relevant associations.

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

Limitations and caveats

LFPs recorded in the experiments are thought to be mainly created by transmembrane currents in neurons located around the electrode and depend on several factors, including the morphology of the arborization of contributing neurons and the location of AMPA and GABA boutons (Katzner et al., 2009; Lindén et al., 2011; Łęski et al., 2013; Mazzoni et al., 2015). Since our model has no spatial extension, we used an LFP proxy; this proxy was shown to reflect the rhythmic output of the network, which we believe to be the essential result (for more details see Results “Increased low-theta frequency is a biomarker of fear learning”, and Supplementary Information “A higher low theta power increase emerges in LFP approximated with the sum of the absolute values of the currents compared to their linear sum”).

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 of each type, 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

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.

Additional information

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 NK; 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 NaP 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). All the currents are expressed in units per area, rather than absolute units, to avoid making assumptions about the size of the neuron surface.

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_L) 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). 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⁹. To make Fig. S3, we also considered a variation of the model with PV interneurons affected by US, instead of CS, as reported in (Krabbe et al. 2019).

Network connectivity and synaptic currents

We modeled the network connectivity as presented in Fig. 2B, derived from the prominent functional, instead of structural, connections reported in (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, as common in the modeling literature (Jensen et al., 2005; Traub et al., 2005; Chartove et al., 2020). 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. S2 and Fig. S4). The update is described as follows. At each postsynaptic neuron (F) spike, M(t) is decremented by an amount A_N = 0.005, i.e., M(t) = M(t) − A_N. At each presynaptic neuron (ECS) spike, P(t) is incremented by an amount A_O = 0.005, i.e., P(t) = P(t) + A_O. This Hebbian plasticity rule is depression-dominant when A_O = A_{_} and τ_O < τ_N, as formulated in this work. The exponential decay is described by the following equations:

with τ_N = 28 ms and τ_O = 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. S2A, 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_O = 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. We considered as an LFP proxy as the linear sum of all the AMPA, GABA, NaP-, D-, and H-currents in the network. The D-current is in the VIP interneurons, and NaP-current and H-current are 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)