Offline memory replay in recurrent neuronal networks emerges from constraints on online dynamics

During spatial exploration, neural circuits in the hippocampus store memories of sequences of sensory events encountered in the environment. When sensory information is absent during “offline” resting periods, brief neuronal population bursts can “replay” sequences of activity that resemble bouts of sensory experience. These sequences can occur in either forward or reverse order, and can even include spatial trajectories that have not been experienced, but are consistent with the topology of the environment. The neural circuit mechanisms underlying this variable and flexible sequence generation are unknown. Here we demonstrate in a recurrent spiking network model of hippocampal area CA3 that experimental constraints on network dynamics such as population sparsity, stimulus selectivity, rhythmicity, and spike rate adaptation enable additional emergent properties, including variable offline memory replay. In an online stimulus-driven state, we observed the emergence of neuronal sequences that swept from representations of past to future stimuli on the timescale of the theta rhythm. In an offline state driven only by noise, the network generated both forward and reverse neuronal sequences, and recapitulated the experimental observation that offline memory replay events tend to include salient locations like the site of a reward. These results demonstrate that biological constraints on the dynamics of recurrent neural circuits are sufficient to enable memories of sensory events stored in the strengths of synaptic connections to be flexibly read out during rest and sleep, which is thought to be important for memory consolidation and planning of future behavior.


Introduction
In mammals, the hippocampus is a brain region involved in the storage and recall of spatial and episodic memories. As an animal explores a spatial environment, different 25 subpopulations of hippocampal neurons known as "place cells" are selectively activated at different positions in space, resulting in sequences of neuronal spiking that are on the seconds-long timescale of locomotor behavior (O'Keefe and Conway, 1978).
The synchronous firing of these sparse neuronal ensembles is coordinated by population-wide oscillations referred to as theta (~4-10 Hz) and gamma (~30-100 Hz) 30 rhythms (Colgin, 2016). Within each cycle of the theta rhythm (~125 ms), the spiking of active neurons is organized into fast timescale sequences such that neurons selective for just-visited positions spike first, then neurons selective for the current position, and finally neurons selective for the next and future positions spike last in the sequence (Drieu and Zugaro, 2019;Foster and Wilson, 2007). These order-preserving fast 35 timescale "theta sequences" are thought to be involved in planning and learning of event order through associative synaptic plasticity (Jensen et al., 1996;Kay et al., 2020).
When an animal stops running, theta and gamma oscillations decrease, and neuronal circuits in the hippocampus instead emit intermittent synchronous bursts of 40 activity that are associated with high-frequency oscillatory activity detectable in local field potential recordings in hippocampal area CA1. These ~100-200 ms long events are referred to as "sharp wave-ripples" (SWRs) (Colgin, 2016), and they occur during non-locomotor periods of quiet wakefulness, during reward consumption, and during slow-wave sleep, when sensory information about the spatial environment is reduced 45 or absent. During SWRs, sparse subsets of neurons are co-activated, with a tendency for neurons that fire sequentially during exploratory behavior to also fire sequentially during SWRs, either in the same order, in reverse, or with a mixture of both directions (Davidson et al., 2009;Pfeiffer, 2020;Stella et al., 2019;Wu and Foster, 2014). The hippocampus can also activate sequences of place cells during SWRs that correspond 50 to possible paths through the environment that have not actually been experienced, suggesting a possible role for offline sequence generation in behavioral planning (Gupta et al., 2010;Igata et al., 2021;Ólafsdóttir et al., 2015Wu and Foster, 2014). Manipulations that disrupt neuronal activity during SWRs result in deficits in memory recall (Girardeau et al., 2009;Jadhav et al., 2012), supporting a hypothesis 55 that offline reactivation of neuronal ensembles during SWRs is important for the maintenance and consolidation of long-term memories (Buzsáki, 1989;Joo and Frank, 2018).
Hippocampal SWRs are thought to be generated by the synchronous firing of subpopulations of neurons in the CA2 and CA3 regions of the hippocampus (Csicsvari 60 et al., 2000;Oliva et al., 2016), which are characterized by substantial recurrent feedback connectivity (Duigou et al., 2014;Guzman et al., 2016;Okamoto and Ikegaya, 2019). Recurrent networks have long been appreciated for their ability to generate rich internal dynamics (Amit and Brunel, 1997), including oscillations (Ermentrout, 1992). It has also been shown that associative plasticity at recurrent connections between 65 excitatory neurons can enable robust reconstruction of complete memory representations from incomplete or noisy sensory cues (Amit and Brunel, 1997;Griniasty et al., 1993;Guzman et al., 2016;Hopfield, 1982;Marr, 1971;Treves and Rolls, 1994). However, this "pattern completion" function of recurrent networks requires that similarly tuned neurons activate each other via strong synaptic 70 connections, resulting in sustained self-activation, rather than sequential activation of neurons that are selective for distinct stimuli (Lisman et al., 2005;Pfeiffer and Foster, 2015;Renno-Costa et al., 2014). Previous work has shown that, in order for recurrent networks to generate sequential activity, some mechanism must be in place to "break the symmetry" and enable spread of activity from one ensemble of cells to another 75 (Sompolinsky and Kanter, 1986;Tsodyks et al., 1996). During spatial navigation, feedforward sensory inputs carrying information about the changing environment can provide the momentum necessary for sequence generation. However, during hippocampal SWRs, sensory inputs are reduced, and activity patterns are thought to be primarily internally generated by the recurrent connections within the hippocampus. 80 In this study we use a computational model of hippocampal area CA3 to investigate the synaptic, cellular, and network mechanisms that enable flexible offline generation of memory-related neuronal sequences in the absence of ordered sensory information.
A number of possible mechanisms for sequence generation in recurrent networks have been proposed: 85 1) Winner-take-all network mechanism (Almeida et al., 2007(Almeida et al., , 2009aLisman and Jensen, 2013): Within this framework, the subset of excitatory neurons receiving the most strongly weighted synaptic inputs responds first upon presentation of a stimulus. This active ensemble of cells then recruits feedback inhibition via local interneurons, which in turn prevents other neurons from firing for a brief time 90 window (e.g. the ~15 ms duration of a single gamma cycle). This highlights the important roles that inhibitory neurons play in regulating sparsity (how many cells are co-active), selectivity (which cells are active), and rhythmicity (when cells fire) in recurrent networks (Almeida et al., 2009b;Rennó-Costa et al., 2019;Stark et al., 2014;Stefanelli et al., 2016). However, while oscillatory feedback inhibition 95 provides a network mechanism for parsing neuronal sequences into discrete elements, additional mechanisms are still required to ensure that distinct subsets of excitatory neurons are activated in a particular order across successive cycles of a rhythm (Lisman et al., 2005;Ramirez-Villegas et al., 2018).
2) Heterogeneous cellular excitability (Luczak et al., 2007;Stark et al., 2015): If the 100 intrinsic properties of neurons in a network are variable and heterogeneous, when a stimulus is presented, those neurons that are the most excitable will fire early, while neurons with progressively lower excitability will fire later, resulting in sequence generation. This mechanism can explain the offline generation of stereotyped, unidirectional sequences, but cannot account for variable generation 105 of sequences in both forward and reverse directions.
3) Asymmetric distributions of synaptic weights (Sompolinsky and Kanter, 1986;Tsodyks et al., 1996): During learning, if changes in synaptic weights are controlled by a temporally asymmetric learning rule, recurrent connections can become biased such that neurons activated early in a sequence have stronger 110 connections onto neurons activated later in a sequence (Levy, 1989;Malerba and Bazhenov, 2019;McNaughton and Morris, 1987;Reifenstein et al., 2021). This enables internally generated activity to flow along the direction of the bias in synaptic weights. While this mechanism accounts for offline replay of specific sequences in the same order experienced during learning, it cannot account for 115 reverse replay or the flexible generation of non-experienced sequences (Gupta et al., 2010;Igata et al., 2021;Ólafsdóttir et al., 2015Wu and Foster, 2014). 4) Cellular or synaptic adaptation: It has also been proposed that short-term adaptation of either neuronal firing rate (Itskov et al., 2011;Treves, 2004) or synaptic efficacy (Romani and Tsodyks, 2015) can enable neuronal sequence 120 generation in recurrent networks without asymmetric synaptic weights. According to this scheme, recently-activated neurons initially recruit connected partners with high efficacy, but continued spiking results in either a decrease in firing rate, or a decrease in the probability of neurotransmitter release. This causes connections to fatigue over time, and favors the sequential propagation of activity to more 125 recently activated cells. These mechanisms do allow for the stochastic generation of neuronal sequences in either the forward or reverse direction, though they do not prescribe which or how many neurons will participate in a given replay event.
In this study, we sought to understand how neuronal sequence generation in hippocampal area CA3 depends on the structure and function of the underlying 130 network. To do this, we constructed a computational neuronal network model comprised of recurrently connected excitatory and inhibitory spiking neurons, and tuned it to match experimental constraints on the spiking dynamics of CA3 during spatial navigation, including sparsity, selectivity, rhythmicity, and spike rate adaptation.
We then analyzed the direction and content of neuronal sequences generated both 135 "online" during simulated navigation, and "offline" during simulated rest. We found that when the network was driven by ordered sensory information in the online state, it generated forward-sweeping "theta sequences" that depended on the structure of recurrent connectivity in the network. In the offline state driven by noise, the network generated heterogenous memory replay events that moved either in forward, reverse, 140 or mixed directions, and depended on network sparsity and rhythmicity, and neuronal stimulus selectivity and spike-rate adaptation. Model degeneracy analysis and network perturbations indicated that offline memory replay does not occur in networks with disrupted recurrent connectivity, or in networks lacking sparsity, selectivity, rhythmicity, or spike rate adaptation. Finally, when particular spatial locations were 145 over-represented by the network, as occurs in the hippocampus at sites of reward (Lee et al., 2006;Turi et al., 2019;Zaremba et al., 2017), memory replay events were biased towards trajectories that included those salient positions (Gillespie et al., 2021;Ólafsdóttir et al., 2015;Singer and Frank, 2009).

Results
To investigate how sequential activity in the hippocampus generated "online" during spatial exploration can be recapitulated "offline" in the absence of sensory cues, we constructed a simple spiking neuronal network model of rodent hippocampal area CA3 (Materials and Methods). Neural circuits in the hippocampus and cortex typically 155 comprise a majority of excitatory neurons that project information to downstream circuits, and a minority of primarily locally connected inhibitory interneurons. We included populations of excitatory (1000) and inhibitory (200) neurons in proportion to experimental observations (Pelkey et al., 2017;Tremblay et al., 2016) (Figure 1A). Cell models were single-compartment, integrate-and-fire neurons with saturable, 160 conductance-based excitatory and inhibitory synapses (Carnevale and Hines, 2006;Izhikevich, 2007;Izhikevich and Edelman, 2008). Excitatory neurons were endowed with spike-rate adaptation to support punctuated bursting behavior during theta oscillations (O'Keefe and Recce, 1993;Scharfman, 1993), and inhibitory neurons exhibited fast-spiking dynamics to sustain continuous high frequency firing during 165 from one example network instance). Cells in each population are sorted by the location of maximum firing. (F) Average stimulus selectivity of each cell population. Trial-averaged activity of each cell was centered around the location of maximum firing, and then averaged across cells. (G) The average activity of each population on a single trial (top row) was bandpass filtered in the theta (middle row) and gamma (bottom row) frequency bands. Colored traces show filtered signals (theta: green, gamma: purple).

180
Traces derived from one example network instance. (H) Power spectrum of average population activity indicates dominant frequency components in the theta and gamma bands (one-sided paired t-tests: theta: E vs. FF,p=0.00001;I vs. FF,p<0.00001;gamma: E vs. FF,p<0.00001;I vs. FF,p<0.00001). In (C), (D), (F), and (H), data were first averaged across 5 trials per network instance. Mean (solid) ± SEM (shading) were computed across 5 independent instances of each network model. p-values reflect FDR 185 correction for multiple comparisons.
gamma oscillations (Csicsvari et al., 2003;Ylinen et al., 1995)  To simulate the sensory experience of locomotion in a spatial environment, we 190 provided both excitatory and inhibitory neurons with external afferent inputs from a population of 1000 excitatory neurons, each of which was selectively activated at distinct but overlapping positions within a simulated circular track that took 3 seconds to traverse ( Figures 1A -1C, Materials and Methods). Recurrent connections within and between excitatory and inhibitory cell populations were also included ( Figure 1A), 195 as they are hallmark features of hippocampal area CA3, and have been shown to support rich network dynamics (Renno-Costa et al., 2014;Stark et al., 2014).
Specifically, inhibitory feedback connections have been shown to regulate the number of simultaneously active neurons (sparsity) (Stefanelli et al., 2016), and to contribute to the generation of theta and gamma network oscillations (Bezaire et al., 2016;Geisler et 200 al., 2005;Rennó-Costa et al., 2019;Stark et al., 2014;Wang, 2010). Plastic excitatory connections between excitatory neurons have long been implicated in stimulus selectivity and the storage and recall of memories (Almeida et al., 2007;Hopfield, 1982;Lisman and Jensen, 2013). It has been proposed that strong connections between ensembles of co-active neurons could arise through a combination of biased 205 connectivity during brain development (Buzsáki et al., 2021;Dragoi and Tonegawa, 2013;Farooq and Dragoi, 2019;Grosmark and Buzsáki, 2016), and experience-driven synaptic plasticity during learning (Bittner et al., 2015Brunel and Trullier, 1998;Káli and Dayan, 2000;Milstein et al., 2020;O'Neill et al., 2008). While here we did not simulate these dynamic processes explicitly, we implemented the structured 210 connectivity that is the end result of these processes by increasing the strengths of synaptic connections between excitatory cells that share overlapping selectivity for spatial positions in the environment (Table 1, Supplementary Figure S1B, Materials and Methods) (Arkhipov et al., 2018).
Despite the relatively simple architecture of this network model, a wide range of 215 networks with distinct dynamics could be produced by varying a number of parameters, including 1) the probabilities of connections between cell types (Káli and Dayan, 2000), 2) the kinetics and strengths of synaptic connections between cell types (Brunel and Wang, 2003), and 3) the magnitude of the above-mentioned increase in synaptic strengths between neurons with shared selectivity (Brunel, 2016;Dorkenwald 220 et al., 2019). To calibrate the network model to produce dynamics that matched experimentally-derived targets, we performed an iterative stochastic search over these parameters, and optimized the following features of the activity of the model network: This procedure identified a model with dynamics that met all of the above constraints. Given sparse and selective feedforward inputs during simulated navigation (Figures 1B and 1C), the excitatory neurons in the network responded with a fraction of 230 active cells ( Figure 1D) and with average firing rates comparable to the those of the feedforward input population ( Figure 1E). The majority of inhibitory neurons were activated continuously ( Figures 1C and 1D) at high firing rates ( Figure 1E). While excitatory neurons received random connections from feedforward afferents and from other excitatory neurons with heterogeneous spatial tuning, excitatory cells exhibited a 235 high degree of spatial selectivity ( Figures 1C and 1F). This selective increase in firing rate at specific spatial locations within the "place field" of each excitatory neuron was supported by enhanced synaptic connection strengths between excitatory neurons with overlapping tuning (Supplementary Figure S1B). While substantial background excitation occurred in all cells at all spatial positions, firing outside the place field of 240 each cell was suppressed by sufficiently strong inhibitory input (Bittner et al., 2015;Grienberger et al., 2017). Interestingly, inhibitory neurons also exhibited spatial selectivity, albeit to a weaker degree and with a higher background firing rate ( Figures   1C and 1F). This feature of the network dynamics was an emergent property that was not explicitly designed or optimized. While excitatory connections onto inhibitory cells 245 were random and not weighted according to shared selectivity (Supplementary Figure   S1B), the total amount of excitatory input arriving onto individual inhibitory cells fluctuated across spatial positions, and predicted a small degree of spatial selectivity (Supplementary Figure S1C). Inhibitory inputs received by inhibitory cells reduced their average activity, effectively enabling fluctuations in excitation above the mean to stand 250 out from the background excitation (Supplementary Figure S1C and S1D). This mechanism of background subtraction by inhibitory synaptic input may explain the partial spatial selectivity previously observed in subpopulations of hippocampal inhibitory neurons (Ego-Stengel and Wilson, 2007;Geiller et al., 2020;Hangya et al., 2010;Marshall et al., 2002;Wilent and Nitz, 2007). 255 The tuned network model also exhibited oscillatory synchrony in the firing of the excitatory and inhibitory neuron populations, despite being driven by an asynchronous external input ( Figures 1G and 1H). The requirement that the network self-generate rhythmic activity in the theta band constrained recurrent excitatory connections to be relatively strong, as this input provided the only source of rhythmic excitation within the 260 network (Supplementary Figure S1E). Interestingly, as the firing rates of inhibitory cells increased within each cycle of the theta rhythm, their synchrony in the gamma band increased, resulting in an amplitude modulation of gamma paced at the theta frequency ( Figure 1G and Supplementary Figure S1F). This "theta-nested gamma" is a well-known feature of oscillations in the hippocampus (Soltesz and Deschenes, 1993;265 Ylinen et al., 1995), and here emerged from fundamental constraints on dual band rhythmicity without requiring additional mechanisms or tuning.

Position decoding reveals "theta sequences" during simulated navigation
Next, we analyzed neuronal sequence generation within the network during simulated 270 navigation. First, we simulated multiple trials and computed trial-averaged spatial firing rate maps for all neurons in the network ( Figure 1C). We then used these rate maps to perform Bayesian decoding of spatial position given the spiking activity of all cells in the network from individual held-out trials not used in constructing the decoding template ( Figure 2A, Materials and Methods) (Davidson et al., 2009;Zhang et al., 1998). 275 For the population of feedforward excitatory inputs, the underlying spatial firing rates were imposed, and the spikes of each cell were generated by sampling from an inhomogeneous Poisson process. Thus, decoding position from the activity of this population served to validate our decoding method, and indeed simulated position  Figure 1, in which recurrent excitatory connections between E cells were structured such that neurons with shared selectivity have elevated synaptic weights. Top row: spike times of all neurons in each cell population on a single trial of simulated "online exploration" are marked. A separate set of 5 trials was used to 285 construct a spatial firing rate template for each neuron (shown in Figure 1C). Cells in each population are sorted by the location of maximum average spatial firing rate. Bottom row: the spatial firing rate templates for all neurons were used to perform Bayesian decoding of spatial position from the single trial spiking data shown in the top row. For each cell population, the likelihood of each spatial position in each time bin (20 ms

295
The absolute value of decoded position error is expressed as a fraction of the track length (  Materials and Methods). Interestingly, we found that position could also be accurately decoded from the moderately spatially-tuned activity of inhibitory cells in the network (Figures 2A and 2D), and that the spiking activity of the inhibitory population was also organized into theta sequences (Figures 2A and 2E).
A number of possible mechanisms have been proposed to account for theta 330 sequence generation in vivo, including synaptic, cell-intrinsic, and network-level mechanisms (Chadwick et al., 2015(Chadwick et al., , 2016Drieu and Zugaro, 2019;Foster and Wilson, 2007;Grienberger et al., 2017;Kang and DeWeese, 2019;Mehta et al., 2002;Skaggs et al., 1996). That theta sequences in the model emerged in both excitatory and inhibitory neuron populations implicates recurrent interactions within the network 335  Figure 2E). These results indicate that structure in the synaptic strengths 355 of recurrent excitatory connections is required for the generation of fast timescale (~125 ms) neuronal sequences when network activity is driven by behavioral timescale (> 1 s) sequences of sensory inputs, as occurs during spatial exploration.

Emergence of offline memory replay 360
The above results show that the same network structure that enables population dynamics in CA3 to exhibit sparsity, selectivity, and rhythmicity also supports neuronal sequence generation in the online state when ordered sensory information is present.
We next sought to understand how neuronal sequences consistent with the sensory environment are generated offline when sensory inputs are reduced. To mimic the

405
( Figures 3A and 3B). We then used the same decoding templates as above, constructed from the trial-averaged activity during simulated run, to decode spatial position from spiking activity during these transient offline events (Materials and Methods).
Given that the place field locations of the stimulated neurons in the feedforward 410 input population were heterogeneous and unordered, the spatial positions decoded from their spiking were typically discontiguous across adjacent temporal bins ( Figures   3A, 3B and 3F). This input pattern evoked spiking in sparse subsets of both the excitatory and inhibitory populations in the network ( Figures 3A and 3B). In contrast with the feedforward population, the activity evoked in excitatory neurons was 415 structured such that neurons with nearby place field locations spiked in adjacent temporal bins, resulting in decoded spatial trajectories that were continuous (Figures 3A,3B and 3F). Inhibitory neuron activity during these events was organized into highfrequency oscillations (Figures 3A,3B and 3G). This procedure was repeated to produce thousands of offline events evoked by stimulation of different random 420 ensembles of inputs (Materials and Methods). Across these events, each position along the track was decoded with equal probability ( Figure 3C). For each event, the length and mean velocity of the decoded trajectory was calculated from the differences in decoded positions between adjacent bins ( Figures 3D and 3E). A mean velocity of zero corresponds to events with equal steps in the forward and reverse directions, while 425 positive velocities correspond to net forward-moving trajectories, and negative velocities correspond to net backwards-moving trajectories. While the trajectories decoded from the random feedforward input population were comprised of large, discontiguous steps that traced out large path lengths with an average velocity near zero, the excitatory neuron population generated shorter, more continuous sequences 430 that progressed in either the forward or reverse directions ( Figures 3D -3F). These trajectories on average covered ~0.5 the length of the track in the short (~150 ms) duration of the offline event. Compared to the run trajectory, which took 3 seconds to cover the full track length, this corresponded to a ~10-fold temporal compression ( Figure 3E), similar to experimental data (Davidson et al., 2009). Spatial trajectories 435 decoded from the inhibitory neuron population were intermediate in length, but with little forward or reverse momentum, similar to the feedforward inputs. However, the inhibitory cells exhibited high-frequency synchrony (Figures 3A, 3B and 3G), similar to experimentally recorded CA3 interneurons during hippocampal SWRs (Csicsvari et al., 2000;Tukker et al., 2013). 440 These data demonstrate that random, unstructured input can evoke sequential activity in a CA3-like recurrent spiking network, with sequences corresponding to forward, reverse, or mixed direction trajectories through an experienced spatial environment. This self-generated memory-related activity implicates information stored in the synaptic weights of the recurrent connections within the network as being 445 important for offline replay of experience. However, in most previous models, sequence generation was unidirectional, and was enabled by an asymmetric bias in the strengths of recurrent connections such that neurons encoding positions early in a sequence formed stronger synapses onto neurons encoding later positions (Levy, 1989;Malerba and Bazhenov, 2019;McNaughton and Morris, 1987;Reifenstein et al., 450 2021;Sompolinsky and Kanter, 1986;Tsodyks et al., 1996). In contrast, the current network flexibly generated sequences in forward, reverse, or mixed directions, and had symmetric recurrent connections such that synaptic strengths between pairs of excitatory neurons depend only on overlapping selectivity, not on sequence order.
Does sequence generation in the present network still depend on recurrent 455 connectivity? To test this, we first verified that including an asymmetric bias in the strengths of excitatory connections produced offline replay events that were biased towards forward sequences (Supplementary Figure S3). We next analyzed the sequence content of offline events generated in the variants of the network model with random ( Figures 3H -3N Figures S4A -S4F). Still, these networks exhibited high-465 frequency oscillatory synchrony during these offline events ( Figure 3N and Supplementary Figure S4G).

Exploration of model diversity and degeneracy
The above results strongly supported the hypothesis that recurrent connectivity is 470 important for offline sequence generation. During optimization of each of the alternative network model configurations shown above to meet the multiple objectives of sparsity, selectivity, and rhythmicity, we evaluated 30,000 variants of each model with different parameters (Materials and Methods). For each model configuration, the model with the lowest overall multi-objective error was chosen as the "best" model for further analysis, 475 as shown in Figures 2 and 3 and Supplementary Figures S1, S2 and S4. However, we noted that the parameter values that specified these "best" models were variable across the multiple network configurations (Table 1). This raised the possibility that while many models with diverse parameters may produce networks with similar online dynamics (referred to as model degeneracy (Marder and Taylor, 2011)), perhaps only a 480 smaller subset of models would additionally support the emergence of offline sequence generation. The fact that sequence generation was not observed for either of the "best" models with disrupted recurrent excitatory synaptic weights could reflect an incomplete sampling of the parameter space and a failure to identify models that both meet the objective criteria for online dynamics and produce offline sequences. 485 Therefore, in order to explore the diversity and degeneracy of models evaluated during optimization, we devised a method to identify models that performed similarly with respect to multiple optimization objectives, but were specified by divergent sets of parameters (Materials and Methods). For each model configuration, all model variants evaluated during parameter optimization were sorted by their Euclidean distance from 490 the "best" model in the space of model parameters. This resulted in an error landscape (e.g. Figure 4A) in which models with similar parameters resulted in similar multiobjective error scores. We then identified models located at local minima in this error landscape, which formed a group of models that were distant from each other in parameter space, but similar to each other in terms of overall multi-objective error. We 495 termed a set of such models as a "Marder group" after pioneering work characterizing degeneracy in biological systems (Marder and Taylor, 2011). For each alternative network model configuration, we selected the 5 members of this "Marder group" with the lowest multi-objective error (labeled "best" and "M1" -"M4" in Figure 4A), and evaluated their network dynamics during simulations of both online exploration and 500 offline rest. We first verified that for all model configurations with and without structured recurrent excitatory connectivity analyzed above (Figures 1 -3 and Supplementary Figures S1, S2 and S4), model variants within a "Marder group" exhibited considerable diversity across model parameters ( Figures 4B and 4C), and met all objective criteria for population sparsity ( Figure 4D), neuronal stimulus 505 selectivity ( Figure 4E), and rhythmogenesis in the theta and gamma bands (Figures 4F and 4G). However, only model variants with synaptic weights structured by shared and degeneracy, for each network model configuration, a subset of model variants termed a "Marder group" were selected based on their large distance from each other in the space of parameters, but their similar performance with respect to multiple optimization objectives (see Materials and Methods). This selection procedure is illustrated here for the model with random E ← E weights as an example. The 5 "Marder group" members with the lowest multi-objective error (labeled "best" and "M1" -"M4") were 515 selected for further evaluation. (B) For the network model configuration with structured E ← E weights, the range of parameter values across 5 distinct "Marder group" models are shown. (C) Same as (B) for the model with random E ← E weights. In (D -I), features of the simulated network dynamics produced by distinct model variants within a "Marder group" are compared across network model configurations. Each data point (grey circles) depicts one "Marder group" model. (D) Spatial selectivity of the excitatory 520 neuron population during simulated "online exploration" is computed as a ratio of maximum to mean activity (two-sided t-tests vs. data from model with structured E ← E weights: Random E ← E weights: p=1; Shuffled E ← E weights: p=0.00224). (E) The fraction of the excitatory neuron population that is synchronously active during simulated "online exploration" is shown (two-sided t-tests vs. data from model with structured E ← E weights: Random E ← E weights: p=1; Shuffled E ← E weights: 525 p=0.12079). (F) Gamma rhythmicity of the excitatory neuron population is computed as the integrated power spectral density in the gamma frequency band (two-sided t-tests vs. data from model with structured E ← E weights: Random E ← E weights: p=0.03552; Shuffled E ← E weights: p=0.16364). (G) Theta rhythmicity of the excitatory neuron population is computed as the integrated power spectral density in the theta frequency band (two-sided t-tests vs. data from model with structured E ← E 530 weights: Random E ← E weights: p=0.00272; Shuffled E ← E weights: p=0.01944). (H) Theta sequence score (see Figure 2 and Materials and Methods) (two-sided t-tests vs. data from model with structured E ← E weights: Random E ← E weights: p=0.02817; Shuffled E ← E weights: p=0.01473). (I) Fraction of events during simulated "offline rest" that met criterion for sequences consistent with continuous spatial trajectories (see   stimulus selectivity exhibited theta sequences during online run ( Figure 4H) and generated offline sequences consistent with continuous spatial trajectories ( Figure 4I). This analysis demonstrated that generation of memory-related neuronal sequences by 545 recurrent networks requires that information about the topology of the sensory environment is stored in the strengths of recurrent excitatory connections between excitatory neurons.

Constraints on online sparsity, selectivity and rhythmicity enable offline memory 550 replay
Our above findings suggest that experimental constraints on the online dynamics of hippocampal area CA3 during spatial exploration are sufficient to enable the emergence of offline memory replay. We next sought to determine whether all or only a subset of these constraints were required for generation of memory-related sequences. 555 To determine the importance of rhythmicity, we removed the optimization criteria that excitatory and inhibitory neuron populations synchronize in the theta and gamma bands, and instead added an objective to minimize power density across the full frequency spectrum (Supplementary Figure S5E). Following optimization, this alternative network model exhibited reduced rhythmicity, but still met objectives related 560 to sparsity and selectivity ( Supplementary Figures S5A -S5D). However, when challenged with random stimuli during simulated offline rest, this model with suppressed rhythmicity failed to generate continuous forward or reverse sequences (Supplementary Figures S5F -S5L). Rather, spatial trajectories decoded from offline population activity contained large discontiguous jumps in position, and most events 565 had zero net velocity in either the forward or reverse direction. This indicated that the same reciprocal interactions between excitatory and inhibitory neurons that support rhythmogenesis in the theta and gamma bands also contribute to the sequential organization of neuronal activity during offline memory replay.  Figure S6E). During simulation of offline rest, this network generated highly synchronous population bursts that tended to either hover at one decoded position, or make large discontiguous jumps between 580 positions (Supplementary Figure S6F -S6K). These results suggest that the network connectivity parameters that support highly sparse and selective neuronal activity in the online stimulus-driven state, also enable sparse reactivation of neuronal sequences in the offline state. We also repeated the model degeneracy analysis described above ( Figure 4) for multiple model variants with compromised sparsity, selectivity or 585 rhythmicity, which corroborated these findings (Supplementary Figure S8).

Role of neuronal spike rate adaptation in forward and reverse offline memory replay
Above we showed that structure in the synaptic connectivity of the CA3 network is 590 important for neuronal sequence generation. However, unlike previous models of sequence generation where asymmetry in connection strengths biased the direction of neuronal sequences (Levy, 1989;Malerba and Bazhenov, 2019;McNaughton and Morris, 1987;Reifenstein et al., 2021;Sompolinsky and Kanter, 1986;Tsodyks et al., 1996), here synaptic connectivity was symmetric, and yet variable neuronal sequences 595 were flexibly generated in both forward and reverse directions ( Figure 3A, 3B and 3E).
If this symmetric connectivity enables recurrent networks to generate either forward or backward steps, what "breaks" this symmetry and produces sequences that make net progress in either the forward or backward direction? We next wondered if directionality of offline sequences in our network model was facilitated by our choice of 600 "bursty" excitatory cell model, which exhibited spike-rate adaptation ( Figure 5A). As mentioned before, use-dependent decreases in either firing rate or synaptic transmission over time can provide momentum to neuronal sequences by favoring the recruitment of new neurons that have not yet been activated during a network population event (Itskov et al., 2011;Romani and Tsodyks, 2015;Treves, 2004). To test 605 a possible role for cellular adaptation in sequence generation in our model network, we replaced the "bursting" excitatory cell model with a "regular spiking" model without spike rate adaptation ( Figure 5A). This cell model did not support the high instantaneous firing rates of the bursting cell model, which compromised the peak firing rates of excitatory cells in the network and their entrainment by higher frequency  Figure   S8E). These data show that adaptation of neuronal spiking provides a cellular-level 620  Figure 3D: E, p<0.00001; I, p=0.99997). (F) Offline sequence velocity (two-sided K-S tests: E vs. FF, p=0.15282; I vs. FF, p=0.73515; two-sided K-S tests vs. data from model with structured E ← E weights in Figure 3E: E, p=0.03099; I, p=0.91017). (G) Fraction of events that met criterion for sequences consistent with continuous spatial trajectories (see Materials and Methods) (one-sided paired t-tests: E vs. FF, p=0.00038, I vs. FF, p=0.00119; two-sided t-tests vs. data from model with structured E ← E 635 weights in Figure 3F: E, p<0.00001; I, p=0.03860). (H) Offline high-frequency rhythmicity (one-sided paired t-tests: 75 Hz -300 Hz frequency band: E vs. FF, p=0.00114; I vs. FF, p=0.00003; two-sided ttests vs. data from model with structured E ← E weights in Figure 3G: E, p<0.00001; I, p=0.00027). In (D -F) and (H), mean (solid) ± SEM (shading) were computed across 5 independent instances of each network. In (G), box and whisker plots depict data from 5 independent instances of each network model 640 (see Materials and Methods). p-values reflect FDR correction for multiple comparisons. mechanism for flexible and reversible sequence generation in recurrent spiking networks.

Preferential replay of reward location 645
Thus far, we have simulated network activity during spatial navigation, and identified features of the network that enable offline replay of behavioral sequences stored in memory. However, in these simulations all spatial positions were visited with equal occupancy, and considered to be of equal salience or relevance to the virtual animal.
This resulted in all positions being replayed with equal probability offline ( Figure 3C), 650 mimicking experimental conditions where all spatial positions contain discriminable sensory cues, and opportunities for reward are provided at random times and positions (Turi et al., 2019;Zaremba et al., 2017). However, it has been shown that when reward is provided at a fixed location that the animal must discover through learning, offline memory replay events become biased towards sequences of place cells that encode 655 positions nearby and including the site of reward (Gillespie et al., 2021;Ólafsdóttir et al., 2018;Pfeiffer, 2020;Singer and Frank, 2009). In parallel with the development of this bias in offline memory replay during learning, it has been shown that the fraction of hippocampal pyramidal cells that selectively fire along the path to reward increases (Lee et al., 2006;Turi et al., 2019;Zaremba et al., 2017). Here we sought to test the 660 hypothesis that this network-level over-representation of reward location is sufficient to bias the content of offline memory replay.
We chose a position along the virtual track to be the fixed location of a simulated reward, and biased the allocation of place field locations such that an increased fraction of excitatory neurons were selectively activated at positions near the 665 reward ( Figure 6A). As before, feedforward and recurrent synaptic connection strengths were increased between neurons with overlapping selectivity (Supplementary Figure   S9A). Aside from an enhanced fraction of active excitatory neurons near the reward site (Supplementary Figure S9B), this produced network dynamics during simulated navigation that conformed to experimental constraints for sparsity, selectivity, and 670 rhythmicity ( Figure 6A and Supplementary Figures S9B -S9E). During simulated offline

680
S tests: E vs. FF, p<0.00001; I vs. FF, p=0.41580; two-sided K-S tests vs. data from model with structured E ← E weights in Figure 3C: E, p<0.00001; I, p=0.84875). (E) Offline sequence length (twosided K-S tests: E vs. FF, p<0.00001; I vs. FF, p<0.00001; two-sided K-S tests vs. data from model with structured E ← E weights in Figure 3D: E, p=0.25150; I, p<0.00001). (F) Offline sequence velocity (twosided K-S tests: E vs. FF, p=0.02742; I vs. FF, p=0.73515; two-sided K-S tests vs. data from model with 685 structured E ← E weights in Figure 3E: E, p=0.16847; I, p=0.62023). (G) Fraction of events that met criterion for sequences consistent with continuous spatial trajectories (see Materials and Methods) (onesided paired t-tests: E vs. FF, p<0.00001, I vs. FF, p=0.00027; two-sided t-tests vs. data from model with structured E ← E weights in Figure 3F: E, p=0.05397; I, p=0.00754). (H) Offline high-frequency rhythmicity (one-sided paired t-tests: 75 Hz -300 Hz frequency band: E vs. FF, p=0.00067; I vs. FF, 690 p=0.00006; two-sided t-tests vs. data from model with structured E ← E weights in Figure 3G: E, p=0.00003; I, p=0.00001). In (D -F) and (H), mean (solid) ± SEM (shading) were computed across 5 independent instances of each network. In (G), box and whisker plots depict data from 5 independent instances of each network model (see Materials and Methods). p-values reflect FDR correction for multiple comparisons.

695
rest, the excitatory neurons in the network responded to random feedforward inputs by generating neuronal sequences corresponding to forward, reverse, and mixed direction trajectories through the environment (Figures 6B -6G), paced by high frequency oscillations in the inhibitory cells ( Figure 6H), as before ( Figures 3A -3G). However, 700 now positions near the simulated reward site were replayed in a higher proportion of replay events ( Figure 6D). This preferential replay of locations over-represented by the network recapitulated experimental findings and supported the hypothesis that nonuniform place cell allocation and biased memory replay are causally linked (Levy, 1989). 705

Discussion
In this study we used a simple recurrent spiking network model of hippocampal area CA3 to investigate the structural and functional requirements for offline replay of spatial memories. We optimized synaptic, cellular, and network parameters of the network to 710 produce population dynamics that match experimentally observed sparsity, selectivity and rhythmicity. We found that networks that fit these constraints exhibit additional emergent properties, including the ability to generate fast timescale memory-related neuronal sequences. During simulated spatial navigation, when ordered sensory information was provided on the seconds-long timescale of locomotion behavior, the 715 network produced neuronal sequences that swept from past to future positions on the faster timescale (~125 ms) of the theta rhythm ("theta sequences"). During simulated offline rest, the network responded to transient noisy activation of random, sparse inputs by generating neuronal sequences that corresponded to forward, reverse, or mixed direction trajectories through the spatial environment. 720 Both online and offline sequence generation depended on structure in the strengths of excitatory synaptic connections such that pairs of neurons with overlapping spatial tuning were more strongly connected. In the online phase, different sparse subsets of excitatory neurons were activated at different spatial positions due to structure in the strengths of connections from spatially-tuned feedforward afferent 725 inputs. The constraint that recurrent excitation must drive rhythmic synchronization in the theta band resulted in relatively strong recurrent connections. During each cycle of the theta rhythm, when the firing rates of the excitatory neurons were at their maximum, synaptic excitation from recurrent connections exceeded that from the feedforward afferents (Supplementary Figure S1E). This caused the sum of forward-730 moving feedforward inputs and symmetric mixed-direction feedback inputs to favor activation of cells encoding positions at or ahead of the current position. This generated forward-sweeping sequences that outpaced the speed of locomotion.
However, at the opposite phase of the theta rhythm, when the firing rates of the excitatory cells reached their minimum, the non-rhythmic feedforward input became 735 greater than recurrent excitation (Supplementary Figure S1E), causing theta sequences to reverse direction and relax back towards the current position encoded by the feedforward inputs.
In the offline phase, the feedforward inputs were not activated in a sequence, so momentum had to be entirely internally generated by the network. In this case, the 740 particular subset of active feedforward inputs initially selected a sparse subset of excitatory neurons to begin to fire, which set a starting position for the replayed trajectory. Slight biases in the feedforward input could then influence whether the active ensemble of excitatory neurons next recruited neurons encoding spatial positions in either the forward or reverse direction. Once activity began moving in one 745 direction, spike-rate adaptation facilitated continued sequence movement along that direction. However, depending on fluctuations in the feedforward inputs, sequences were also generated that included changes in direction. Interestingly, this process is akin to interpolation or smoothing -the recurrent connections within the network served to bridge large, discontinuous jumps in position encoded by the noisy 750 feedforward inputs with smaller, more continuous steps. This produced offline sequences that were consistent with the topology of the spatial environment, but did not necessarily replay exact experienced trajectories. These findings are consistent with a recent report that neuronal sequences activated during hippocampal SWRs in vivo resembled Brownian motion, or a random walk through the sensory space, rather 755 than precise replay of experience (Stella et al., 2019). This suggests that, rather than serving mainly to consolidate specific episodic memories of ordered sensory experiences, neuronal sequences during SWRs could also explore possible associations within the environment that had not been fully sampled during experience.
Our modeling results showing that increased population representations of goal sites 760 bias the content of offline memory replay also corroborate recent findings that previously rewarded locations are replayed more readily than immediate past or immediate future trajectories (Gillespie et al., 2021). Within this framework, synaptic plasticity during offline replay could modify connection strengths to increase the chance that a new path will be taken that is likely to lead to a desired outcome 765 (Ólafsdóttir et al., 2015).
In summary, our modeling results identified a minimal set of elements sufficient to enable flexible and bidirectional memory replay in neuronal networks: spike rate adaptation, and recurrent connectivity between excitatory and inhibitory neuron populations with strengths and kinetics optimized for rhythmogenesis and sparse and 770 selective stimulus representations. In previous models of neuronal sequence generation, additional network components were proposed to enable unidirectional sequences stored in memory to be reversed during offline recall, including neuromodulation (Gauy et al., 2018), excitability of neuronal dendrites (Gauy et al., 2018;Jahnke et al., 2015), coordinated plasticity at both excitatory and inhibitory 775 synapses (Ramirez-Villegas et al., 2018), and functional specialization of diverse subpopulations of inhibitory interneurons (Cutsuridis and Hasselmo, 2011). While these mechanisms may regulate and enhance memory replay, our results suggest that they are not necessarily required.
This model also makes some experimentally-testable predictions. First, it implies 780 that ion channel mutations that disrupt neuronal spike rate adaptation may also degrade neuronal sequence generation and memory consolidation (Peters et al., 2005).
Secondly, while the direction and content of offline sequences may be largely controlled by internal dynamics and information stored in the synaptic weights within a recurrent neuronal circuit, the model network still required a small amount of random 785 feedforward afferent input to evoke an offline population burst, suggesting that experimental manipulations of afferent projections to hippocampal area CA3 may alter the frequency or content of memory replay events (Chenani et al., 2019;Sasaki et al., 2018). Recent work has also begun to explore the advantages of generative replay for learning in artificial neural networks (Roscow et al., 2021). In addition to better 790 understanding the biological mechanisms of memory consolidation and flexible planning of behavior, characterizing the minimal mechanisms of memory replay could facilitate the engineering of artificial systems that can refine their internal representations of the environment during periods of offline rest (Buzsáki, 1989).

Acknowledgements
We are grateful to Kristopher Bouchard at LBNL for sharing large-scale computing resources provided by the National Energy Research Scientific Computing Center, a Department of Energy Office of Science User Facility (DE-AC02-05CH11231). This work was also made possible by computing allotments from NSF (XSEDE Comet, 800 NCSA Blue Waters, and TACC Frontera) and supported by NIH BRAIN U19 grant NS104590 and NIMH grant R01MH121979.

Materials and Methods
Simulations of a recurrent network of excitatory and inhibitory spiking neurons were 805 executed using the python interface for the NEURON simulation software (Hines et al., 2009). Cell models were single-compartment integrate-and-fire neuronal cell models, as defined by Izhikevich (Izhikevich, 2007), and as implemented for the NEURON simulator by Lytton et al (Lytton et al., 2016). Previously calibrated cell models were replicated from those previous reports without modification -the "intrinsically bursting 810 cell" model was used for excitatory neurons (E) with spike rate adaptation, the "regular spiking pyramidal cell" model was used for excitatory neurons without spike rate adaptation, and the "fast-spiking interneuron" model was used for inhibitory neurons (I) (Izhikevich, 2007;Lytton et al., 2016). Individual spikes in presynaptic neurons activated saturable conductance-based synapses with exponential rise and decay 815 kinetics after a constant delay of 1 ms to emulate axonal conduction time (Carnevale and Hines, 2006). Excitatory synapses had a reversal potential of 0 mV (like AMPA-type glutamate receptors), and inhibitory synapses had a reversal potential of -80 mV (like GABA(A)-type receptors). In addition to the excitatory (E) and inhibitory (I) neuron populations, a population of feedforward afferent inputs (FF) provided a source of 820 external excitatory synaptic drive to the model network.
The baseline weights of excitatory synapses onto E cells were sampled from a log-normal distribution (Almeida et al., 2009b;Buzsaki and Mizuseki, 2014), while the weights of excitatory synapses onto I cells, and all inhibitory synapses were sampled from a normal distribution . In addition to the random baseline 825 synaptic weights assigned to excitatory synapses onto E cells, input strengths were increased by a variable additive factor that depended on the distance between the place fields of cells with overlapping spatial selectivity (Supplementary Figure S1B).
The place field locations of the FF and E populations were assigned by distributing locations throughout the circular simulated track at equal intervals, and randomly 830 assigning them to cells within each population. Random connectivity resulted in each E neuron receiving inputs from many FF and E neurons with heterogeneous selectivity, which produced substantial out-of-field excitation at all positions along the track.
For each of six types of connections between the three cell types (E <-FF, E <-E, E <-I, I <-FF, I <-E, I <-I), a number of parameters were varied and explored during 835 optimization to identify model configurations that produced dynamics comparable to experimental observations. These parameters included: the mean and variance of the synaptic weight distribution for each connection type, the decay time constants of the synaptic conductances, the mean number of synapses made by one presynaptic cell onto one postsynaptic cell for each pair of cell types, and the maximum increase in 840 synaptic weight due to shared selectivity, as mentioned above. Self-connections were not permitted.
Optimization was performed using a population-based iterative multi-objective algorithm. During each of 50 iterations, a population of 600 models with different parameters were simulated for one trial of simulated online run, and for 5 trials of 845 simulated offline rest. During offline rest trials, a random subset of 25% of feedforward inputs were active with a mean rate of 12.5 Hz for an event duration of 160 ms (8 bins of 20 ms each). Different trials were implemented by using a distinct random number stream to sample unique spike times of the feedforward inputs from an inhomogeneous Poisson process. The following features of the network dynamics were 850 evaluated for each model: average minimum and maximum firing rates of E cells during run, average mean firing rates of I cells during run, average fraction of active E and I cells during run, mean firing rates of E cells during rest, average fraction of active E cells during rest, and finally, features related to the frequency and power of theta and gamma band oscillations in E and I cells during run. These features were compared to 855 target values to obtain a set of multiple objective error values. Models within a population were compared to each other and ranked with a non-dominated sorting procedure (Deb, 2011). Then, a new population of models was generated by making small perturbations to the parameter values of the most highly-ranked models from the previous iteration. This algorithm effectively identified model configurations that 860 satisfied multiple objective criterion. Below, the final optimized parameter values (Table   1) and measured features of the network dynamics (Table 2) are compared for various model configurations discussed in this study:

865
In the above table, theta and gamma amplitudes were quantified as follows: average population firing rates were band-pass filtered, and the envelopes of the filtered traces were computed from the Hilbert transformation. Then power was expressed as a ratio of the average envelope amplitude to the average population firing rate. To quantify theta and gamma frequency, bandpass filtered traces were subject to 870 frequency decomposition, and the frequency corresponding to the centroid or centerof-mass of the power spectral density distribution was taken as the dominant frequency within the band. The area of the power spectral density distribution was also used to compute a "frequency tuning index" which quantified how concentrated the power distribution was around the centroid frequency. This metric was akin to a signal-875 to-noise ratio in the frequency domain instead of the time domain, and was computed as follows: where is the average power at frequencies within the center of mass quartile containing the centroid frequency (signal), is the average power at frequencies in the 880 extreme high and low quartiles outside the center of mass quartile (noise), is the standard deviation of the power distribution, and is the half-width of the power distribution in the frequency domain normalized to the width of the bandpass filter. This metric has values near zero when power is distributed uniformly within the filter band, and values larger than one when power is concentrated around the centroid frequency. 885 Following parameter optimization, each variant of the network was evaluated by simulating 5 trials of online run, and 1000 trials of offline rest for each of 5 independent network instances. For a given set of model parameters, independent instances of each network variant were constructed by using distinct random number streams to assign place field locations, input spike times, synaptic connections, and synaptic 890 weights for all cells in the network.
Bayesian decoding of spatial position from spike times recorded during a single trial (Figures 2, 3, 5 and 6, and Supplementary Figures S3 -S6) was performed using the procedure described in (Davidson et al., 2009). The spatial firing rates of all cells were averaged across 5 trials of simulated online run to compute the spiking 895 probabilities of each neuron in 20 ms bins. Then, spiking data was taken from either a held-out set of 5 trials of simulated run (Figure 2 In Figure 4 and Supplementary Figure S8, the diversity and degeneracy of various model configurations was explored as follows: for each model configuration, 30,000 models were evaluated during parameter optimization, and the model with the 920 lowest multi-objective error score was considered the "best" model. The remaining models were sorted by their Euclidean distance from the "best" model in the space of model parameters. This resulted in an error landscape (e.g. Figure 4A) in which models with similar parameters resulted in similar multi-objective error scores. We then identified models located at local minima in this error landscape, which as a group 925 comprised models that were distant from each other in parameter space, but similar to each other in terms of overall multi-objective error. We further enforced that selected models had to be a minimum distance of 0.15 from each other in parameter space, and selected 5 such models with the lowest error score to be included in a "Marder group" of models for further analysis. For each alternative network model configuration (i.e. 930 network models with and without structured recurrent excitatory connections), each of 5 "Marder group" model variants with different parameters were evaluated for offline sequence generation by simulating 1000 trials for each of 5 independent network instances.

Data and Code Availability 940
All code necessary to reproduce the data and analysis presented in this work are available here: Network simulation and analysis code: https://github.com/neurosutras/optimize_simple_network Network optimization code: https://github.com/neurosutras/nested 945 than expected in both E and I cells due to suppression of background activity by synaptic inhibition (one-sided paired t-tests: Expected vs. Actual: E, p<0.00001; I, p<0.00001). (D) The distance between the expected location of maximum firing, and the actual location is quantified as a fraction of the circular track (one-sample t-tests: E, p<0.00001; I, p=0.00005). (E) Traces depict average population firing rates of E and FF cells during an example trial for one instance of the network. E cells and FF cells dominate at different phases of the population theta oscillation. (F) In Figure 1G (bottom row), the amplitudes of the gamma-filtered (purple) population firing rates for E and I populations vary in time. Here, the amplitudes or envelopes of the gamma-filtered population firing rates were subject to frequency decomposition. The resulting frequency distributions show peaks in the theta band (one-sided paired t-tests: 4 Hz -10 Hz frequency band: E vs. FF, p=0.00002; I vs. FF, p=0.00002). Mean (solid) ± SEM (shading) were computed across 5 independent network instances. In (C) and (D), box and whisker plots depict data across cells for one example instance of the network (see Materials and Methods). Statistics were computed across 5 independent instances of the network. p-values reflect FDR correction for multiple comparisons.  Figure S1B for an alternative network model with random E ← E weights. (B -F) Same as Figures 1C -1F and 1H for alternative network model. (F) Rhythmicity (one-sided paired t-tests: theta: E vs. FF, p=0.00003; I vs. FF, p<0.00001; gamma: E vs. FF, p<0.00001; I vs. FF, p<0.00001; two-sided ttests vs. data from model with structured E ← E weights in Figure 1H: theta: E, p<0.00001; I, p<0.00001; gamma: E, p<0.00001; I, p<0.00001). (G -L) Same as (A -F) for an alternative network model with structured, but shuffled E ← E weights. (H) Rhythmicity (one-sided paired t-tests: theta: E vs. FF, p=0.00009; I vs. FF, p=0.00003; gamma: E vs. FF, p=0.00002; I vs. FF, p<0.00001; two-sided t-tests vs. data from model with structured E ← E weights in Figure 1H: theta: E, p<0.00001; I, p<0.00001; gamma: E, p=0.00001; I, p<0.00001). In (C -F) and (I -L), data were first averaged across 5 trials per network instance. Mean (solid) ± SEM (shading) were computed across 5 independent network instances. pvalues reflect FDR correction for multiple comparisons.  Figure 3C: E, p=0.99974; I, p=0.99974). (F) Offline sequence length (two-sided K-S tests: E vs. FF, p<0.00001; I vs. FF, p<0.00001; two-sided K-S tests vs. data from model with structured E ← E weights in Figure 3D: E, p<0.00001; I, p=0.31278). (G) Offline sequence velocity (two-sided K-S tests: E vs. FF, p<0.00001; I vs. FF, p=0.00026; two-sided K-S tests vs. data from model with structured E ← E weights in Figure 3E: E, p<0.00001; I, p=0.00026). (H) Fraction of events that met criterion for sequences consistent with continuous spatial trajectories (see Materials and Methods) (one-sided paired t-tests: E vs. FF, p<0.00001, I vs. FF, p=0.00003; two-sided ttests vs. data from model with structured E ← E weights in Figure 3F: E, p<0.00001; I, p=0.00112). (I) Offline high-frequency rhythmicity (one-sided paired t-tests: 75 Hz -300 Hz frequency band: E vs. FF, p=0.00100; I vs. FF, p=0.00001; two-sided t-tests vs. data from model with structured E ← E weights in Figure 3G: E, p=0.00001; I, p=0.02509). In (E -G) and (I), mean (solid) ± SEM (shading) were computed across 5 independent instances of each network. In (H), box and whisker plots depict data from 5 independent instances of each network model (see Materials and Methods). p-values reflect FDR correction for multiple comparisons.  Figures 3A -3G for an alternative network model with structured, but shuffled E ← E weights. (A) and (B) correspond to two example trials from one example network instance. (C) Decoded positions (two-sided K-S tests: E vs. FF, p=0.99974; I vs. FF, p=0.99974; two-sided K-S tests vs. data from model with structured E ← E weights in Figure 3C: E, p=0.99974; I, p=0.99974). (D) Offline sequence length (two-sided K-S tests: E vs. FF, p<0.00001; I vs. FF, p=0.11015; two-sided K-S tests vs. data from model with structured E ← E weights in Figure 3D: E, p<0.00001; I, p<0.00001). (E) Offline sequence velocity (two-sided K-S tests: E vs. FF, p=0.96717; I vs. FF, p=0.73515; two-sided K-S tests vs. data from model with structured E ← E weights in Figure 3E: E, p=0.00272; I, p=0.24955). (F) Fraction of events that met criterion for sequences consistent with continuous spatial trajectories (see Materials and Methods) (one-sided paired t-tests: E vs. FF, p=0.06713, I vs. FF, p=0.99967; two-sided t-tests vs. data from model with structured E ← E weights in Figure 3F: E, p<0.00001; I, p<0.00001). (G) Offline high-frequency rhythmicity (one-sided paired t-tests: 75 Hz -300 Hz frequency band: E vs. FF, p=0.00092; I vs. FF, p=0.00004; two-sided t-tests vs. data from model with structured E ← E weights in Figure 3G: E, p=0.00006; I, p=0.00059). In (C -E) and (G), mean (solid) ± SEM (shading) were computed across 5 independent instances of each network. In (F), box and whisker plots depict data from 5 independent instances of each network model (see Materials and Methods). p-values reflect FDR correction for multiple comparisons.