Impact of optogenetic pulse design on CA3 learning and replay: A neural model.

Optogenetic manipulation of hippocampal circuitry is an important tool for investigating learning in vivo. Numerous approaches to pulse design have been employed to elicit desirable circuit and behavioral outcomes. Here, we systematically test the outcome of different single-pulse waveforms in a rate-based model of hippocampal memory function at the level of mnemonic replay extension and de novo synaptic weight formation in CA3 and CA1. Lower-power waveforms with long forward or forward and backward ramps yield more natural sequence replay dynamics and induce synaptic plasticity that allows for more natural memory replay timing, in contrast to square or backward ramps. These differences between waveform shape and amplitude are preserved with the addition of noise in membrane potential, light scattering, and protein expression, improving the potential validity of predictions for in vivo work. These results inform future optogenetic experimental design choices in the field of learning and memory.

Introduction

The development of optogenetic techniques allowing millisecond timescale manipulation of neural dynamics in genetically targeted populations has fundamentally advanced the field of systems neuroscience (Boyden et al., 2005; for a review, see Fenno et al., 2011). These tools are now readily available and well tested, even under the variability inherent to in vivo conditions from sources such as protein expression levels and light scattering by tissue (Zhang et al., 2006; Zhao et al., 2008; Gradinaru et al., 2008). However, the optogenetic research community has yet to fully characterize the impact of stimulation parameter choice in order to optimize experimental circuit manipulations. Understanding the significance of pulse selection on circuit behavior will inform future optogenetic experiments and perhaps eventually therapeutic strategies (Okonogi and Sasaki, 2021).

The use of varying optogenetic pulse designs in the hippocampus (HPC) has shed new light on learning, memory, spatial navigation, and prospective planning. Continuous, long-duration square pulses in the HPC inhibit behavioral expressions of memory (Denny et al., 2014) and disrupt spatial working memory relying on HPC-prefrontal cortex transmission (Spellman et al., 2015). Closed-loop optogenetic inhibition of CA1 at different phases of the theta rhythm alters spatial working memory, supporting theories of HPC encoding and retrieval dynamics (Siegle and Wilson, 2014). However, only one recent in vivo study sought to directly assess the impact of optogenetic waveform design on HPC-mediated behavior. The 8 Hz sinusoidal stimulation of ventral HPC terminals in prefrontal cortex drives avoidance behavior in mice, but not 2, 4, or 20 Hz sines or square pulses of any frequency (Padilla-Coreano et al., 2019). These studies highlight a clear need to better characterize waveform design options, a question to which computational modeling is suited.

The CA3 subregion of the HPC has been used extensively as a model for learning and pattern recollection processes due to its sparse recurrent excitatory feedback architecture (McNaughton and Morris, 1987; Treves and Rolls, 1994; Hasselmo et al., 1995; Malerba et al., 2017; Malerba and Bazhenov, 2019). This anatomical structure is hypothesized to allow rapid encoding of distinct memories and stable recall through sharp-wave ripples (SPW-Rs). SPW-Rs are high-frequency oscillations initiated from combined activity in CA3 pyramidal neurons projecting to the CA1 subregion of the HPC that appear to contribute to recall and consolidation of recent memories as well as forward planning (Buzsaki, 1986; Girardeau et al., 2009; Ego-Stengel and Wilson, 2010; Pfeiffer and Foster, 2013; Roux et al., 2017; for a review, see Buzsaki, 2015). Ensembles of CA1 pyramidal neurons show consistent sequential reactivation during SPW-Rs that recapitulate recently experienced spatial trajectories in maze environments in a phenomenon termed “replay” (Lee and Wilson, 2002; Foster and Wilson, 2006). CA3 ensemble activity is similarly capable of ordinal spike sequences (Karlsson and Frank, 2009; Carr et al., 2012), and recent evidence suggests this may bias sequence direction of CA1 replays (Wang et al., 2020; Ishikawa and Ikegaya, 2020). Therefore, understanding the impact of optogenetic pulse design on circuit dynamics of HPC learning and memory will improve experimental design for future in vivo studies of memory.

The conversion of optical pulses to neuronal spiking in vivo is non-linear. CA3 and CA1 pyramidal cells exhibit spike frequency adaptation during prolonged input of sufficient strength (Madison and Nicoll, 1984; Lancaster and Nicoll, 1987; Scharfman 1993; Buckmaster et al., 1993; Hemond et al., 2008). The use of optogenetics to excite cells also results in adaptation in the opsin response (Boyden et al., 2005), although this form of adaptation has not been extensively studied. In voltage clamp, adaptation of the inward current induced by opsin activation has first-order dynamics with a time constant similar to the adaptation of firing rate caused by current injection. These considerations inform predictions of membrane activity during optical pulses.

We investigated the importance of optogenetic pulse parameters using rate-based computational models of ordinal sequence generation in hippocampal CA3 and CA1. To constrain the potential search space, we focused on single pulses of varying shape, hereafter termed waveform class (WFC). We elected to test square, single-ramp, and double-ramp pulses based on prior in vivo studies investigating ripple generation, ripple prolongation, and place-field manipulation, respectively (Stark et al., 2015; Fernandez-Ruiz et al., 2019; McKenzie et al., 2021). To account for spike frequency adaptation, all models presented here included an adaptation parameter that decreased membrane responsivity to prolonged input. Various contributions of noise and optogenetic response heterogeneity were also tested. Last, we investigated the influence of waveform shape on novel sequence learning in CA3 inspired by studies of artificial place-field induction (Diamantaki et al., 2016, 2018; Bittner et al., 2017). We found that different waveforms had variable impacts on circuit behavior, depending on the mnemonic feature investigated in each model. In particular, forward ramps were particularly suited to extending ongoing sequences, while square pulses disrupted these dynamics. Our results support the need for careful consideration of optogenetic pulse shape, amplitude, and duration to optimize desired circuit manipulations in vivo.

Results

Forward ramps prolong CA3 replay while preserving temporal sequence

Prior studies have employed non-square pulses to evoke or prolong mnemonic behavior at the circuit level (Padilla-Coreano et al., 2019; Fernandez-Ruiz et al., 2019). To directly investigate the impact of waveform shape on an ensemble-scale phenomenon, we modeled a network capable of generating consistent temporal sequences, or replays, similar to those seen during SPW-Rs (Foster and Wilson, 2006). Pyramidal CA3 neurons were connected with a forward bias, allowing activity of earlier nodes to propagate forward, but not in reverse (Figure 1A). The stimulation of pyramidal cell units included a calcium-mediated potassium current, IK(AHP), that resulted in activity rate adaptation, as has been noted in pyramidal cells of the HPC (Figure S1A; Madison and Nicoll, 1984; Lancaster and Nicoll, 1987; Scharfman 1993; Buckmaster et al., 1993; Hasselmo et al., 1995; Hemond et al., 2008) and other cortical structures (Barkai and Hasselmo, 1994). Adaptation to non-square waveforms matching those used in subsequent manipulations of replay was also simulated (Figures S1B and S1C). To model in vivo optogenetic prolongation of SPW-R activity, a square, 20 ms cue pulse (“cue”) to the first pyramidal node of the network was delivered to elicit a short series of replay activations. Table 1 summarizes the parameter set that was chosen such that fewer than half of the nodes reached activity threshold under the cue pulse alone (Figure 1B).

Figure 1

Forward ramps prolong CA3 replay while preserving temporal sequence

(A) Left: network connectivity diagram of the CA3 ripple extension model. For clarity, only recurrent connections from the first node are diagrammed. See STAR Methods for full description and Table 1 for parameters. Right: summary of the synaptic weight matrix, W, which employs a limited forward bias such that activity is able to propagate to nearby forward nodes. Weight strength decays linearly on the diagonal to reflect weaker involvement of more distal sequence nodes.

(B) Control replay event evoked by 20 ms square cue pulse (black) without an additional waveform. Blue through magenta lines indicate successive pyramidal node activations. Black triangles indicate threshold crossing points.

(C) Example optogenetic inputs with 50% ramp percentage delivered to the network, each starting with the same 20 ms cue pulse (yellow). (i) FR IMA, (ii) DR IMA, (iii) BR IMA, (iv) FR IP, (v) DR IP, (vi) BR IP.

(D) Example simulations using 100 ms pulses and 50% ramped waveforms to prolong ongoing ripple activity. Red through yellow lines represent successive pyramidal node activations. Black triangles indicate threshold crossing points after the onset of the secondary optogenetic pulse (gray shade). Each waveform corresponds to the example pulse delivered in (C).

(E) Parameter space characterization of ramp percentage and input duration of different WFC. Each panel corresponds to the WFC pictured in (C). Heatmaps represent degree of temporal disruption (effect size) of the interthreshold-interval (IThI) away from the IThI of the control sequence.

(F–I) Quantifications of WFC impact on ripple prolongation across ramp percentages. WFC lines as follows: red, FR; dark blue, DR; light blue, BR; solid, IMA; dashed, IP. Gray bands indicate 95% confidence intervals around mean of a 1,000-sample shuffle in (F) and (I).

(F) Mean replay sequence length (number of suprathreshold nodes) obtained from a shuffle of the simulations in (E).

(G) Minimum temporal disruption.

(H) Shortest pulse duration achieving the minimum temporal disruption displayed in (G).

(I) Mean temporal disruption.

Abbreviations: FR, forward ramp; DR, double ramp; BR, backward ramp; IMA, iso-max amplitude; IP, iso-power; P, pyramidal node; I, interneuron node, WFC, waveform class.

See also Figures S1–S3.

Table 1

CA3 ripple extension model parameter summary

Parameter	Meaning	Value
Acue	cue pulse strength	1
A	afferent pulse max amplitude	0.09 or ∗
P	pyramidal membrane potential	∗
W	recurrent weight matrix	0.036 max
W′	feedforward excitation	0.05
I	interneuron membrane potential	∗
H	feedback inhibition	0.035
H′	recurrent inhibition	0.003
η	passive decay rate	0.01
θP	pyramidal linear threshold	4
θI	interneuron linear threshold	4
θA	activation threshold	10
T	simulation duration (ms)	1,000
C	Ca2+-current strength	∗
μ	K+-current strength	0.01
γ	voltage-dependent Ca2+-current strength	0.001
ω	weight matrix ceiling	0.001
θC	calcium linear threshold	4
Ek	reversal potential of K+	−10

Related to Figures 1, 2, and S1–S6. An asterisk indicates a parameter that varied during the simulation, such as the instantaneous estimation of the neuron’s membrane potential. For the auto-recurrent weight matrix, W, the value indicated is the highest strength for auto-recurrent excitation. Forward connection to the next node was half as strong and one-quarter as strong as the second following node. All other connections were 0 (see Figure 1A). The values for C, μ, γ, ω, θC, and E were used as the adaptation parameters in all models (see Equation 3).

Figure 2

Differential effects of waveform shape on CA3 replay extension remain despite heterogeneity of optogenetic stimulation

(A) Illustration of different optogenetic heterogeneity components included in the simulation. (i) Random neuron starting positions in 3D space relative to a light source lead to suboptimal irradiance and lower gain on afferent optogenetic pulses for distal neurons. (ii) Opsin protein expression heterogeneity leads to inefficiency of translating optogenetic pulses into membrane voltage depolarization. (iii) Example of low-amplitude white noise applied to each neuron throughout the simulation.

(B) Example replays in CA3 pyramidal neurons with noise and optogenetic response heterogeneity outlined in (A). (i) Control replay cued by a 20 ms pulse. (ii) Replay extension by a 100 ms, 50% ramped FR IMA or (iii) FR IP pulse.

(C) Example parameter space searches of ramp percentage and input duration for the FR IMA WFC without noise (top row) and with all noise (bottom row). Left column: heatmap color indicates temporal disruption (effect size) relative to control replay (Bi). Maximum value capped at 5 for comparison purposes. Right column: heatmap color indicates total sequence length before and after optogenetic pulse.

(D) Same as (C) but for the FR IP WFC.

(E) Comparison of the mean temporal disruption for each WFC with all noise sources, across ramp degrees. Gray bands indicate 95% confidence interval of a 1,000-sample shuffle of input duration.

(F) Same as in (E) but for mean sequence length. Abbreviations: FR, forward ramp; DR, double ramp; BR, backward ramp; IMA, iso-max-amplitude; IP, iso-power; WFC, waveform class.

See also Figures S4–S6.

The replay sequence was perturbed 150 ms after the cue offset by a single pulse of external input delivered across all pyramidal neurons, modeling optogenetic manipulations common to closed-loop in vivo studies. Many combinations of input durations and waveform shapes yielded a full sequence retrieval of 15 nodes (Figure 1F). The specific 150 ms post-cue delay was found to elicit the longest replay sequences with minimal disruption across a number of input durations (Figure S2). Two general types of waveform shapes were tested: squares (0% ramp) or ramped (Figure 1C). Ramps were further divided between forward (FR; Figures 1Ci and 1Civ), double (DR; Figures 1Cii and 1Cv), and backward ramps (BR; Figures 1Ciii and 1Cvi). Finally, two different amplitudes were tested for each waveform, iso-max-amplitude (IMA; Figures 1Bi–1Biii) and iso-power (IP; Figures 1Biv–1Bvi), relative to the template square pulse.

The parameters of waveform duration and ramp percentage were varied systematically across each of the waveform classes of square, FR, BR, and DR and IMA and IP pulses (Figure 1C). Figure 1D (top) shows an example of a square pulse (0% ramp) capable of prolonging the replay sequence by bringing all later nodes above the threshold. Ramp waveforms were also capable of prolonging the replay sequence in full (Figures 1Di–1Dvi). When delivering “optogenetic” pulses across the network to prolong replay activity, early threshold-crossing nodes were not reactivated due to feedback inhibition from the recruited interneurons (Figure S3). A full search of the parameter space along the axes of waveform duration (0–250 ms) and ramp percentage (0%–100% ramp) was made for each waveform class (Figure 1E).

To determine the efficacy of different waveforms at prolonging replay, the average sequence length generated under each waveform shape and ramp degree was found using a shuffling procedure (Figure 1F; n = 1,000 samples). On average, all pulses prolonged the sequence by recruiting additional neurons beyond the seven active nodes recruited by the cue pulse alone (Figure 1B). Most pulses fully extended the sequence to the 15 node maximum. There was a small improvement in sequence length using IP waveforms compared with matching IMA ramps (Figure 1F [dashes versus solids]). The DR and BR (dark and light blue) outperformed the FR waveforms by producing slightly longer sequences on average.

Replay sequences are important for memory consolidation and expression (Girardeau et al., 2009; Ego-Stengel and Wilson, 2010; Roux et al., 2017). We hypothesized that reproducible ordinal relationships between neuronal activity are an important component to memory recall. Therefore, the average temporal difference between consecutive nodes reaching the activity threshold, or interthreshold interval (IThI), was used as a metric to quantify temporal reproducibility of the replay. For each simulation, a measure of the temporal disruption (effect size, Cohen’s d) was calculated between the IThI of the control and the IThI of the waveform class manipulation. Lower effect size indicates lower disruption of the natural sequence timing by the optogenetic pulse.

Waveform shape, pulse duration, and ramp degree all had an impact on the level of temporal disruption away from control (Figure 1E). Below a certain input duration, no pulse was capable of prolonging the replay activity (lowest input duration rows in Figure 1E). Darker bands indicate optogenetic pulse duration/ramp percentage combinations that yielded minimal sequence timing disruption. Figure 1G summarizes the smallest temporal disruption at each ramp degree for each waveform class. A non-significant trend for higher ramp percentage and lower minimum disruption was found for each IMA pulse shape (Figure 1G [solid lines]). Among the IP pulses, there was a small, negative correlation between ramp degree and minimum disruption under the FR (red dashed: r = −0.58, p = 0.006), a positive correlation under the DR (red dashed: r = 0.46, p < 0.038), and no correlation under the BR. These results (Figure 1G) indicate that parameters could be found across ramp degrees and input durations that supported extended sequences similar to the control activation timing, with only a minor effect of some waveform shapes on this outcome.

The pulse duration yielding this lowest temporal disruption was also calculated (Figure 1H). Among the IMA pulses, all had a positive correlation between ramp percentage and pulse length of minimum disruption (FR, r = 0.83; DR, r = 0.91; BR, r = 0.82; p < 0.001 for all), indicating that longer pulses were necessary to achieve better temporal sequence prolongation. IP ramp waveforms maintained this trend for forward and backward ramps (FR, r = 0.86, p < 0.001; BR, r = 0.62, p = 0.002) and identical lowest pulse duration for the DR condition (Figure 1H [dashed]). These results suggest that longer duration ramping pulses are required to achieve more natural sequence prolongation.

To test the average capacity for each pulse shape and ramp percentage to prolong replay behavior with minimal IThI disturbance, an effect size curve was calculated. The average of a random shuffle with replacement across input durations was taken for each ramp percentage (Figure 1I; n = 1,000 samples). IMA ramps (Figure 1I [solid]) tended toward lower average sequence disturbance (lower effect size) than their IP counterparts (Figure 1I [dashed]), especially at higher ramp percentage. FR waveforms also tended to preserve sequence dynamics better than DR waveforms. Both FRs and DRs disturbed the sequence timing far less than BR waveforms. These timing disturbance trends are consistent with the visualization of the search space in Figure 1D, which demonstrates dark, “wedge-like” regions of low perturbation in the IMA waveform class panels (i–iii), especially the FR and DR shapes. Similar patterns of timing disruption as a function of waveform shape and ramp degree were found for the interneuron populations (Figure S3).

These manipulations of replay-like sequence activation in a CA3-circuit model demonstrate the impact of optogenetic pulse shape on the temporal reliability of ensemble-scale phenomena. All waveform classes could dramatically extend the sequence beyond the control. IP waveforms of any shape tended to disrupt sequence timing except at very short input durations (Figure 1E). Highly ramped waveforms with lower amplitude have larger duration windows to achieve natural replay timing (Figure 1E).

Differential effect of waveform shape on CA3 replay extension remains despite heterogeneity of optogenetic stimulation

One of the greatest strengths of optogenetics is its utility in vivo, but real experimental conditions are characterized by heterogeneity, which leads to variability of neural responses. We therefore sought to challenge the above model by implementing several potential sources of variability or “noise,” including light scattering by tissue, opsin protein expression efficiency, and variable membrane state of the neurons (Figure 2A).

First, the position of the pyramidal units in the replay extension model was randomly determined in a 3D CA3 pyramidal layer (STAR Methods; Hussein and George, 2009; Jinno and Kosaka, 2010; Öz and Saybaşılı, 2017). The distance from a simulated laser source above the layer determined whether each node received sufficient light or whether light scattering imposed a distance penalty on the optogenetic drive to individual neurons. Without other noise sources, similar outcomes between ramp degree and mean temporal disruption and replay sequence length were found in these simulations for all waveform classes, assuming a source output of 10 mW (Figure S4D; Redondo et al., 2014; Chen et al., 2019). Lower source outputs produced similar outcomes in most waveform classes.

Modern opsins have been designed to optimize high expression efficiency without overexpression, which can lead to deleterious intracellular aggregations (Gradinaru et al., 2008; Zhao et al., 2008). Even so, some variability in optimal protein expression, dependent on factors such as the promoter, must be assumed (Zhang et al., 2006). We therefore calculated an expression efficiency for each neuron based on a random normal distribution (see STAR Methods). In the absence of other noise sources, protein expression variability yielded shorter replay sequences and higher temporal disruption as a function of protein expression variance, although the overall correlations to ramp degree were preserved for all waveform classes (Figure S5D).

Another source of variability, fluctuations in the membrane voltage, was modeled by the addition of white noise inputs of varying amplitude. This noise could depolarize or hyperpolarize the membrane voltage of a node prior to the replay cue and throughout the replay, allowing examination of replay extension stability under noisy membrane conditions. Without other variability sources, replay extension by various optogenetic waveform classes was affected by membrane voltage noise as a function of the white noise amplitude (Figure S6). At low noise amplitude, minor increases in the temporal perturbation were observed for all waveform classes, while at very high amplitude, the effect of waveform class and ramp degree on sequence reliability was abolished (Figure S6D). Replay sequence length was largely unchanged by noise amplitude.

Combining these sources of response variability, we tested the same optogenetic waveform classes on replay extension in the CA3 model under more experimentally relevant conditions. Despite the addition of light attenuation through tissue, variability in protein expression efficiency, and membrane activity fluctuations, the control replays and outcomes of replay extension appeared minimally affected under most simulations (Figure 2B). A full characterization of the search space of ramp degree and input pulse duration was conducted for each waveform class (examples of FR IMA and IP waveforms; Figures 2C and 2D). The addition of noise yielded a correlation between ramp degree and minimum temporal disruption in the FR IMA (r = −0.48, p = 0.027) and DR IP (r = 0.44, p = 0.045) shapes, but no other significant correlations. However, the mean of the shuffled temporal disruptions revealed similar dissociations between waveform classes across ramp degrees as in the noiseless model (Figure 2E versus Figure 1I). Lower-amplitude IMA waveforms caused less temporal perturbation (Figure 2E), but shorter sequences (Figure 2F), compared with IP waveforms of the same shape, as a function of ramp degree. Under noisy conditions, DR and BR IMA waveforms performed more similarly. As in the noiseless model, FR IMA waveforms yielded minimum temporal disruption with modest replay extension length, especially at longer ramp degrees. Therefore, the model predictions for different waveform classes remain broadly similar under the tested sources of variability.

Low-amplitude forward ramps optimally prolong replay in CA1

CA3 activity during sharp waves contains sequential information (Karlsson and Frank, 2009; Carr et al., 2012), but the preponderance of studies investigating SPW-Rs have been carried out in the CA1 region. To confirm the impact of waveform shape on ripple prolongation, we repeated the above waveform class, ramp degree, and input duration searches in a model of region CA1 architecture (Figure 3A; see Table 2 for parameters). Note that CA1 pyramidal cells had weak intra-areal connectivity to one another but no auto-recursion (Deuchars and Thomson, 1996; Shepherd, 2004; Malerba and Bazhenov, 2019). A 15 node long replay was cued by a 20 ms pulse in CA3 and propagated to CA1 by synaptic connections (Figure S7). Parameters were selected such that a partial, eight-node replay was read out in the CA1 under the cue pulse alone (Figure 3B). After a 150 ms delay from cue offset, a single waveform pulse was delivered to CA1 pyramidal cells to prolong the sequence by enhancing the response to sequential CA3 inputs (Figure 3C).

Figure 3

Low-amplitude forward ramps optimally prolong replay in CA1

(A) Network connectivity diagram of the CA1 replay extension model. For clarity, only recurrent connections from the central node are diagrammed. See STAR Methods for full description and Table 2 for parameters.

(B) Control replay event across regions CA3 and CA1 evoked by cue pulse. CA3 pyramidal nodes (hot colors) generate a full, suprathreshold replay. CA1 pyramidals (cool colors) and interneurons (not shown) receive input from pre-synaptic CA3 pyramidals, generating a decaying replay sequence of eight nodes (threshold crossing times marked by black triangles).

(C) Example CA1 replay extension by various classes of optogenetic pulse delivered 150 ms after CA3 replay onset. Only suprathreshold nodes after pulse onset are marked. Black line below each plot illustrates ramp input type and timing, not to scale.

(D) Heatmaps of parameter space searches of ramp degree and input pulse duration for IMA waveforms. Top: heatmap color indicates temporal perturbation. Bottom: heatmap color indicates total replay sequence length of suprathreshold nodes.

(E) As in (D) but for IP waveform classes.

(F) Summary of mean temporal disruption across WFCs for CA1 replay extension model.

(G) Same as (F) but for mean sequence length (all suprathreshold nodes). Gray bands indicate 95% confidence interval around a 1,000 sample shuffle of input duration.

Abbreviations: OL, overlap; FR, forward ramp; DR, double ramp; BR, backward ramp; IMA, iso-max-amplitude; IP, iso-power; P, pyramidal node; I, interneuron node; WFC, waveform class.

See also Figure S7.

Table 2

CA1 ripple extension model parameter summary for both CA3 and CA1 regions

Parameter	Meaning	Value
Acue	cue pulse strength	1
A	afferent pulse max amplitude	0.1 or ∗
Pi	CA1 Pyr membrane potential	∗
PCA3r	CA3 Pyr membrane potential	∗
W	CA3 Pyr to Pyr weight	0.0331 max
W′	CA3 Pyr to CA3 IN weight	0.05
WZ	CA3 Pyr to CA1 Pyr weight	0.02
WQ	CA3 Pyr to CA1 IN weight	0.02
ZQ	CA1 Pyr to CA1 IN weight	0.05
ZZ	CA1 Pyr to CA1 Pyr weight	0.002
QZ	CA1 IN to CA1 Pyr weight	0.045
Ik	CA1 IN membrane potential	∗
H	CA3 IN to CA3 Pyr weight	0.034
H′	IN recurrent inhibition	0.003
η	passive decay rate	0.01
θP	pyramidal linear threshold	4
θI	interneuron linear threshold	4
θ	activation threshold	10
T	simulation duration (ms)	1,000

Related to Figures 3 and S7. An asterisk indicates a parameter that varied during the simulation, such as the instantaneous estimation of the neuron’s membrane potential. Abbreviations: Pyr, pyramidal; IN, interneuron.

A parameter search on input duration and ramp degree was performed for each waveform class (Figures 3C and 3D). Negative correlations between ramp degree and minimum temporal disruption in CA1, indicating better timing preservation by longer ramps, were found for all IMA waveforms (r for FR, −0.80, p < 0.001; DR, −0.47, p = 0.031; BR, −0.49, p = 0.024), but for only the FR IP waveform (r = −0.74, p < 0.001) A 1,000 sample shuffle of input durations revealed the impact of waveform class and ramp degree on CA1 replay extension. IP waveforms of all shapes caused more temporal disruption to the optogenetically prolonged sequence and recruited more nodes than their IMA counterparts (Figures 3F and 3G), as found in the CA3 model. Each of the IMA waveforms had a different degree at which the average temporal disruption was minimized and sequence length maximized (Figures 3F and 3G). The FR IMA waveform optimal ramp degree occurred at 45% ramp, while the DR IMA waveform performed best at 90%–100% ramp. The BR IMA waveform disrupted the sequence least at 100% ramp, but minimally extended the sequence beyond control at any ramp degree. These findings in CA1 are consistent with the observations in the CA3 model that lower-amplitude FR and DR waveforms best balance the trade-off of sequence reproducibility and sequence prolongation length.

Low-amplitude waveforms induce novel sequence learning capable of stable replay

Finally, inspired by recent investigations of place-field induction (Bittner et al., 2017; Diamantaki et al., 2018; McKenzie et al., 2021) and the growing capacity of all-optical approaches to observe and manipulate circuit function (Rickgauer et al., 2014; Packer et al., 2015; Mardinly et al., 2018; Marshel et al., 2019), we investigated the outcome of waveform class on inducing novel sequence learning. A CA3 network was used as before, but this time with no pre-formed weight matrix (Figure 4A; see Table 3 for parameters). During the learn phase of the following simulations a 15 element sequential optogenetic pattern was presented once to the pyramidal cells of the network (Figure 4A [right, yellow pulses]). Activity induced by the afferent optogenetic input pattern, A, drove auto-associative learning between nearby pyramidal cells, but not more distal ones (Figure 4A, W post-learn).

Figure 4

Low-amplitude waveforms induce novel sequence learning capable of stable replay

(A) Left: simplified connectivity diagram of the CA3 learning model. For clarity, only the auto-recurrent connections of the central node are pictured. See STAR Methods for full details and Table 3 for parameters. Right: model phases and example optogenetic input patterns leading to weight change. Learn phase: the uniform zero auto-recurrent matrix, W, is presented with a single pattern of temporally sequenced inputs and learns new synaptic weights. Test phase: pattern replay is probed with a 20 ms square cue pulse.

(B) Left: control replay event evoked by cue pulse (black) using (right) a pre-formed weight matrix. Blue through magenta lines indicate successive pyramidal node activations. Black triangles indicate threshold crossing points.

(C) Examples of learn phase optogenetic input patterns, A. Each pattern element had a total 80 ms duration, with variable shape, amplitude, and pulse OL (60% OL pictured) with neighboring elements (red hatched box).

(D) Example test phase activity evoked by cue pulse. Red through yellow lines represent successive pyramidal activity resulting from the pattern encoded during the learn phase (60% pulse OL), summarized in the W weight matrix (inset).

(E and F) Parameter space characterization across learn phase OL and ramp percentage. (E) Heatmaps of replay sequence length (number of suprathreshold nodes) for each WFC in (C). (F) Heatmap of temporal disruption of the interthreshold interval (IThI) away from the control sequence IThI (effect size) for each WFC in (C).

(G–J) Quantification of the WFC impact on pattern learning across ramp percentages. WFC lines are as follows: red, FR; dark blue, DR; light blue, BR; solid, IMA; dashed, IP. Gray bands indicate 95% confidence intervals around mean of a 500 sample shuffle in (G), (H), and (J).

(G) Mean replay sequence length (number of suprathreshold nodes).

(H) Mean IThI (ms).

(I) Minimum temporal disruption.

(J) Mean temporal disruption.

Abbreviations: OL, overlap; FR, forward ramp; DR, double ramp; BR, backward ramp; IMA, iso-max-amplitude; IP, iso-power; P, pyramidal node; I, interneuron node; WFC, waveform class. See also Figure S8.

Table 3

CA3 learning model parameter summary

Parameter	Meaning	Value
Acue	cue pulse strength	1
Ai	learn pulse max amplitude	0.5 or ∗
P	Pyr membrane potential	∗
B	weight matrix ceiling	0.035
Wij	auto-recurrent weight matrix	0.035 max
W′	feedforward excitation	0.05
I	IN membrane potential	∗
H	feedback inhibition	0.05
H′	recurrent inhibition	0.003
Η	passive decay rate	0.01
ε	learning rate	0.001
ψLearn	acetylcholine level learn	0.1
ψTest	acetylcholine level test	1.0
θP	Pyr linear threshold	4
θI	IN linear threshold	4
θA	activation threshold	10
TLearn	learn phase duration	1,500
TTest	test phase duration	1,500

Related to Figures 4 and S8. An asterisk indicates a parameter that varied during the simulation, such as the instantaneous estimation of the neuron’s membrane potential. The value of each connection in the auto-recurrent weight matrix, W, was allowed to potentiate up to a ceiling set by the parameter B during the learn phase (Figure 4A). Abbreviations: Pyr, pyramidal; IN, interneuron.

To study the impact of waveform shape on pattern learning, square, FR, BR, and DR and IMA and IP pulses were used as the input pattern elements (Figure 4C). The timing overlap (OL) between one pattern element and the next was tested between 0% and 100% OL. Another set of simulations examined the length of these pulses and found that 80 ms pulses yielded new sequence learning over a range of OLs (Figure S8). After one-trial learning, the capacity of the network to recall the pattern was probed in the test phase using a square, 20 ms cue pulse delivered to the first node of the network. Example test phase sequences using 80 ms pulses of 60% OL are presented in Figure 4D for each waveform class, including square (top), along with the learned weight matrix that produced them (inset). Very low pattern OL yielded minimal learning and singular node activations during the test phase across all pulse shapes. Conversely, very high pattern OL caused oversaturation of synaptic weights, leading to excessive self-excitation during the test phase. We therefore restricted our analysis to the band of OL percentage in which a dynamic range of learning occurred.

A number of metrics were used to analyze potential differences in sequence learning across ramp percentage and waveform shape. First, the average number of nodes above threshold, or sequence length, for each ramp percentage was calculated. Sequences of all lengths, between 1 and 15, appeared in the test phase for all shapes, across all ramp percentages (Figure 4E). The average sequence length generated in the test phase by a range of OLs was found using a shuffling procedure (Figure 4G; n = 500 samples). IP ramps tended to produce longer sequences compared with matching IMA waveforms. IMA pulses managed to perform as well as the square pulse at all ramp percentages in terms of sequence length. We also observed that IP waveforms of any shape tended toward shorter mean IThI than IMA waveforms (Figure 4H; n = 500 shuffle). These findings of long sequence but short IThI interval raised the possibility that those recalled sequences were temporally compressed, as was observed in the ripple extension model.

We therefore assessed test phase replay activity against a control sequence to gauge temporal sequence reproducibility (Figure 4F). The control activations were generated using a network with the same parameters as in the test phase, but using a pre-formed weight matrix. The synaptic strength of this matrix was chosen to enable a seven-node-long sequence activation after cue pulse (Figure 4B). Our weight selection allowed for combinations of ramp degree and OL that yielded high similarity (low temporal disruption) to control replay timing across all waveform shapes (Figure 4I). This choice was further validated by in vivo literature indicating that the majority of SPW-Rs include a finite number of place fields that tile less than the total potential length of the learned spatial trajectory (e.g., the stem of a T maze; Fernandez-Ruiz et al., 2019).

The IThI between serial node activations was again used to quantify the temporal reproducibility of the test phase replay sequence. First, the degree of IThI disruption between each waveform manipulation and the control replay was calculated (Cohen’s d effect size). The OL degrees achieving the minimum sequence perturbation for each ramp percentage and waveform class are summarized in Figure 4I. We observed a significant negative correlation (Pearson’s r) between ramp percentage and smallest sequence timing disruption in the DR IP waveform (r = 0.57, p = 0.006) only. Finally, the average sequence timing disruption for each ramp percentage was calculated using a shuffle of OL percentage (Figure 4J), again restricted to OL percentages where dynamic response to learning input occurred. IP waveforms tended to produce a higher average perturbation regardless of OL percentage compared with IMA waveforms. This finding complements our previous results in the ripple extension model, once again emphasizing more disruption to standard sequence timing with higher amplitude waveform choices (Figures 4I and 4J).

Taken together, these findings support the hypothesis that pulsed waveforms with longer ramps and lower power can induce synaptic plasticity and allow for naturalistic recall after single-trial learning. However, the range of pattern element OL may be quite restricted in order to produce recognizable sequences during later recall.

Discussion

This study employed a computational approach to investigate the importance of optogenetic waveform parameter selection, or waveform class, on neural response in a hippocampal model of memory. We employed a non-spiking, continuous rate model of neural dynamics (Wilson and Cowan, 1972, 1973) to show proof of principle that waveform class can have an impact on ensemble activity within a constrained search space of waveform shape, amplitude, and duration. We further narrowed the question by focusing on the use of optogenetic pulses to elicit mnemonic sequence activity in the hippocampal CA3 and CA1 subregions. We tested our waveform design options in a model of replay sequence extension with and without noise and in a model of novel sequence learning. As a main result of our simulations, we find that sequence extension models robustly showed more temporally reliable sequence replay when using forward or double ramp waveforms, in contrast to square or backward ramp waveforms.

Hippocampal ripple extension models

In vivo replays are associated with SPW-Rs (Lee and Wilson, 2002; Foster and Wilson, 2006), a characteristic 140–220 Hz oscillatory activity in CA1 with distinct signatures compared with other oscillations (Buzsaki et al., 1992; Ylinen et al., 1995; Sullivan et al., 2011). Our replay extension models employed pre-formed, asymmetric weight matrices to produce forward-biased sequential activity, or replay, in CA3 alone or in a dual-region CA3-CA1 simulation (Figures 1A, 1B, and 3B). This synaptic weight structure is analogous to well-trained animals capable of recalling prior learning to plan future forward trajectories through an environment (Diba and Buzsaki, 2007; Carr et al., 2011; Pfieffer and Foster 2013; for a different weight architecture, see the discussion of the learning model). Recent research highlights a gradient of replay directions, which favors forward replays in CA1 as animals become familiar with a route (Shin et al., 2019; Igata et al., 2021). Delivering a cue pulse alone yielded only a partial sequence, with later ranked neurons falling below threshold. Delivery of a single optogenetic pulse across the whole of CA3 (Figure 1D) or CA1 (Figure 3C) prolonged the replay sequence by raising late-ranked nodes above threshold, consistent with in vivo manipulations (Fernandez-Ruiz et al., 2019). The optimal timing of the pulse relative to the event onset was found at 150 ms (Figure S2), but this should be interpreted only insofar as there likely is some optimal window of stimulation delivery, given successful ripple prolongation in vivo with closed-loop detection and stimulation delays of tens of milliseconds. The recruitment of interneurons prevented early sequence cells from reactivating during long optogenetic pulses (Figure S2).

Spike frequency adaptation refers to a gradual slowing of action potential rate during firing to prolonged input, usually due to activation of potassium currents, which may act as an important influence on cortical processing efficiency (Wainwright 1999; Gutierrez and Denève, 2019). Many hippocampal pyramidal units in both CA3 and CA1 display this property of spike frequency adaptation during sustained depolarization (Madison and Nicoll, 1984; Lancaster and Nicoll, 1987; Scharfman 1993; Buckmaster et al., 1993; Hemond et al., 2008), as do neurons activated with optogenetics (Boyden et al., 2005). While the spike adaptation rate is heterogeneous (Barkai and Hasselmo, 1994), in HPC it generally strengthens with a time course on the order of about 50–100 ms (Madison and Nicoll, 1984; Buckmaster et al., 1993 Verma-Ahuja et al., 1995). The time course of this neural adaptation (Barkai and Hasselmo, 1994) is similar to the time course of the adaptation of the current induced by opsin activation, as shown in voltage clamp recordings (Boyden et al., 2005), and is mirrored by the time constant achieved in our isolated CA3 pyramidal cell model (Figure S1). This time course of adaptation is far longer than the participation of a single neuron in an SPW-R, and may therefore have less impact on a transient dynamic such as a replay. Under our circuit model conditions, we observed that even strongly activated pyramidal units spent less than 50 ms above the high-activity threshold regardless of pulse duration or shape, due to the adaptation current and strength of feedback inhibition from interneurons (Figures 1 and S3). Due to a lack of in vivo data characterizing the nature of adaptation to optogenetic stimulation separate from other forms of adaptation, our model lumps opsin and spike frequency adaptation into a single, conservative time constant.

Waveform class on ripple prolongation

In agreement with our central hypothesis that waveform shape affects circuit behavior, some shapes prolonged ongoing ripples better than others. In both CA3 and CA1, square pulses robustly recruited more distal sequence elements (Figures 1D and 3C [top]). However, the sharp onset of square optogenetic drive disrupted natural sequence timing severely (Figures 1E, 3D, and 3E [leftmost column of each]), resembling studies using square pulses to inhibit ongoing ripples (Fernandez-Ruiz et al., 2019). Only short-duration square pulses, around 20 ms duration, preserved natural sequence dynamics, more analogous to studies using short optogenetic square pulses to identify tagged cells without altering excitability or tuning properties (Tanaka et al., 2018). By contrast, pulses initiated with a forward ramp (FR, DR) better preserved sequence timing relative to square-onset pulses (Figures 1D and 3C; examples of 50% ramp simulations). In CA1, the temporal disruption and extended sequence length were generally lower compared with CA3 alone. In addition, the modeled CA3-CA1 architecture resulted in a different pattern of waveform optimization compared with extending replay in CA3 alone (Figures 1F and 1I versus Figures 3F and 3G). For example, in CA3, the FR IMA waveform at 100% ramp yielded the least temporal disruption but lowest sequence length, while in CA1, a 45% FR IMA shape yielded the least disruption and longest replay extension. These results must be interpreted cautiously within the limitations of the model, but the predictions are consistent with currently available in vivo data on SPW-R manipulation using moderate (20%) ramp pulses to prolong SPW-Rs in CA1 (Fernandez-Ruiz et al., 2019).

These findings suggest a higher tolerance of the circuit to respond to gentler optogenetic input (Figures 1F and 3F). Differences in optimal ramp degree and overall sequence timing perturbation emerged between the FR and the DR waveforms, possibly due to lower total power delivered. We therefore hypothesize that under conditions of ongoing circuit activity, gentle excitatory ramps preserve existing dynamics better. This predicted effect of gentle ramp-up is further validated by the finding that IP waveforms (with steeper ramps) reduced the benefit of pulse shape on sequence timing preservation (Figures 1F and 3F, dashed lines). In CA3, strong IP pulses showed a narrow horizontal band of acceptable input durations agnostic to waveform class or ramp degree (Figures 1Diii and 1Div). We link this finding to the overrecruitment of the simulated interneuron population under high-power, steep-onset waveforms, which caused a wave of inhibition to blur or disrupt the timing of the ongoing replay. Therefore, while past experiments have used square pulses to evoke ordinal sequences in CA1 (Stark et al., 2015), prolonging ongoing network events with behavioral relevance may benefit more from lower-amplitude, shallow ramped waveforms (Fernandez-Ruiz et al., 2019). Future in vivo work explicitly testing SPW-R initiation versus prolongation under different pulse shapes will be needed to validate these predictions.

Optogenetic pulse response variability

Many factors determine the outcome of optogenetic circuit manipulations in vivo. We further tested the CA3 replay extension model by including opsin expression heterogeneity, irradiance reduction due to light scattering through tissue, and membrane voltage fluctuations from other sources. Testing each source of heterogeneity alone (Figures S4–S6) indicated a tolerance of the model for each type of noise over some levels of variability and a breakdown of replay reliability at more extreme levels. Irradiance drop-off through tissue can be accounted for either by increasing source power (which may lead to unwanted heating artifact; Arias-Gil et al., 2016; Shin et al., 2016) or by using multiple sources to cover larger areas (Wu et al., 2015; McKenzie et al., 2021). Much work has been done to develop more efficient opsins that traffic preferentially to the membrane and avoid under- or overexpression in vivo (for examples, see Gradinaru et al., 2008; Zhao et al., 2008). Finally, only very high white noise amplitudes relative to the optogenetic drive strength of the model were sufficient to disrupt the replay extension outcomes (Figure S6). We therefore multiplicatively combined the lower range of variability of each type (Figure 2) in a model that bore out overall findings similar to those of the noiseless model (Figure 1). Thus, the predictions of waveform class impact on replay extension are more interpretable for future experimental design. Further, this model provides a framework for testing new predictions of memory readout under specific experimental constraints and noise parameters. Modeling more specific constraints prior to in vivo optogenetic work could improve experimental yield and maximize efficient animal use.

CA3 learning model

Innovations in simultaneous multi-photon imaging and optogenetics have enabled real-time identification and single-cell-specific stimulation of neural activity (Rickgauer et al., 2014; Packer et al., 2015; Mardinly et al., 2018; Marshel et al., 2019). These so-called “all-optical” approaches have been used to induce new sequences in visual cortex (Carrillo-Reid et al., 2016) and manipulate CA1 place fields in vivo to bias behavior (Robinson et al., 2020), but novel place-field sequence induction has yet to be thoroughly investigated in the HPC. We therefore sought to model whether waveform class would have an impact on the artificial creation of new CA3 sequences to predict the results of future in vivo circuit manipulations in the HPC using all-optical techniques. We delivered sequential pulses across the model network to induce synaptic weight formation. Under an auto-associative learning rule, we found that a wide variety of waveform shapes could be used to induce sparse associations between nodes stimulated close in time (Figure 4C). During a test phase, we examined the sequence recall produced by one-trial learning. The timing OL between pulses largely determined whether sufficient synaptic weights would be formed to allow sequential replay during the test phase (see Figures 4E and 4F).

The mechanisms behind novel place-field formation have only recently come under investigation. Juxtacellular stimulation in mice during freely moving behavior using single square pulses in the dentate gyrus (Diamantaki et al., 2016), or trains of square pulses in CA3 and CA1 (Diamantaki et al., 2018), can induce novel spatial tuning in previously silent neurons. Our results are in line with these findings, as sequential stimulation of CA3 neurons induced novel synaptic weight change sufficient to support subsequent sequence replay (see Figure 4D for examples). Lower-amplitude IMA pulses of any waveform and ramp percentage induced shorter replays (Figure 4G [solid]) with lower temporal disruption of sequence timing (Figure 4J [solid]). Higher-amplitude IP pulses of any waveform allowed for longer sequences (Figure 4G [dashed]) but caused more perturbation, especially at longer ramps (Figure 4I [dashed]). This supports the conclusion that modest-power waveforms of any shape could be used to induce learning between previously unconnected nodes. However, care should be taken in interpreting these findings for a number of reasons, including the imprecision of spike times to prolonged stimuli (Mainen and Sejnowski, 1995), endogenous bursting in CA3 (Hablitz and Johnston, 1981), irrespective of stimulus-tuning (Frerking et al., 2005), and the sparsity of CA3 auto-recursion (Miles et al., 2014). This is likely reflected in the oversaturation effect at high OL percentages when too many nodes are simultaneously recruited at high activity by optogenetic drive. Finally, the behavioral relevance for randomly stimulated subsets of cells, as was investigated previously in V1 (Carrillo-Reid et al., 2016), may be harder to explore further from the sensory input interface.

In CA1, induction of new place cells during behavior from silent cells is accompanied by complex, bursting spike trains similar to previously reported plateau potentials that reflect integration of information in CA3 and the entorhinal cortex (Bittner et al., 2015; Diamantaki et al., 2018). Induction of plateau potentials by 300 ms square current injections is reported to open a multi-second plasticity window for novel place-field formation that satisfies the requirement of integrating sensory information at a behavioral timescale (Frerking et al., 2005; Bittner et al., 2017). These findings support a sensory-driven model of place-field formation onto the tabula rasa or “blank slate” of the HPC. Our CA3 learning model, which allowed for synaptic enhancement between any co-active nodes, can be interpreted within this framework. However, the propensity of any CA3 cell to initiate complex spike bursts is heterogeneous (Raus Balind et al., 2019), and the burst-potentiated inputs may be outside the typical window of Hebbian plasticity (Bittner et al., 2017).

Intrinsic connectivity and pre-existing learning are important factors in natural place-field formation and sequence encoding. In the absence of external sensory cues, hippocampal networks in CA1 naturally form recurring sequences with consistent order that “chunk” voluntary run distance (Villette et al., 2015). During periods of rest, small assemblies of CA1 pyramidal cells fire in consistent patterns that predict the rank order of place-field tuning on later exposure to novel mazes in a phenomenon termed “preplay” (Dragoi and Tonegawa, 2011, 2013). Differences in neuronal firing, excitability, and synaptic plasticity in CA1 pyramidal cells contribute to a heterogeneous population that supports long-term sequence stability and allows flexible incorporation of new information (Grosmark and Buzsaki et al., 2016). Recent optogenetic studies have found that local stimulation induced heterogeneous remapping of spatial tuning to maze locations with weak, pre-existing drive (McKenzie et al., 2021) or transient unmasking of subthreshold place-field activity (Valero et al., 2022). In contrast to the direct-stimulation studies above, these findings support a hypothesis that new learning is constrained and mapped onto “pre-wired” internal sequences that flexibly incorporate new experience. While our learning model assumed a uniform zero weight matrix with unconstrained potential for novel associations, this architecture is unlikely in vivo. Future modeling work should investigate the interaction between new sequence induction and existing network architecture.

Waveform power

The greatest commonality to waveform manipulation across the models of the HPC circuit was the impact of waveform power. In the ripple extension models of CA3 and CA3-CA1, IMA forward and double ramps proved more effective at prolonging the ripple relative to square or backward ramp stimuli. IP pulses induced closest similarity to the control ripple at lower input durations than IMA pulses (Figure 1H), but on average disrupted the sequence timing more than IMA pulses (Figure 1I). In the learning model, IP waveforms delivered during learning allowed longer replay-like sequences during the test phase (Figure 4G), but simultaneously caused more disruption in node activation timing compared with the same IMA pulses (Figure 4J). In vivo optogenetic design must account for the limitations of equipment and optical setups to achieve higher laser power delivery. More importantly, long-duration, high-amplitude pulses can result in tissue heating (Arias-Gil et al., 2016; Shin et al., 2016), even in the absence of opsins (Owen et al., 2019). Therefore, to optimize circuit manipulation with minimal tissue heating, our results predict the use of short square pulses or forward ramping, lower-power waveforms with longer duration.

Prospects for modeling waveform class optimization

Much future work will be required to guide optogenetic waveform design optimization for in vivo experimentation. More biophysically plausible models could examine the impact of pulse shape on opsins expressed preferentially in dendritic or somatic compartments. Future models of sequence learning could also incorporate a range of pre-existing network connectivity (from uniform zero as in our model to highly heterogeneous connectivity prior to new learning) not only as a means of predicting pulse impact on learning in vivo, but also to address fundamental questions surrounding the tabula rasa versus constrained network debate.

Many of the in vivo studies employing optogenetics to elicit or inhibit mnemonic behavior rely on pulse trains. For example, 20 Hz square pulses have been used to activate sparse dentate gyrus ensembles to elicit freezing (Liu et al., 2012; Ramirez et al., 2013) or induce place avoidance and place preference, depending on the activated ensemble (Redondo et al., 2014; Chen et al., 2019). Recent work has also investigated pulse trains of varying frequency and waveform shape, demonstrating preferential recruitment of hippocampal outputs on danger avoidance under 8 Hz sinusoidal, but not other frequencies or using square pulses (Padilla-Coreano et al., 2019). Sinusoidal stimulation of the basolateral amygdala differentially drove freezing (4 Hz) or non-freezing (8 Hz) in mice depending on the context of stimulation (Ozawa et al., 2020). These results strongly indicate that both the frequency and the shape of pulsatile stimulation have an impact on behavioral expression. Future modeling should therefore address the search spaces of frequency and waveform shape on different circuit architectures. The work presented here offers a useful preliminary exploration of the search space to demonstrate proof of principle for further investigation.

Conclusion

We explored how to optimize optogenetic waveform design for desirable mnemonic circuit activity. In our rate-based models of memory circuit behavior we provide evidence that forward and double ramps of moderate power are much more effective for prolonging existing network events such as SPW-Rs and replays in both CA3 and CA1. These findings were robust to the inclusion of biological detail to the models, including spike-frequency adaptation to excitation and several sources of variability inherent to in vivo optogenetic implementations. Developing a theoretical basis for optogenetic waveform design choices will advance the technical capabilities of this important tool set for research, and perhaps for clinical applications.

Limitations of the study

This work focuses solely on investigating the outcomes of single-pulse waveforms in the learning and memory system of the HPC. However, the approach could easily be adapted for use in other circuit architectures or for studying pulsatile stimulation outcomes. The predictions are drawn from a rate-based, rather than spiking, model of neuronal activity. Implementing this approach with spiking models would require redefining some of the response variables but would provide similar, systematic insight into spike timings evoked by optical manipulation.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information on methodology and code should be requested and will be fulfilled by the lead contact, Lucius Wilmerding (lwilmerd@bu.edu).

Materials availability

This study did not generate new unique reagents.

Method details

Model justification

The CA3 subregion of HPC was chosen as a model system for examining circuit behavior under differential input waveform conditions. The extensive recurrent collateral architecture has been exploited in past models of learning and recall (Wilson and Cowan, 1973; McNaughton and Morris, 1987; Hasselmo et al., 1995; Hasselmo and Wyble, 1997; Elfman et al., 2014) and more recently for its ability to spontaneously generate sharp wave population activity (Malerba et al., 2017; Malerba & Bazhenov, 2019). It is further known to produce sequential patterns (Karlsson and Frank, 2009; Carr et al., 2012) and coherent burst activity (Csicsvari et al., 2000) to cue downstream CA1 sequential neuronal reactivations in vivo (Wang et al., 2020). Therefore, the pattern recall, temporal sequence reliability and synaptic weight modification properties of the CA3 region were used to investigate mnemonic circuit outputs under different waveform shapes.

CA3 ripple extension model

We modeled self-sustained, excitatory activity in CA3 neurons based on previous work (Hasselmo et al., 1995) using a continuous firing rate model (Wilson and Cowan, 1973). To examine the impact of waveform shape on a model of forward sequenced memory activation (Diba and Buzsaki, 2007; Fernandez-Ruiz et al., 2019), we created a 30 unit network with equal excitatory and inhibitory units (Figure 1A). The auto-recurrent weight matrix, W, was initialized with a forward bias such that neurons had nonzero connectivity to themselves and linearly decreasing synaptic strength to the next two nodes, and zero connectivity in the reverse direction (Figure 1A). Connections between inhibitory and excitatory nodes were maintained on a one-to-one basis. The membrane potential equations for the pyramidals (P) and interneurons (I) were as follows:Where A represents the excitatory afferent optogenetic input to the pyramidal cell, η is the passive membrane decay toward rest, W is the weight of auto-recurrent excitation, W′ is the weight of P to I, H is the weight of I to P, H′ is the weight of auto-recurrent inhibition (I to I). Each pyramidal is indexed by i, each interneuron is indexed by k, and the auto-recurrent input to each neuron is calculated by the sum of the outputs of the rectified linear thresholding function, θ, multiplied by the weight matrix. Thus, a neuron received input only from those other neurons with activity above threshold θ and with nonzero weight connectivity. The parameters used in the isolated CA3 model are presented in Table 1.

Cells representing more distal components of the replay sequence have lower participation in SPW-Rs (Fernandez-Ruiz et al., 2019). To reflect this, the strength of the recurrent weight matrix, W, decreased linearly along the diagonal (Figure 1A, right, shift from light yellow to darker). Thus, earlier nodes benefitted from a slightly stronger connectivity than later nodes.

Given the phenomenon of spike frequency adaptation to prolonged depolarizing stimulation in CA3 (Buckmaster et al., 1993; Scharfman 1993; Verma-Ahuja et al., 1995; Hemond et al., 2008) and CA1 (Madison and Nicoll 1984; Lancaster and Nicoll, 1987), we included an intracellular-calcium-concentration-dependent potassium current in the pyramidal neuron. The update to the calcium concentration of each cell, C, was calculated as follows:Where E is the reversal potential of potassium, μ is the strength of the Ca2+-dependent K+-current, γ is the strength of the voltage dependent Ca2+-current, ω is the diffusion constant for intracellular Ca2+ and θ is a threshold for the rectified linear activation threshold of the membrane-voltage dependent Ca2+-current. The parameters used for adaptation throughout all models and regions are listed at the bottom of Table 1 and based on previous work (Hasselmo et al., 1995). A simulation of an isolated CA3 pyramidal cell with blocked recurrent connectivity exhibited similar adaptation to the time course of adaptation currents induced by opsin activation in voltage clamp (Figure S1; Boyden et al., 2005). In circuit, these parameters produced units that spent less than 50 ms above the high activity rate threshold and decayed from maximum activity rate to below threshold in under 30 ms representing a strong decrease to inferred firing rate.

CA3 ripple extension model input waveforms

We examined the effect of ramped and square pulse stimulation on prolonging an ongoing sequence, hereafter termed ‘replay’ activity. First, a forward propagating series of activation was elicited by delivery of a strong 20 ms-long cue pulse to the first node of the network. In a control simulation the replay was allowed to continue uninterrupted until relaxing to baseline (Figure 1B). In the waveform pulse simulations, a second waveform of varying duration and shape was delivered 150 ms after the cue offset to all excitatory neurons. In one set of simulations, the input duration was tested between 0 ms and 250 ms (Figure 1) while in another, the pulse delay after cue from 0 ms to 500 ms was examined (Figure S2). Figure 1C demonstrates the waveform shapes delivered to the network including Forward (FR), Backward (BR), and Double Ramp (DR). The degree of ramp, ‘Ramp Percentage’, varied linearly between 0% (square waveform) and 100% (triangular waveform) to characterize the response of the network under different waveforms. One set of simulations used an iso-maximum-amplitude (IMA) to match final ramp amplitude to the Square pulse amplitude (Padilla-Coreano et al., 2019). Taking advantage of the in silico model, a second set of simulations used a higher maximum amplitude to match the area-under-the-curve of the template Square pulse (Iso-Power, IP).

CA3 ripple extension model analysis

To examine the potential differences in CA3 sequence prolongation by different pulse shapes, a number of quantifications were used. First, the sequence length (from 0 to 15 nodes) was calculated for each simulation. The mean and 95% confidence interval of sequence length dependent on ramp degree was found using a 1000-sample shuffle with replacement of the sequence length outcome at different pulse durations (Figure 1F). Note that many combinations of ramp percentage and input duration parameters of the secondary pulse yielded full sequence recall.

Next, the time of each node reaching an activation threshold was obtained. An activation threshold of 10 was chosen as the first point where a high firing probability in each node occurred. Based on this, the Inter-Threshold-Interval (IThI) time between consecutive neurons quantifies the implied temporal sequence stability of the modeled sequence activation or ‘replay’ by the model. For each simulation an effect size (Cohen’SD; https://www.mathworks.com/matlabcentral/fileexchange/62957) was calculated between the waveform manipulation IThI time and the IThI of a control (Figure 1B; no secondary pulse) condition. Smaller effect size indicated a low impact of the optogenetic pulse on the temporal structure of the ongoing replay.

The effect of pulse shape on ripple extension behavior in the circuit was quantified as follows. First, the waveform that achieved the lowest temporal disruption (effect size) was calculated for each case (Figure 1G). The pulse duration at which this minimal effect size was also measured (Figure 1H). In cases with multiple lowest disruptions within the same ramp degree, the pulse with shortest input duration was used. To determine the generalized impact of waveform shape across tested input durations, a shuffle with replacement of effect size scores was calculated (Figure 1I). The mean effect size and a 95% confidence interval was calculated from a 1000-sample shuffle distribution for each waveform category.

The activity and timing disruption after the optogenetic pulse was analyzed in a similar fashion for the CA3 interneurons (Figure S3).

Optogenetic response variability model

To model sources of heterogeneity to the simulated optogenetic pulse, three additions to the model were made. First, the position of the fifteen pyramidal units in the replay extension model was uniformly randomly distributed in a 0.5 mm by 0.5 mm by 0.1 mm area, analogous to the dimensions of a rodent CA3 pyramidal layer (Hussein and George, 2009; Jinno and Kosaka, 2010; Öz and Saybaşılı, 2017). A point positioned in the center of the layer, 0.2 mm away, was used as the laser source of maximum irradiance (Ramirez et al., 2013). Irradiance decay as a function of distance from source was fit to the brain tissue using light transmission data available at https://web.stanford.edu/group/dlab/cgi-bin/graph/chart.php. Two different curves for 8 mW and 10 mW sources were fit. A theoretical irradiance value reaching each randomly positioned pyramidal unit was determined using this function and the 3D distance from the laser source. A gain constant between 0 and 1 was applied to simulated nodes distal to the light source receiving less than 5 mW/mm2 (Figure S4A; Zhang et al., 2006). The gain constant for all closer nodes receiving greater than 5 mW/mm2 was thresholded at 1 (Figure S4A). The optogenetic input pulse strength was multiplied by this gain factor.

Second, protein expression efficiency was modeled to address inefficient translation of light into membrane voltage changes when opsins are under or over-expressed (Gradinaru et al., 2008; Zhao et al., 2008). For each pyramidal unit, an expression efficiency value was drawn from the lower half of a random normal distribution center at 100% efficiency (Figure S5A). Two different variances, σ0.05 or σ0.1, were used to model a better or worse optimized opsin expression distribution, respectively. The randomly drawn constant of each neuron was applied to the optogenetic drive for that neuron during simulations.

Third, a white noise of varying amplitude was applied to the membrane voltage of each neuron at every time step (Figure S6A). An amplitude of +/ 0.1 was used to approximate low levels of noise. An amplitude of +/− 0.5 was used to approximate high levels of noise. This resulted in neurons being more or less depolarized relative to resting potential at the time of the replay cue onset or the optogenetic pulse onset.

Each of these 3 sources of variability were modeled alone (Figures S4–S6) or combined (Figure 2). In this case, the opsin efficiency was multiplied by the received irradiance gain factor as a combined gain applied to the optogenetic drive for each neuron. For the combined simulations, the higher light power (10 mW at source), lower opsin efficiency variability (σ0.05) and lower voltage noise amplitude (+/− 0.1) were used. For all parameter space search simulations, a different pseudorandom kernel was used for each individual simulation with matched kernel for control vs optogenetic pulse.

CA1 replay extension model

To test the outcome of optogenetic pulses extending ongoing replay in CA1, we implemented a CA1 model of 15 Pyramidal and 15 Interneurons connected to the CA3 model described above (Figure 3A). The CA3 pyramidal cells sent feedforward excitation with symmetric narrow spread to CA1 interneurons and asymmetric spread to CA1 pyramidals (Figure S6). CA1 pyramidals strongly excited CA1 interneurons which imposed strong inhibition back to the pyramidals. Finally, CA1 pyramidals were very weakly connected to one another but not themselves (Deuchars and Thomson, 1996; Shepherd, 2004; Malerba & Bazhenov, 2019). The membrane potential of CA1 Pyramidals and Interneurons was updated as follows:where WZ, WQ, ZZ, ZQ and QZ are weight matrices determining the connectivity between neurons (Figure S7) and r is the index of the pre-synaptic CA3 pyramidal unit. The weights from CA3 to CA1 pyramidal nodes to were biased such that earlier nodes received stronger connections. Similarly, the innervation of CA1 interneurons was weaker from CA3 for earlier CA1 interneurons. The adaptation variables used are not included for brevity, but were identical to those in Equation 1 and applied to both CA3 and CA1 pyramidal units.

The same optogenetic pulse waveforms were applied to the CA1 model 150 ms after a 20 ms cue to the CA3 network. Similar shuffling procedures for analysis of the impact of ramp degree on temporal disruption and sequence length were used as in the CA3 replay extension model.

CA3 learning model

Finally, a synaptic learning rule was implemented to examine the impact of pulse shape on the formation of synaptic weight matrices. The same network as in the isolated CA3 model (15 excitatory, 15 inhibitory) was employed with two exceptions (Figure 4A). First, the weight matrix, W, was instantiated as uniform zeroes indicating no pre-existing connectivity between the nodes. Second, an acetylcholine (Ach) term was added to modify the strength of local recurrent connections and network plasticity, consistent with prior experimental observations and computational models (Hasselmo et al., 1995, 1998; Hasselmo 2006). Simulations were divided into a Learn Phase and a Test Phase.

CA3 learning model - learn phase

During the Learn Phase, an input pattern composed of 15 elements was presented as a sequence through time, one element for each excitatory node of the network (Figure 4A, right). The equation estimating activity of pyramidals was modified to include an Acetylcholine (Ach) component as follows:where ψ is the level of Ach present in the system, between 0 and 1. During the Learn Phase ψ was set to 0.9 to represent strong neuromodulatory attenuation of local recurrent CA3 collaterals and prevent over-generalization of learning. The adaptation variables used are not included for brevity, but matched those in Equation 1 and used the same parameters as in Table 1.

A learning rule governed the modification of synaptic weight between co-active nodes:where W is a matrix of weights with pre-synaptic columns and post-synaptic rows, ε is the learning rate for the system, the presence of Ach (ψ) modulates the plasticity of the synapses, B acts as a potentiation ceiling, and P is the instantaneous activity vector of the CA3 excitatory nodes. Table 3 summarizes the values used for each parameter in this model, modified slightly from the Ripple Extension model to enable better 1-trial pattern learning.

During the Learn Phase each pattern element was set at a fixed 80 ms duration, chosen based on a set of simulations varying learn pulse duration at different overlap values (Figure S8). As in previous simulations, the ramp degree was varied between 0% (Square) and 100% ramped for Forward, Backward, and Double Ramp, and Iso-Max-Amplitude and Iso-Power pulses. Another set of simulations manipulated the duration of pattern element overlap (OL), between 0% (no overlap of consecutive elements) and 100% overlap (simultaneous pulses across nodes).

CA3 learning model - test phase

During the Test Phase a square, 20 ms cue pulse was delivered to the first node to elicit replay-like sequential node activations and assess the strength of retrieval. The final state of the weight matrix W after the Learn Phase was used as the auto-recurrent architecture representing the outcome of single-trial encoding of pattern information. The same equations and parameters as in the Learn Phase were used, with the exception that (Ach level) was set to 0 to prevent further modifications to the weight matrix and allow full auto-recurrent excitation (Hasselmo, 2006).

CA3 learning model analysis

To quantify the recall capacity of the network under different optogenetic waveform paradigms, the same threshold crossing metric as in the CA3 Ripple Extension Model was used. Any nodes during the Test Phase that crossed an activity threshold of 10 were counted toward participating in the sequence. The IThI was used as a metric of sequence replication. For each simulation varying ramp shape, percentage and pattern overlap, an effect size (Cohen’SD) was calculated between the test phase mean IThI and a control IThI.

The control data was generated using a pre-allocated, symmetric weight matrix with a max auto-recurrent strength equal to 0.035. Connectivity to nearest neighbors was 0.0255 (Figure 4B). All other parameters were equal to those in the learning simulations. This control was chosen because it elicited a seven-node sequence length, allowing for potentially longer or shorter sequences to be achieved by varying the learning input paradigm. This sequence length variability is consistent with in vivo work on shorter and longer sequence replays (Fernandez-Ruiz et al., 2019). This control also yielded high similarity (low-effect size) relative to the learned replays across ramp percentages (Figure 4I).

Low overlap percentage generally resulted in suprathreshold activation of only the first node, elicited by the cue pulse directly. Very high overlap percentage generally resulted in over-saturation of the weight matrix, causing excessive node activation in the Test Phase. We therefore restricted the following analysis to a band of overlap percentages which resulted in a dynamic range of sequence lengths. The cutoff for excessively high activity was set at 100 for any single unit in the network.

We quantified the effects of waveform shape on sequence learning a number of ways. First, the lowest overlap percentage that elicited the smallest effect size was calculated for each waveform shape over each ramp degree. For each Waveform Class, a 500-sample shuffle with replacement was performed on the effect size, mean IThI and sequence length data. In order to examine the contribution of waveform shape independent of learning sequence timing, each shuffle maintained ramp percentage identity while shuffling overlap percentage.

Quantification and statistical analysis

All models were simulated with custom code on a PC desktop running Windows 10 using Matlab 2019a (Mathworks Inc.). Analyses and statistics were run using the same resources. Ramp percentage and memory metric correlations were performed using Pearson’s r with α = 0.05 and are reported in the Results text. Shuffle permutations of N = 1000-samples for replay extension models (Figures 1, 2, and 3) and N = 500-samples for the learning model (Figure 4) are displayed with a 95% confidence interval and are described further in each subsection of the method details.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Parameter variables and shuffled data	This paper	Zenodo: https://doi.org/10.5281/zenodo.6415338
Code	This paper	Zenodo: https://doi.org/10.5281/zenodo.6415338
Software and algorithms
Matlab 2019a	MathWorks	https://www.mathworks.com