Introduction

The balance between excitatory and inhibitory drives (EI balance) in neural networks has a wide range of computational functions. These include increased storage capacity¹, optimal coding ², fast temporal encoding and gating ^3–5, and higher dynamic range ^6,7. EI balance applies over tight time scales and even combinations of inputs at a single neuron level, termed precise balance⁶. However, little is known about how EI balance itself evolves dynamically due to activity-driven plasticity in sparsely active networks.

Active synapses undergo STP, during which neurotransmitter release varies upon receiving successive action potentials at the boutons ^8–11. STP emerges from the interaction and competition between pre-synaptic processes of vesicle depletion, calcium influx, calcium build-up, and neurotransmitter reuptake. Together, these processes create a temporal filter between activity in presynaptic and postsynaptic neurons and networks ^8,12,13. Excitatory and inhibitory synapses show distinct short-term modifications in their synaptic weights ^13–15.

Activity in the CA3 is believed to converge to ensembles of cells representing distinct stimuli and contexts ^16,17. In abstract terms, these ensembles form instantaneous patterns, or vectors of activity, which are decoded in a heteroassociative manner by the CA1¹⁸. At a cellular physiology level, this static view implies differential responses by individual CA1 neurons to different patterns. However, precise EI balance in the feedforward excitatory-inhibitory (FFEI) network of CA3-CA1 normalizes these differences⁶. Given that activity is not static, it may be more salient for neurons to detect transitions between input patterns. Differential changes on E and I synapses may provide a mechanism for single neurons to extract interesting features of spatiotemporal inputs through STP¹⁹.

In the current study, we used in vitro whole cell patch clamp recordings and patterned optical stimulation of the CA3-CA1 in the mouse hippocampus, to probe how EI balance evolves over spatiotemporal input trains at CA3. We incorporated this STP data into mass-action stochastic models of presynaptic release for E and I synapses, and embedded a CA1 neuron receiving hundreds of such synapses into an abstract model of the hippocampal FFEI network. This model predicted that the postsynaptic CA1 neuron could detect pattern changes, which we confirmed experimentally. Finally we showed that change detection is enhanced when gamma frequency inputs trigger CA1 neuron spiking.

Results

Optogenetic patterned stimulation in the CA3-CA1 network

We obtained whole-cell patch-clamp recordings from CA1 pyramidal cells (PCs) while stimulating channelrhodopsin (ChR2) expressing CA3 pyramidal cells using optical patterns generated by a digital micromirror device (Fig 1A-B, see Methods). A pattern consisted of either 1, 5, or 15 spots (spot size 13 x 7 µm, power 14.5 µW/spot) of blue light (470 nm). We illuminated the CA3 layer in a series of eight pulses at a fixed frequency of 20, 30, 40, and 50 Hz (Fig 1D). A probe pulse preceded this train by 300 ms to provide a reference baseline response. We measured the dynamics of E-I balance in a three-dimensional stimulus space consisting of CA3 activation size (number of spots), activation frequency (stimulation frequency), and short-term plasticity dependence (pulse number) in the train (Fig 1C).

Fig 1:

To monitor the strength of the total resultant optogenetic activation of the CA3 layer (Extended Data Fig 1), we used an extracellular field electrode in the CA3 stratum radiatum (Fig 1A, methods). The field response correlated well with optically-driven CA3 PC depolarization (Fig 1E-G), and scaled with the size of the pattern (Fig 1F). The field response showed a slight decrease over the course of the pulse train of approximately 2% per pulse (regression fit slope=0.02, r²=0.05). We attribute this to ChR2 desensitization.

As expected, patterns with more spots elicited larger responses in CA1 in current clamp (Fig 1 H-I) and voltage clamp (Fig 1 K-L). We separated E and I components by holding the cell at ࢤ70 mV (GABAR1 reversal potential) and 0 mV (gluR reversal potential) respectively. As reported previously⁶, different activation patterns in CA3 produced proportional E and I post-synaptic currents in the recorded neuron (Fig 1K). However, over the pulse train, E and I underwent distinct STP profiles (Fig 1 M).

Thus our patterned optical stimulation provided multiple input combinations to probe the dynamics of monosynaptic excitation and disynaptic inhibition arriving at CA1 pyramidal cells.

EI balance evolves over the pulse train

We next asked how precise EI balance at the level of CA1 PCs dynamically evolves due to short-term plasticity processes across a range of frequencies and input strengths (Fig 2A). All the PSPs of an 8-pulse train were normalised to the probe pulse. The PSP varied significantly on all the tested three axes of parameters: pattern size (p<0.001), stimulus frequency (p<0.001), and pulse index (p<0.001) (ANOVA). Notably, the PSP of the first four pulses was significantly larger than that in the last four. We tested this by comparing the pulse responses divided between the two halves of the train against a shuffled pulse order (p<0.001, Wilcoxon rank test).

Fig 2.

While PSPs are a measure of subthreshold EI balance dynamics, neuronal firing is the network readout of escape from EI balance. Spiking was rare for the stimulus strengths in our experiments. We observed spiking in only 7% of trials (n=154/2201). We measured the likelihood of escape from EI balance by counting the fraction of total spiking trials that elicited for each combination of parameters (Fig 2B). Similar to the PSP response amplitudes, spike likelihood was higher in the first half of the pulse train (p < 0.001, ANOVA) Similarly, the spike likelihood varied significantly across the frequencies and number of spots (stimulus frequency p<<0.001, number of spots, pulse index p << 0.001, ANOVA).

We next measured E and I currents under voltage clamp across pattern size, pulse frequency, and pulse number. The EPSC showed a trend of early potentiation followed by depression (Fig 2 Ei-Eii), while the inhibition underwent depression from the start (Fig 2 Fi-Fii). This effect was significant at p < 10^-6 across all the tested dimensions (ANOVA).

Thus EI balance over a burst briefly tilts in favour of E, after which both E and I undergo depression. This may result in a brief opportunity of escape from EI balance two or three pulses into a burst.

Multi-synapse summation and divisive normalization evolve over a pulse train

How might a neuron sum distinct input ensembles if they have distinct histories? We used the pattern-size axis of our stimulus protocol to probe history dependence of summation and subthreshold divisive normalization. We have previously shown that summed PSPs are not direct outcomes of EI ratios, because the delay between E and I onset is also a function of stimulus strength⁶. We therefore measured each of these terms: postsynaptic potentials, postsynaptic excitatory and inhibitory currents, and their onset and peak delays (Fig 3 A-D, extended data figure 2) over the pulse series, and obtained gamma (γ) as a measure of linearity ⁶ (Fig 3 E-G). Gamma (γ) relates the observed response of a CA1 cell to the expected response obtained by summing the inputs linearly in the following manner:

The smaller the gamma, the larger the deviation from the linear combination of inputs due to the divisive effects of inhibition.

Figure 3:

Using responses from the reference pulse, we replicated earlier observations^6,7 showing divisive normalisation, and obtained a median gamma of 7.16 (95% CI = 4.76 - 10.2)(Fig 3G). For comparison, we also performed a linear regression and obtained a median slope of 0.57 (95% CI = 0.45 - 0.72). We found that both gamma and slope values show a rise and then a fall in a history-dependent manner (Fig 3G). Thus the differential STP profiles for E and I synapses yield a history-dependence of pattern-specific summation, with greater linearity around pulse two or three in a train.

Presynaptic release models are bounded by observations from pulse-train sequences

Having characterised the dynamics of plasticity of excitatory and inhibitory synaptic currents over trains of pulses, we next used this data to develop models of presynaptic signalling. Datasets with similar pulse trains have been used to obtain phenomenological models (e.g.,²⁰). Instead, we chose to build a chemical kinetics-based multi-step neurotransmitter vesicle release model similar to previously published models ^21,22 (Figure 4A). This incorporated both calcium buffering, which contributes to short-term potentiation, and depletion of docked vesicles, which underlies short-term depression. We implemented a ball-and-stick cellular model having AMPA and GABA synapses on the dendrite, with current or voltage-clamp recordings from the soma to compare to experiments. This chemical-electrical model was implemented in MOOSE ²³, and parameters were fit using a multi-objective optimization pipeline to match waveforms and dynamics over all frequencies and the entire pulse train (Methods, Extended Fig 3) ²⁴. Each cell had somewhat distinct STP characteristics (Supplementary Fig 1). We independently fit presynaptic signaling models for E and I to their respective voltage-clamp waveforms (methods, Supplementary Fig 2,3).

Figure 4:

While the E fit applied directly to the CA3->CA1 pyramidal neuron glutamatergic synapse, the I fit was a composite of CA3->Interneuron and Interneuron->CA1 projections (Fig 4B, C). Based on the peak and onset times, we made the further simplifying assumption that only one interneuron class was involved, most likely PV interneurons. An example of CA3-triggered presynaptic signalling and postsynaptic currents is presented in Extended Fig 4. Note that the synaptic chemistry was modelled using the Gillespie Stochastic Systems Algorithm ²⁵, hence the responses are noisy and synaptic release is probabilistic.

We used parameters for Cell 7492 for subsequent model-building (Extended Fig 3).

Burst EPSP responses constrain single-cell parameters

We next embedded our presynaptic signalling models as synaptic boutons in a conductance-based electrophysiology model of a neuron, and used burst-EPSP responses to constrain its properties. 100 glu synapses were placed on dendritic spines, and 200 GABA synapses on the dendrite (Figure 4B). In the CA1 the E:I synapse ratio is between 12:1 to 25:1 ²⁶, but we implemented a subset of E synapses targeted by the CA3 pattern-responsive neurons in the interests of computational efficiency, as most of the E inputs would be silent.

We implemented an abstract hippocampal network to drive activity on the synapses. CA3 was implemented as a 16x16 array of integrate-and-fire neurons. Optical patterned input was mapped onto this array (Extended Fig 5) and included charge buildup from successive light pulses as well as ChR2 desensitisation (methods, Extended Fig 6). The interneuron layer was a 16x16 array of binary neurons. The CA3 array projected sparsely to the 100 glutamatergic synapses on the CA1 model and to the 256 interneurons. The interneuron array projected onto the 200 GABAergic synapses (Fig 4C).

As a first-pass test of our model and comparison with experimental data, we modelled pulse train experiments performed using current clamp. We ran the simulations deterministically to correspond to the averaged experimental traces (Fig 5 A-D). Since simulations permit access to internal state variables not available to the experiments, we also looked at how E and I currents and conductances balanced out during the pulse train (Figure 5F-O). Notably, the simulations predict that the E current builds up, especially on the second pulse of the train. This corresponds well with the observations in Figure 5 A-D and in Figures 2 and 3.

Figure 5:

Overall, the synapse and single neuron components of our model worked well together to replicate burst inputs in current clamp.

Poisson train stimuli constrain network parameters

Following the synapse and cell-level analysis of signal summation and short-term plasticity in the CA3-CA1 system, we next tuned our model to network-level readouts. Here we explicitly included an analysis of scatter in the experimental data. We modeled the origin of this scatter as stochastic synaptic chemical kinetics and transmitter release (methods). We delivered a ‘frozen’ Poisson-train sequence to experiment and model, to test the capacity of the model to replicate data and to further constrain network-level parameters using the statistics of the responses. Figure 6 A and B show example EPSP traces for a single trial of an experiment and simulation, respectively. We used simple parameter sweep runs to tune the model network parameters.

Figure 6:

We first compared the probability that each optical stimulus would elicit an EPSP (Fig 6 C and D). As expected, 15-square patterns (yellow dots) almost always gave an EPSP, while 5-square patterns frequently failed. In the simulation there was no instrumentation noise, but we found that setting the detection threshold for an EPSP to 0.5 mV yielded similar scatter plots to experiment. In Fig 6 E and F we compare scatter plots of EPSP amplitudes over time. Note that the initial high EPSP responses declined over the first two seconds (Fig 6E). We asked if ChR2 desensitisation of the stimulus might account for this. We found that the field response does indeed show a decline over multiple pulses (Fig 1, Extended Fig 7). We also checked for STD as a mechanism by comparing the decline between the reference and the no-STP model (Extended Fig 8). This suggests that there is a small additional contribution due to presynaptic depression.

We then plotted each EPSP against the inter-stimulus interval (ISI) preceding it (Fig 6 G, H). We observed that the EPSP was high for very short intervals and then declined to a lower steady value, both in experiments and simulations. We initially assumed this was due to increased STP for closely succeeding input pulses. Unexpectedly, the no-STP model also showed the same profile of high EPSP at short ISIs (Extended Fig 9). Removal of NMDA receptors in the model also did not change the EPSP-vs-ISI profile (Extended Figure 9C). Thus, we conclude that the high EPSP was simply due to charge accumulation during rapid synaptic input.

We then compared the distributions of EPSP amplitudes (Figure 6 I,J). While the large sample size meant that the histograms were statistically different at p<1e-50, the K-S test statistics were small (<0.18) indicating that the qualitative profiles were similar. As a further comparison between experiment and model, we estimated the power spectral density of the EPSP waveforms (Figure 7 K, L). These correlated very closely between experiment and simulation (Spearman’s coeff. = 0.92 and 0.94 for 5 and 15 square, p<1e-80 for both).

Figure 7:

Overall, we were able to replicate numerous features of the experimental dataset in a multiscale model incorporating presynaptic signalling, postsynaptic electrophysiology, and abstracted network connectivity and responses. While there was considerable scatter in both datasets, we obtained semi-quantitative matches between experiment and model for several readouts.

STP and sparse coding yield pattern mismatch detection

We next used our model to predict a single-cell computation emerging from patterned input and STP. As a heuristic argument, consider two patterns of activity, A and B, comprising different ensembles of CA3 pyramidal cells. Let these overlap to a small extent, say AՈB < 20% (Figure 7A). Since the CA3 projects directly to CA1 spines, the overlap onto excitatory synapses will also be around AՈB. We ignore inhibitory effects for simplicity, but these are considered in the detailed simulations below. Now, we consider STP. Repeated inputs of pattern A result in the initial strengthening of excitatory inputs, then depression (Figure 7B). At the transition from pattern A to B, most of the excitatory inputs have not been stimulated before, so the fresh synapses are stronger. This gives a larger excitatory drive (Figure 7B rightmost cell). From this heuristic argument, we predict two things: First, at transitions between patterns (mismatches), we expect a transient elevated response in the postsynaptic cell. Second, the size of the transient should be larger if the input patterns are sparse in E, because this lessens the overlap between them.

To test this in the model, we simulated a sequence of patterns of the form “AAAAAAAABBBBBBBBCCCCCCCCDDDDDDDD”, that is, eight repeats of each of the four patterns. To calibrate the basal response, we ran the model without STP and without stochasticity (Figure 7C). We observed that the output EPSP amplitudes were similar for the successive input patterns. When we introduced STP, each of the transitions generated a positive transient (Figure 7D), even if the basal response led to a small downward step. With denser connectivity, the transients were smaller (Figure 7E).

Based on these predictions, we asked if we could experimentally detect pattern transitions in our slice preparation. We delivered the same 8x4 stimulus sequence at repeat rates of 8, 20 and 50 Hz to span a range from theta to gamma frequencies (Figure 7 F). We did this for patterns with 5 and 15 squares, respectively. Figure 7 Fiv shows a strong response to the third transition at 50 Hz for the 15-square case. We then asked if our model with full stochastic synaptic kinetics could replicate these results. Our sparse connectivity model produced multiple strong mismatch responses, consistent with sparse overlap of synapses (extended Fig 9), and these were reduced in a denser model (Fig 7 G, H, Supplementary Movie 1). This agrees with our heuristic analysis above. Like in the experiments, the 50 Hz pattern, in particular, showed strong mismatch responses (Figure 7 G, H). We used the Wilcoxon test to compare experimental EPSPs immediately before and after pattern transitions (Fig 7I-K). More mismatch responses were seen for 15 square patterns and gamma-frequency pulses (Fig 7 K).

We highlight the observation that mismatch detection was selective for specific transitions. This occurred both in experiments (e.g, Fig 7 Fiv for the third transition) and in simulations (e.g., Fig 7 Giv for transitions one and three).

We then explored a range of connectivity parameters in our model to predict how mismatch computations might be affected (Table 1, Figure 7L-Q).

Table 1:

We first confirmed that sparse patterns had stronger mismatch responses (Figure 7L). There was a sharp peak in selectivity for low values of wtGlu, and a plateau for higher values (Fig 7M). Mismatch was small for low values of wtGABA, supporting the role of inhibition in mismatch detection (Fig 7N). Greater values of excitatory synaptic projection probability, pCA3_CA1 led to lower selectivity (Fig 7O). This is complementary to the trend for sparseness (Fig 7L), which is expected since they reciprocally affect the number of activated synapses for each stimulus pattern. Both the inhibitory projection terms pCA3_Inter and pInter_CA1 showed increasing selectivity with increasing connectivity (Fig 7P,Q). We interpret this as due to increasing overlap between input patterns, so that depression of inhibitory inputs carries over between patterns.

Thus, at this stage we had predicted STP-dependent single-neuron detection of mismatches in sequences of inputs, confirmed it experimentally, and then explored network conditions which supported it. Notably, sparse connectivity was a prerequisite.

Mismatch detection is effective in spiking neurons

Finally, we utilized the model to predict how mismatch responses might manifest in spiking neurons. To do this we introduced voltage-dependent Na and K_DR channels in the soma of the CA1 neuronal model and increased the wtGlu and wtGABA each by 67%. We obtained spiking responses for 50 stochastic runs with patterned input at 50 Hz (Figure 8 A,B), and averaged the resulting firing rates with a 10ms window, as a proxy for population activity (Figure 8 C-K). To our surprise, the spiking model yielded strong mismatch detection (Figure 8C). Furthermore, spiking could track each optical pulse. As a consistency test, we replicated the finding from Figure 7 that high sparsity of connections impaired selectivity (Figure 8D, E).

Figure 8.

We next explored a range of stimulus conditions for functional implications. First, we considered the frequently used behavioural paradigm of ‘oddball’ detection, which produces a rapid ∼100-250ms evoked potential in several modalities ^27,28. We delivered a constant, repeating pattern interspersed every 8 stimuli with the oddball (deviant) stimulus. Under these conditions the mean spiking approached zero, but the oddballs stood out with brief, brisk spiking (Figure 8 F).

Gap stimuli, where the deviant stimulus is an absence of stimulus, are also known to be strong behavioural triggers for oddball responses (refs). However, in our model the gap stimulus sequence led to nearly zero spiking (Fig 8 G), suggesting that a more complete CA1 connectivity network model may be required to detect the gap.

We then asked how differential STP in E and I might impact mismatch detection. We first removed STP from glutamatergic synapses (Fig 8 H), then from GABA synapses (Fig 8 I), then from both (Fig 8 J). Remarkably, mismatch detection was impaired in all three cases. Lack of Glu STP, either alone or with GABA STP, eliminated the primary mechanism of mismatch detection through depression of repeated excitatory inputs. Lack of just GABA STP kept the inhibitory tone high, while the Glu synapses became weaker. Hence the response almost shut down though some selectivity remained (Fig 8 I). When we removed GABA input entirely we observed elevated spiking but selectivity was retained for the first transition (Figure 8 K).

We next performed two controls to confirm the role of transitions in eliciting mismatch responses. In the first control we delivered the same spatial pattern for all 32 pulses (Figure 8J). As expected, this led to a rapid decline in spiking. Next, we randomly selected one pattern among our set of five for each pulse (Fig 8 K). Surprisingly, this also gave a decline in spiking, though it was not as steep as the uniform pattern. Upon further analysis, it turned out that random sequences had sufficient overlap to trigger a subset of the same synapses (Extended Fig 10). The net effect was that excitatory synapses were stimulated at ∼21 Hz for random patterns, as compared to the 50 Hz for uniform patterns. This led to slower synaptic depression.

How might the mismatch response manifest in more natural stimulus contexts than a train of 32 pulses? To investigate this, we considered the well-known phenomenon of gamma bursts modulated by a theta background²⁹. To set the gamma burst parameters, we confirmed that mismatch selectivity persisted broadly over the gamma frequency range (50-120 Hz, Figure 8 N). We then provided theta-modulated gamma bursts, where each theta cycle had five pulses of input at 100Hz, with the start of each burst separated by 130 ms to give a 7.69 Hz theta rhythm (Fig 8 O). We compared the case where each burst was a different pattern with the control where each burst was the same pattern (Fig8 P, Q). The CA1 neuronal response was much larger for the case with changing input patterns between theta cycles. Thus mismatch detection also works when the ensemble of neurons comprising a gamma burst changes between theta cycles. We also checked the case where each cycle of the gamma burst had a different ensemble, which could precess, as in place cells²⁹. Both non-precessing and precessing bursts were efficiently transmitted (Fig 8, R,S).

Thus, in this final set of tests of our model, we predicted that mismatch detection remains sharp in spiking neurons. The mismatch response also occurs for single oddball patterns, but not gaps. We confirmed that STP on both E and I synapses are required, and showed that mismatch detection also works for 100 Hz gamma-frequency bursts riding on the theta rhythm. Together, these support possible functional relevance of these mechanisms for spatio-temporal pattern selectivity in in-vivo-like network conditions.

Discussion

We probed CA3-CA1 input-output properties for functional outcomes of pattern-specific short-term plasticity using spatially and temporally patterned optical stimuli. EI balance and summation properties tilted briefly to excitation on the second or third pulse of a burst due to differential STP profiles for excitation and inhibition. We parameterized a molecule-to-network multiscale model of CA3-CA1 network and plasticity using a series of burst, summation, and Poisson train inputs of optically defined input patterns. The model predicted that STP dynamics of E and I inputs provide a mechanism for single cells to detect transitions in input pattern sequences, which we then confirmed experimentally. We used the model to explore network configurations which could detect such transitions, and showed that spiking is strongly mismatch-tuned in continuous as well as theta-modulated gamma bursts.

EI balance is dynamically tuned through differential STP of E and I inputs

EI balance is a ubiquitous mechanism for controlling neural excitability. Numerous processes, including inhibitory plasticity, homeostatic plasticity, and developmental co-tuning have been proposed to underlie the origin and maintenance of EI balance ^30–38. EI balance is dynamic, due in part to differential effects of STP on excitatory and inhibitory inputs to the postsynaptic neuron ^{13–15,39,40}. Using field electrode stimulation of the Shaffer collaterals, these studies report a sustained shift to excitation during burst input¹³. In contrast, our sparse optical patterned stimuli results in a small window of escape from EI balance shortly after the start of a burst (Figure 3), following which both E and I undergo depression to restore balance (Figure 3, 8). Thus, short-term plasticity adds another layer of timing control through gating of E-I balance ⁴¹.

Adding further detail to the changes in EI balance, we find that short-term plasticity (STP) modulates neuronal summation over the course of a burst. Specifically, the strength of subthreshold divisive normalisation (SDN)⁶ evolves in a frequency-dependent manner, and is most linear around pulse 3 of a burst (Figure 3).

A resource to model the CA3-CA1 feedforward circuit

A major outcome of our study is the development and parameterisation of a multiscale, chemical+electrical+circuit level model of the CA3-CA1 feedforward circuit, which runs rapidly even on laptop hardware (∼20x slower than real time on an AMD 6800HS). Feature-wise, several studies combine signalling (typically postsynaptic spine and dendrite signalling) with electrical signalling^42–47. Some of these combine compartmentalized signalling with standard branched neuronal electrophysiology calculations, as we do. Others include detailed 3-dimensional reaction-diffusion kinetics on spines⁴⁸ or boutons^49,50 and these run much more slowly. Our current model is distinct in that it incorporates the network, as well as modeling a population of stochastic chemical synapses on the dendrite of a spiking CA1 pyramidal neuron.

Our network model is much reduced compared to models with exhaustive cellular and network-level detail⁵¹. Its simplicity enables extensive exploration of the network parameters and comparison with recorded activity under a series of well-controlled stimulus patterns (Figures 4-7).

Our model is modular. One can swap out different cellular geometries, provide distinct synaptic signaling models using SBML or compatible standards, and alter the detail of the network components. In the current study, we employed four different spine kinetic models, with different combinations of STP in E and I, simply by loading in different kinetic definitions. Thus the model is consistent with the ‘Interoperable’ and ‘Reusable’ aspects of FAIR principles ⁵².

STP leads to single-neuron mismatch detection in pattern sequences

Where can temporally precise gating be useful? At the single-synapse level, synaptic depression performs the operation of decorrelation, leading to more efficient coding ^41,53. Our observations can be framed as decorrelation of the spatial pattern (alternatively, activity vector) of multi-synaptic input converging onto a neuron³. A putative example of this is place field formation in novel environments. Dombeck and co-workers ^54,55, found that dendritic inhibition is reduced when a mouse traverses a new environment. Consistent with this, our simulations predict elevated spiking when novel patterns occur, in part due to presynaptic depression of inhibition (Figure 7,8).

There is a vigorous debate on the mechanisms of auditory mismatch negativity (e.g., ^56–58) of which the fresh-afferent model maps to part of our proposed mechanism (Fig 7A,B). Leaving aside the obvious differences between auditory cortex and hippocampus, we frame our model as a transient differential tilt in EI balance (Fig 3, Fig 8 H-K), rather than just afferent adaptation. We note that the time-scale of our mismatch response (∼125 ms) is similar to the 100 to 200 ms of the classic auditory mismatch negativity ⁵⁹. A similar computation of stimulus-specific adaptation in the auditory cortex has been modeled using a recurrent network with short-term depression⁶⁰. Our findings of pattern mismatch detection may also provide a mechanistic interpretation of published results showing behavioural and physiological selective responses to changes in visually delivered patterns⁵⁴. The authors found that a rate model with short-term depression was able to replicate the experiments. In a recent study, signatures of prediction error (implemented by absence of conditioned response) were observed in recurrent cortical networks in organotypic preparation⁶¹. These prediction errors were visible after a 24-hour long training period that achieved the associative conditioning. Thus the study showed the role of inhibition stabilized attractor networks of cortical circuits to produce prediction-error like responses in single neurons. Our model achieves equivalent computations with simple feedforward connectivity at a single-neuron level, without training, and with direct linkages to the underlying physiological mechanisms.

Both our experiments and models show selectivity for specific transitions. This suggests a further computation of sequence disambiguation. Consider sequence S1 AAAABBBBCCCC and sequence S2 FFFFBBBBGGGG, overlapping in ensemble B. We stipulate the presence of neuron S1_AB which detects the A->B transition, and so on. Then S1 is uniquely identified by the activity of S1_AB,S1_BC, and S2 by S2_FB,S2_BG. This is a rapid and parsimonious mechanism for distinguishing two intersecting sequences such as time-cell sequences ⁶², on sub-second time-scales.

Mismatch detection is pronounced in the gamma band

Gamma oscillations have been implicated in memory formation and communication between brain regions ^63,64. Gamma bursts are typically modulated by the theta rhythm²⁹, which is itself closely coupled to sensory intake such as sniffing^65,66 and whisking⁶⁷, and also to movement ⁶⁸. Our proposed mechanism for mismatch detection is enhanced precisely in the gamma range (Figure 8). Further, we find that theta-modulated gamma bursts perform mismatch detection between theta cycles (Figure 8P,Q). We propose that sensory snapshots obtained through sniffs or whisks provide the CA1 with theta-modulated gamma bursts of stimulus-specific activity vectors. When inputs change from one theta cycle to the next, individual cells in CA1 light up to detect specific changes between stimuli, such as odor 1 changing to odor 2, or texture 1 to texture 2 (Figure 7, 8). This implements a dimensionality reduction of a sequence of sensory inputs into a trajectory through a space of stimulus transitions.

Methods

Animals

Experiments were conducted on a mouse line with the expression of cre-recombinase specific to hippocampal CA3 (Grik4Cre) obtained from Jackson laboratory (ID: 006474). The animals were housed in a cage with up to 6 cage mates and under a 14/10 light/dark cycle. Food and water were provided ad libitum. The housing and health were managed as per the National Centre for Biological Sciences Institutional Animal Ethics Committee (IAEC) guidelines (Project: NCBS-IAE-2022/8 (R2ME). Only transgenic F1 mice, identified with genotyping, of both sexes were used for experiments.

Virus Injection

Transgenic mice between the ages 27-35 days were selected for virus injections. These animals were transferred to a separate holding facility where they were kept till the terminal experiment.

Virus

We used an adeno-associated virus AAV-LEX-rev-ChR2-tdtomato (Addgene 18917) to express channelrhodopsin2 in CA3 pyramidal cells in a cre-lox-dependent manner. The virus was diluted using sterilised PBS to a titre of 1-5 x 10¹² GC/ml and stored at −80°C in aliquots. The animals were injected with a dose of 300-500 nl of virus suspension in 10-20ms long pulses of pressurised nitrogen using a Picospritzer (Parker-Hannifin, USA). The animals were then allowed to recover under supervision and veterinary care. Terminal experiments were performed after an incubation period of 35 days or more, resulting in cre-lox-dependent expression of channelrhodopsin2 (ChR2) specific to the CA3 region of the hippocampus on the injected side. The virus AAV-FLEX-rev-ChR2-tdtomato was a gift from Scott Sternson (Addgene plasmid #18917; http://n2t.net/addgene:18917 ; RRID: Addgene_18917)⁶⁹.

Electrophysiology

Slice preparation

Animals aged between 56-120 days and at least 28 days after the virus injection were used for in vitro experiments. The animals were humanely euthanised using isoflurane (Forane) and brain was harvested into ice-cold slush of cutting solution (87mM NaCl, 2.5mM KCl, 7mM MgCl₂.6H₂O, 0.5mM CaCl₂.2H₂O, 25mM NaHCO₃, 1.25mM NaH₂PO₄, 75mM Sucrose) maintained over ice and aerated with carbogen (95% O₂ + 5% CO₂). After a rest period of 2 minutes, the hippocampus from the injected side was dissected out, and 350 µm thick sections were made using a vibrating microtome (Leica VT1200) submerged in the above-mentioned solution. The hippocampal sections were collected in a beaker filled with the recording solution (124mM NaCl, 2.7mM KCl, 1.3mM MgCl2.6H2O, 2mM CaCl2.2H2O, 26mM NaHCO3, 1.25mM NaH2PO4, 10mM D-(+)-glucose, pH 7.3-7.4 and osmolarity of 305-315 mOsm) at room temperature and allowed to recover for one hour in the dark before starting the experiments.

Recording Setup

Electrophysiological recording and optical stimulation was done on an upright DIC Olympus 53WI microscope. A slice hold-down (Warner Instruments WI 64-0246) was used to keep the slice from moving. The recording bath was perfused with the recording solution at a rate of 2 ml/min with a peristaltic pump. The recording solution was kept aerated with carbogen and the influx of the solution into the recording chamber was heated using an inline heater (Warner Instruments TC-324B) to 32-33°C.

Electrical activity was recorded using a headstage amplifier controlled by a patch clamp amplifier (Multiclamp 700B, Axon Instruments, Molecular Devices) and a data acquisition system (Digidata 1550B, Axon Instruments, Molecular Devices) using Multiclamp commander (Molecular Devices) and custom written protocols in Clampex (Molecular Devices). The electrophysiology data was acquired at 20kHz and filtered in the 0-10kHz band.

In some experiments, another recording electrode was used to record extracellular responses in the CA3 cell layer. Whole cell patch clamp recordings from CA3 and CA1 cells were made using borosilicate glass micropipettes pulled using tungsten filament puller (P1000, Sutter Instruments) and filled with an internal solution containing 130mM K-gluconate, 5mM NaCl, 10mM HEPES, 1mM Na₄-EGTA, 2mM MgCl₂, 2mM Mg-ATP, 0.5mM Na-GTP and 10mM Phosphocreatine, pH ∼7.3, and osmolarity ∼290 mOsm. For voltage-clamp recordings, K-gluconate was replaced with Cs-Gluconate. Field response was often biphasic with multiple cycles, hence, a peak-to-peak value was used.

The cells were either recorded in the current clamp or the voltage clamp. In current clamp recordings, the cells were maintained at a membrane potential of −70 mV. Cells were rejected if the membrane potential changed by 5mV over the course of the recording, holding current either was above 100 pA or changed by more than 25 pA to maintain the mentioned membrane potential. In voltage-clamp recordings, the cells were rejected if the series resistance crossed 25 MΩ or changed by 30% of the initial value. Liquid junction potential was calculated to be 16.926 mV according to the stationary Nernst–Planck equation ⁷⁰ using LJPcalc (RRID: SCR_025044) and was not compensated.

In optogenetic preparations, the additional stimulation that comes in due to the excitation of axons and dendrites can not be ruled out. To assess that, we patched a CA3 cell and recorded its optical depolarization caused by a single spot of light in a grid using a whole-cell voltage clamp. The CA3 cells have a large receptive field that arises due to the expression of ChR2 on their dendritic arbour. This receptive field is stronger closer to the soma but spans on both sides of the dendritic tree (Extended Figure 1A). A patched CA3 cell fired action potential with 100% success (Extended Figure 1B) when the cell layer was stimulated with patterns with multiple spots. A single spot, on the other hand, was only sufficient to induce action potentials when directly incident on soma (Extended Figure 1B).

Optical Stimulation

Pyramidal cells in the CA3 cell layer expressing channelrhodopsin were stimulated using a DMD-based patterned projector system (Polygon 400G, Mightex Systems, Canada) using a spot size of approximately 14µm x 7.8µm under a 40x water immersion objective. The spots (or squares) made a hexagonal grid with an inter-spot spacing equal to 32 µm. The grid occupied the middle 2/3rd of the projector frame (336µm x 187.2µm), and the microscope was translated on a stage to make the grid part overlay on the CA3 cell layer under a 40x objective. In 1-square stimulation experiments, the CA3 cell layer was illuminated one spot at a time with an interstimulus time of 3 seconds to avoid desensitisation of channelrhodopsin. In other experimental protocols, 5-spot and 15-spot patterns were used at various stimulation frequencies. There were five different non-overlapping patterns of 5-spots, and three non-overlapping patterns of 15-spots. Each presentation of a pattern was for 2-5 milliseconds. The light source had a total estimated maximum intensity of 6.7 mW, and for the grid size used, each spot had an absolute power of approximately 14.47 µW and a power density of 100 mW/mm². For all experiments, unless mentioned, the intensity of light was kept at 100%. Optical stimulation patterns were made and loaded onto Polygon using Polyscan2 (Mightex, Canada) and triggered using time series generated from custom protocols written in Clampex (Molecular Devices) and communicated to Polygon using Digidata 1440 DAC. A phototransistor (OPT101) was placed behind the dichroic to acquire the light stimulus waveform using the DAq and recorded along with the patch-clamp time series data.

Code and Model accessibility

All scripts and models for performing simulations for figures are presented in the GitHub site https://github.com/BhallaLab/STP_EI_paper_figs. Simulations were run using MOOSE 4.0 compiled using gcc version 11.4.0 and Python version 3.10. Short calculations were run on Ubuntu 22.04 on an 8-core laptop, and large runs were done on a 128-core server also running Ubuntu 22.04. All software used is open-source and hosted on github. MOOSE is available at https://moose.ncbs.res.in and https://github.com/BhallaLab/moose. Optimization of presynaptic models was performed using HOSS 1.0⁷¹ (https://github.com/BhallaLab/HOSS) and FindSim⁷² (https://github.com/BhallaLab/FindSim).

Presynaptic Plasticity Model

Voltage clamp data from the pulse-train experiments (Figure 3) was used to train a model of presynaptic short-term plasticity (Figure 4A, Extended Figure 3, Supplementary fig 2, 3) implemented in MOOSE ²³. The model included presynaptic signalling, activation of postsynaptic glutamate and GABA receptors, and electrophysiology of a ball-and-stick model of the postsynaptic neuron held in the voltage-clamp. The glutamate receptors were placed on dendritic spines implemented as separate electrical compartments, and the GABA receptors were on the dendritic shaft. Model fitting was performed in two stages. First, the raw voltage-clamp experimental data files in hdf5 format were analysed using a Python script and the peaks and valleys for each trial were extracted and stored in an intermediate hdf5 file. Second, the HOSS pipeline (Hierarchical Optimization of Systems Simulations) ²⁴ was used to fit the presynaptic signalling model to this data, uniquely for each cell and separately for excitatory and inhibitory synapses. There were ten to eighteen ‘experiments’ run on the model to perform the fitting for each cell: Two experiments for each frequency, two to four pulse frequencies (depending on how many were feasible in each patch recording), and an additional experiment for time-course of neurotransmitter release. All this was done for 5-square and 15-square stimulation, respectively. For each frequency (typically 20, 30, 40, and 50 Hz pulse rate), we performed peak-to-valley comparisons with experiments for eight pulses in a trial and paired-pulse response experiments focusing on the first couple of pulses. The simulations for the experiments were performed using MOOSE to compute the somatic voltage-clamp response to synaptic input on the dendrite of a ball-and-stick cellular model.

The synaptic input was obtained using mass-action simulations of presynaptic signalling for the model in Figure 4A. This, too, was part of the multiscale MOOSE model. Excitatory synapses were implemented as boutons connecting to a dendritic spine with a compartment for the spine shaft and spine head, respectively. The simulated neurotransmitter release triggered the opening of ligand-gated receptor channels for AMPA and NMDA receptors placed on the spine head. Inhibitory synapses were implemented as presynaptic boutons with the same chemical topology but different kinetics (Figure 5A) and coupled now to GABA receptors placed on the dendritic shaft.

This entire process was repeated for excitatory and inhibitory synapses, respectively. All calculations for the optimisation were run deterministically since, otherwise, massive averaging would be needed for each estimate of voltage-clamp current. Equations, kinetic parameters, and concentrations for the STP and no-STP models are in Supplementary Tables 1 and 2.

Postsynaptic neuron model

The postsynaptic neuron model (Figure 4B) was a ball-and-stick model with a soma compartment of 10 μm length and 10 μm diameter and a single dendritic compartment of 200 μm length and 2 μm diameter. The dendritic compartment was ornamented with 100 evenly spaced spines, each containing a spine shaft and a spine head compartment. The spine head had a model of the AMPA receptor, the NMDA receptor, and calcium dynamics resulting from influx through the NMDA receptor. However, calcium dynamics were not utilised for any further calculations. A separate presynaptic bouton with the glutamate release kinetic model was apposed to each spine head.

The dendritic compartment also had 200 evenly spaced inhibitory synapses, likewise apposed to presynaptic boutons with the GABA release kinetic model.

Under the passive conditions of our model (RM =1.0 Ω.m², RA = 1.0 Ω.m) the length constant of this model was ∼700 μm, which is much longer than the geometrical length. Hence space clamp problems are unlikely. Furthermore, the short electrotonic length means it does not matter if the Glu and GABA receptors are interspersed as opposed to GABA being more proximal.

The soma had voltage-gated Na and delayed rectifier K_DR channels in the Hodkin-Huxley formalism but with kinetics drawn from previous models. For most models, the Na and K_DR conductances were held to small values (6 and 3.5 Siemens/m^2, respectively) as we were interested in linear somatic integration, but for spiking simulations, we increased Na and K to 400 and 450 Siemens/m^2. Equations and kinetics for the voltage-gated and ligand-gated channels in the model are presented in Supplementary Table 3.

Network model

We implemented a simple thresholding matrix connectivity network model. This had a 16x16 array of CA3 neurons and another 16x16 array of inhibitory interneurons (Figure 5C). Excitatory projections from CA3 to the 100 synapses on the CA1 neuron were implemented as a 256x100 matrix. Default entries were zeros, and ones were filled in at random with a probability pCA3_CA1 as indicated in the table. Similarly, projections from CA3 to interneurons were defined as a 256x256 matrix with connection probability pCA3_Inter. Projections from the interneurons to the 200 GABA synapses on the CA1 neuron were a 256x200 matrix with default connection probability pInter_CA1. All told, there were just seven parameters in this model (Table 1).

Each optically stimulated spatial pattern of activity in CA3 was converted into a vector of ones and zeros by thresholding the CA3 membrane potential, calculated as per the next section. Activity on the vector of glutamatergic synapses was:

Glu synapse vector:

Synaptic activity was thresholded at 1 to drive the presynaptic signalling. Similarly,

Interneuron vector :

GABA synapse vector:

Stimulus model

To deliver the stimulus, we defined a set of eight input patterns. These were hand-coded as 8x8 grids of ones and zeros, set up to resemble the letters A through H. We expanded the 8x8 patterns into 16x16 patterns by replicating the pattern twice over the array (Extended Fig 5). For the first five patterns we zeroed out a specified number of entries on the 16x16 pattern using the parameter ‘zeroIndices’ to adjust sparseness. This zeroing out was done irrespective of the original state of the entry. For the default level of sparseness (75%) the first five patterns (A through E) had 18 to 26 nonzero squares. The last three patterns (F through H) were dense, with 144, 152 and 152 squares.

The product term:

sometimes led to very different numbers of active glutamate synapses, differing by as much as 50% between patterns. Hence we carried out a preliminary calculation to fill in entries in the connection matrix to bring the final synaptic activation for all the 5-square patterns to roughly the same level (within 15%). For the default 75% sparseness this came to between 35 and 38 active glutamatergic synapses (Extended Figure 10).

During the simulation runs we provided a time-series of patterns according to frozen Poisson timing series with a mean of 20 Hz (Figure 6) or a time-series of four patterns each repeated 8 times, at 8, 20 and 50 Hz (Figure 7). The time-series was generated with a resolution of 0.5 ms, and the stimulus duration (activation of Ca influx into the presynaptic boutons) was 2 ms. We used a fixed delay of 5 ms between the Glu and GABA synapse activation.

We implemented ChR2 desensitisation and membrane potential buildup over multiple pulses on the CA3 neurons using the following equations:

For each light pulse, we updated G and Vm as follows:

Where

G = conductance of ChR2

V_m = CA3 pyramidal neuron membrane potential

Δt = time since last stimulus

D₀ = basal desensitization of ChR2 = 0.01

D₁ = light pulse induced desensitization of ChR2 = 0.0004

T_r = Recovery time course of ChR2 from desensitization = 1.5 s

T_m = Membrane time course of CA3 neuron = 10 ms

T_C — Time course to charge CA3 neuron when ChR2 is open = 5 ms

T_O = ChR2 channel open time = 1 ms

Example pulse trains with these parameters are illustrated in Extended Figure 6.

Simulation parallelization and Visualization

Model fitting (Figure 4) and data generation for multi-trial stochastic runs (Figures 7 and 8) involved parallel calculations. The stochastic runs, in particular, were embarrassingly parallel in that they required multiple repeats of the same model with different random seeds. This was orchestrated through Python scripts using the multiprocessing library. The calculations were carried out on a 128-core AMD Epyc 7763 server. Data from these runs was organised into pandas data frames and dumped into hdf5 files in the same format as the experimental data. Figures 4, 6, 7 and 8 were generated from these data files using pandas, Matplotlib and Python scripts.

3-D model visualisation (Fig 4 and Supplementary movie) was performed using the built-in interface within MOOSE, to the 3-D graphics library vpython.

Additional information

Funding Information

AA and USB are at NCBS-TIFR which receives the support of the Department of Atomic Energy, Government of India, under Project Identification No. RTI 4006. The study received funding from SERB Grant CRG/2022/003135-G.

Author Contribution

Aditya Asopa: Conceptualization, Data curation, Analysis, Investigation, Visualization, Methodology, Writing original draft, Writing review and editing;

Upinder Singh Bhalla: Conceptualization, Resources, Simulation, Visualization, Analysis, Software, Supervision, Funding acquisition, Writing original draft, Project administration, Writing review and editing

Funding

Department of Science and Technology (SERB Grant CRG/2022/003135-G)

Department of Atomic Energy, Government of India (Project Identification No. RTI 4006)

Tata Institute of Fundamental Research

Additional files

Model parameters used in the study

Extended and supplementary data figures

Video 1. Visualization of mismatch detection in CA1 cell.