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). A recent rodent experimental study (Antonoudiou et al., 2022) suggests 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 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 potential roles of rhythms 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 (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 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 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). 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). 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 some heterogeneity in the last two sections of the Results.

Figure 2.

Fig. 2A (bottom) shows a typical dynamic of the network before and after the US input onset, with US modeled as a Poisson spike train at ∼50 Hz; the network produces all the rhythms originating from the interneurons alone or through their interactions with the excitatory projection neurons (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 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. In our simulations, the 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. 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. CS is described as a Poisson spike train at ∼50Hz, independent of the US input. 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. S2 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 most of the 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. 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. Finally, 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 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).

We note that the CS and US inputs modeled as independent Poisson spike trains represent stimuli with no structure. Although we have not explicitly modeled pulsating pips, as common in auditory fear conditioning (e.g., (Stujenske et al., 2014; Krabbe et al., 2019)), we show in Fig. S4 that potentiation can be achieved over a relatively wide range of gamma frequencies. This indicates that overall potentiation is ensured if the gamma frequency transiently increases after the pip.

Physiology of interneuron types is critical to their role in 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 to lead to potentiation (Fig. S5). In this case, no interneurons are needed (See Discussion “Synaptic plasticity in our model” for the potential necessity of the depression-dominated rule).

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

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 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 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 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 show 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 and no significant variation in the high theta power (Fig. 6 C,D,E). These results reproduce the experimental findings in (Davis et al., 2017), and Fig. 6 F,G show that the low theta increase is due to added excitation provided by the new learned pathway. The additional unresponsive ECS and F cells in the network were included to ensure we had not biased the LFP towards excitation. Nevertheless, although both the AMPA and GABA currents contribute to the power increase in the low theta frequency range (Fig. 6F), 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). This points to the AMPA currents as the major contributor to the low theta power increase. 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. 6G). 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. S6).

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, 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 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 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). 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. S5). 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 show 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 (Cuhna-Reis et al., 2020) for details about VIP and plasticity). We hypothesize that SST peptide is critical in fear extinction and preventing fear generalization, but not to initial fear learning. 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 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.

Other brain structures may be involved in later stages 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). Brain structures such as the 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. 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).

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 come up with 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).

The study in (Grewe et al., 2017) suggests that associative learning resulting from 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, we show how potentiation between coactive cells occurs when strict pre-post spike timing and appropriate pauses in the spiking activity arise. Depression happens when one or both of these components are not present. Thus, in our model, depression represents the absence of successful fear association and does not take part in the reshaping of the ensemble encoding the association, as instead suggested in (Grewe et al., 2017). 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 addition to the pre- and postsynaptic neuron activity; this may lead to both potentiation and depression in establishing an associative memory.

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 under certain conditions, such as reduced inhibitory tone. 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. Our model also supports the idea that intrinsic mechanisms in the BLA can support the generation of the low theta, high theta, and gamma rhythms.

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.

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. Similarly, we endowed the SOM interneurons with NaP- and H-currents, as the OLM cells in the hippocampus. Due to these currents, the VIP and SOM cells are able to show 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, 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 the majority of biologically detailed models. 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 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 validates 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.

The modeling of the interneurons assumes the involvement of various intrinsic currents; the inclusion of those currents can be considered hypotheses of the model. 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.

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. Despite paradigms involving both overlapping (delay conditioning, where US co-terminates with CS (e.g., (Lindquist et al., 2004)), or immediately follows CS (e.g., (Krabbe et al., 2019)) and non-overlapping (trace conditioning) CS/US inputs exist, we hypothesized that concomitant activity in CS- and US-encoding neuron activity should be crucial in both cases. This may be mediated by the memory effect due to metabotropic effects (Whittington et al., 1995) as suggested above, or by the contribution from other brain regions (see section “Involvement of other brain structures” in the Discussion). The fact that plasticity occurs with US memory trace is a consequence of our larger hypothesis that fear learning uses spike-timing-dependent plasticity; such a hypothesis about plasticity is common in the modeling literature.

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, 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 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 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). 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 _K = 80 mSKcm² 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). 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. S2B 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). Finally, for Fig.S4, we considered a range of frequencies for the CS stimulus. To generate the three Poisson spike trains with average frequencies from 48 to 64 Hz in Fig. S4, we set λ = 800, 1000, 1200.

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; Kumar et al., 2011; 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. 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₋. 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. 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)