The essential role of hippocampo-cortical connections in temporal coordination of spindles and ripples

The predominant activity of slow wave sleep is cortical slow oscillations (SOs), thalamic spindles and hippocampal sharp wave ripples. While the precise temporal nesting of these rhythms was shown to be essential for memory consolidation, the coordination mechanism is poorly understood. Here we develop a minimal hippocampo-cortico-thalamic network that can explain the mechanism underlying the SO-spindle-ripple coupling indicating of the succession of regional neuronal interactions. Further we verify the model predictions experimentally in naturally sleeping rodents showing our simple model provides a quantitative match to several experimental observations including the nesting of ripples in the spindle troughs and larger duration but lower amplitude of the ripples co-occurring with spindles or SOs compared to the isolated ripples. The model also predicts that the coupling of ripples to SOs and spindles monotonically enhances by increasing the strength of hippocampo-cortical connections while it is stronger at intermediate values of the cortico-hippocampal projections.


Introduction
During slow wave sleep (SWS), the deepest stage of non-rapid eye movement (NREM) sleep, cortical network generates synchronous transitions between a quiescent down and an active up state with high frequency activity in beta and gamma frequency ranges being nested at the up state ( Cox et al., 2014 ;Le Van Quyen et al., 2010Valderrama et al., 2012 ). These oscillations coexist with sleep fast spindles which are waxing and waning oscillations generated in the thalamus and transmitted to the cortex via thalamocortical projections ( Timofeev and Bazhenov, 2005 ). Sharp wave ripples (SWRs) are another spontaneously generated NREM rhythms that are generated in hippocampus. SWR consists of high frequency oscillations (100-250 Hz) generated in CA1 network due to strongly enhanced excitatory CA3 inputs characterized by slow potential deflections ( Buzsaki, 2015 ). Previous findings confirmed the coupling among these different NREM oscillatory rhythms, especially the nesting of SWRs in the spindles troughs ( Binder et al., 2019 ;Staresina et al., 2015 ) and occurrence of fast spindles at the onset of the up state ( De Gennaro and Ferrara, 2003 ;Klinzing et al., 2016 ;Mölle et al., 2011 ). The coupling between field potentials of hippocampal sharp-wave ripples, thalamo-cortical sleep spindles and slow oscillations is due to long-range connections among these different networks which is essential for bidirectional interaction between thalamus and cortex as well as between hippocampus and cor-neously generated due to local short-range interactions between excitatory and inhibitory neurons of adapting cortical and hippocampal networks and bursting thalamic network respectively. As a result of nonlinear interactions of multiple oscillators, this simple deterministic model shows irregular behavior with variability from cycle to cycle and generates spindles and ripples with variable durations and amplitudes. The SO-spindle-ripple coupling emerges as a result of interplay between thalamus and cortex as well as cortex and hippocampus controlled by longrange connections between these networks. Further to verify the predictions of the model experimentally, we use simultaneous prefrontal cortical (PFC) and hippocampal LFP recordings in drug-free rodents during SWS and find consistent results with the model results in a quantitatively precise fashion.

CA3-CA1 network
The hippocampal network consists of CA3 and CA1 network ( Fig. 2 A). Similar to the models developed in our previous papers ( Ghorbani et al., 2012 ;Hashemi et al., 2019 ), the point neural firing rate model of CA3 (CA1) network consists of one group of identical excitatory neurons with the membrane potential shown by 3 ( 1 ) , and one group of identical inhibitory neurons with the membrane potential shown by V CA3 Here, r 1 ∼∼70 Hz is the maximal firing rate in a completely depolarized state while V * k m is the threshold firing potential and g k m shows the sharpness of the firing rate dependence on the membrane potential.
All the possible recurrent connections within each network except the excitatory to excitatory connection in CA1 and inhibitory to inhibitory connection in CA3 are considered ( Sik et al., 1993 ;Taxidis et al., 2012Taxidis et al., , 2013Aussel et al., 2018 ). The excitatory neurons of CA3 are connected to both excitatory and inhibitory neurons of CA1. The probability of the connection from neuron m to neuron n is given by P k mn where k equals to CA3 or CA1 shows local connections within each network and k = CA3-CA1 corresponds to long-range connections from CA3 to CA1. m and n can be either e or i corresponding to excitatory and inhibitory neurons respectively (Table S1, P CA3 ee = 0.30, P CA3 ei = 0.3, P CA3 ie = 0.50, P CA3−CA1 ei = 0.2, P CA1 ei = 0.3, P CA1 ii = 0.25 are considered similar to a previous study ( Ramirez-Villegas et al., 2018 )). J k mn shows the strength of the connection from neuron m to neuron n . Dendritic spike frequency adaptation (DSFA) is considered for all the excitatory to excitatory connections so that J k ee is a decreasing sigmoid function of an adaptation variable, c representing the mean adaptation level of an excitatory neuron Here, J0 k ee is the maximal synaptic strength corresponding to adaptation level at equilibrium value (set to be zero). * = 10 is the threshold adaptation level and g c = 3 shows the sharpness of the excitatory to excitatory connection dependence on the adaptation variable. When the excitatory neuron receives input from other excitatory neurons, the adaptation level increases by Δc n , and it slowly decays to its equilibrium value in the absence of excitatory inputs. τ c = 500 ms is the adaptation time constant which is much larger than the membrane time constant ( τ e = 20 ms and τ i = 10 ms for excitatory and inhibitory neurons respectively). Hence the CA3-CA1 network is described by 6 rate equa-tions:

Thalamocortical model
The thalamocortical model has been described previously ( Hashemi et al., 2019 ). Briefly, the thalamocortical network consists of two weakly coupled cortical networks ( N CX e = 8000, N CX i = 2000) and one thalamic network consisting of one group ( N TH t = 500) of identical thalamocortical (TC) neurons, V TH t and one group ( N TH r = 500) inhibitory neurons representing thalamic reticular neurons (RE), V TH r ( Hashemi et al., 2019 ). Similar to CA3-CA1 excitatory neurons, the adaptation is considered for cortical excitatory neurons which enables the isolated cortical network to generate up and down states. To model T-type calcium currents which is essential for spindle generation, a bursting variable, u m is defined for thalamic neurons so that they can show both tonic and burst modes. The rate equations for thalamic neurons are given by: If V TH r > 0 mV then b r = 0 otherwise b r = − 200 mA If V TH t > − 0.1 mV then b t = 0 otherwise b t = − 200 mA Where A m is the specific membrane capacitance. When the cell is hyperpolarized, u m starts to decrease towards its new equilibrium value ( b m ) with a time constant τ u m which is much larger than the membrane time constant. The firing rate of the thalamic neurons can show tonic and burst modes depending on L m u m , where L m is a positive parameter that controls the sharpness of switch between the two modes.
show the maximum firing rate, threshold firing potential and the sharpness of the firing rate dependence on the membrane potential during tonic (burst) mode.
The cortical and thalamic networks are connected by long-range connections from the excitatory neurons of one network to both excitatory and inhibitory neurons of the other network. The parameters of the thalamocortical network are almost set to be the same as our recent paper ( Hashemi et al., 2019 ; also see Table S2) so that the thalamocortical model can spontaneously generate up and down states as well as spindles.

Hippocampo-cortico-thalamic model
The existence of a two-way communication between the PFC (anterior cingulate cortex) and hippocampus has been reported previously ( Rajasethupathy et al., 2015 ). Hence similar to a recent modeling study ( Sanda et al., 2021 ) in the full model consisting of CA3-CA1 networks and the thalamocortical networks, the cortical networks are recurrently connected to CA3-CA1 network so that all the anatomical connections from the thalamocortical network to the hippocampus are modeled by one projection from the cortex to the hippocampus ( Fig. 3 A). Our model is a minimal neural mass model of the networks which are essential for the generation of the 3 important oscillatory activities during SWS i.e. SOs, spindles and ripples. Modeling the entorhinal cortex or midline thalamic nucleus reuniens which comprise additional relays between PFC and hippocampus ( Preston and Eichenbaum, 2013 ) would substantially go beyond the scope of our simple model. Hence in the full model, the thalamic network which simply consists of one group of excitatory neurons and one group of inhibitory neurons is only connected to cortical networks similar to Sanda et al., 2021 ( Fig. 3 A). The cortical, thalamic and hippocampal networks are connected by long-range connections from the excitatory neurons of one network to both excitatory and inhibitory neurons of the other network. The probability and strength of long-range connections from neuron m of network k to the neuron n of network h are given by P k−h mn and J k−h mn respectively. k and h can be CA1, CA3, CX (representing cortex) or TH (representing thalamus). n can take t (representing TC neuros), r (representing thalamic RE neurons), e (representing cortical or hippocampal excitatory neurons) and i (representing cortical or hippocampal inhibitory neurons) while m can be only either t or e which are excitatory neurons. The synaptic strengths and variables of cortical network 2 are indicated by primes. The full model is described by 14 rate equations.
The strength and probability values of the connections between the hippocampal and cortical networks are reported in Table S3. The nonlinear equations of the model were solved using the fourth-order Runge-Kutta method with a time step of 1 ms by MATLAB R2017b software.

Experimental data
In this research, two datasets have been used which were recorded in Professor Buzsaki's lab and taken from the public data sharing repository http://crcns.org. The original studies were approved by the local ethics committees (Institutional Animal Care and Use Committee of New York University Medical Center and Rutgers University). In the first dataset ( Fujisawa et al., 2008( Fujisawa et al., , 2015, an adult male (3-5 months old) rat was implanted with silicon probes in the prefrontal cortex and hippocampal wire bundles were implanted above CA1. The probes consisted of 8 or 4 shanks separated by 200 m and each shank had 8 recording sites (160 m 2 each site, 1-3 M Ω impedance). Number of days and recording sessions during sleep was 5 days, 33 sessions, each with 96 recording sites (64 recording sites in PFC, 32 channels in CA1). During the recording sessions, neurophysiological signals were acquired continuously at 20 kHz on DataMax system (16-bit resolution). The broadband signal was down sampled to 1.25 kHz and used as LFP. In the second dataset ( Peyrache et al., 2015a( Peyrache et al., , 2015b, for four male mice weighing ∼∼30 g (3-6 months), hippocampal wires were implanted above CA1, each probe consisted of 8 shanks separated by 200 m and each shank had eight recording sites (160 m 2 each site, 1-3 M Ω impedance). Number of days and recording sessions during sleep is 18 days, 18 sessions. Electrophysiological signals were acquired at 20 kHz on a 256-channel Amplipex system (16-bit resolution). The broadband signal was down sampled to 1.25 kHz and used as LFP. For this dataset, SWS was detected using CA1 LFP spectrogram and animal movement ( Viejo and Peyrache, 2020 ). For the first dataset, epochs of cortical LFP spectrogram (overlapping 10 s time windows) with delta (0.5-4 Hz) to theta (6-10 Hz) power ratio greater than 2.95 were considered as SWS. All the analyses were conducted during SWS.

Ripple, spindle and SO detection
For the model data, the spindles/SOs and SWRs were detected from the membrane potential of cortical and CA1 pyramidal neurons respectively. The algorithm for detecting ripples was similar to a previous study ( Levenstein et al., 2019 ). The hippocampal LFP (membrane potential of pyramidal neurons for the model data) from the CA1 was band pass filtered in the ripple frequency range (100-250 Hz, 4th order zero-phase delay Butterworth) and the RMS signal was smoothed using 50-ms Gaussian window. A ripple event was identified when the smoothed RMS signal was higher than a threshold (median + 4 SD for experimental and median + 2 SD for model data) for a minimum duration of 30 ms. The amplitude of the filtered signal in the frequency range of 5-20 Hz (similar to Ramirez-Villegas et al., 2018 ) at the time of the SWR maximum peak was defined as sharp wave amplitude for the experimental data. For the model data, the difference between the maximum and minimum of the membrane potential of CA3 excitatory neuron during the SWR was considered as the sharp wave amplitude.
The algorithm for detecting spindles was similar to the past paper ( Klinzing et al., 2016 ). In this manuscript we only focused on the fast spindles since there is no strong evidence for the existence of the slow spindles in rodents. Cortical LFP signal (membrane potential of pyramidal neurons for the model data) was first band pass filtered in the spindle frequency range (11-15 Hz for experimental data and 12-16 Hz for model data) using finite-impulse-response (FIR) filters from the EEGLAB toolbox ( Delorme and Makeig, 2004 , FIR band-pass filter, filter order corresponds to 3 cycles of the low frequency cut off). Next instantaneous amplitude was computed via the Hilbert transform and smoothened using a 300 ms Gaussian window. Putative spindles were selected if their amplitude was higher than a threshold (2.5 SD for experimental and 1 SD for model data) and their duration was between the lower and upper thresholds for duration (0.4-3.5 s). Two nearby events were merged if they were closer than 0.5 s.
The FMA toolbox ( Zugaro et al., 2018 ) was used to detect SOs in the experimental data. Method of detecting SOs in model data was similar to the method used by the previous paper ( Klinzing et al., 2016 ). First, the signal is filtered between 0.5 and 4 Hz (FIR band-pass filter, filter order corresponds to 3 cycles of the low frequency cut off). SOs were identified in the LFP when (i) the distance between consecutive positive-tonegative zero crossings was between 0.5 and 2 s (corresponding to 0.5 to 2 Hz), (ii) negative peaks greater than a threshold (2 SD for experimental data and 1.4 of the mean peaks for model data) were identified as SO negative peaks, and (iii) difference between the positive and negative peaks of the signal must be greater than a threshold (3.5 SD for experimental data and 1.4 times the average amplitude for model data).

Time-frequency analysis
The method to calculate the time-frequency representations (TFRs) and phase amplitude coupling was similar to the previous papers ( Staresina et al., 2015 ;Ladenbauer et al., 2017 ). TFRs were calculated for every event by mtmconvol function of the FieldTrip toolbox ( Oostenveld et al., 2011 ) with frequency steps of 0.25 Hz in the frequency range of 5-20 Hz for SO-spindle coupling and 5-300 Hz for spindle-ripple coupling. Sliding (10 ms steps) Hanning tapered window with a variable length including five cycles were used to ensure reliable power estimates. For experimental data, time-locked TFRs of all epochs around an event were then normalized as percentage change from the pre-event baseline ( − 2.5 to − 1.5 s for SOs and spindles, − 1.5 to − 1.0 s for ripples).

Phase amplitude coupling
To quantify synchronization and locking between the lowerfrequency (modulating) oscillatory event and the power fluctuation of the higher-frequency (modulated) oscillation (that is, TFR bins averaged across the respective frequencies), we extracted the phase values of these time series using the Hilbert transform and defined a synchronization index (SI). To ensure proper phase estimation, both lower-frequency and power of the higher-frequency time series were filtered in the range of lower-frequency (modulating) oscillatory event (SO: 0.5-2 Hz; spindle: 11-15 Hz for experimental and 12-16 Hz for model data; two-pass FIR band pass filter, order = 3 cycles of the low frequency cut-off). The SI was then calculated between the two phase value time series for each event epoch.
where n is the number of time points, φ ut represents the phase value of the fluctuations in the upper frequency power time series at time point t , and φ lt represents the phase value of the lower frequency band time series at time point t . The used interval for estimating the preferred phase was − 1 s to + 1 s around the SO peak and − 0.25 s to + 0.25 s around the spindle peak. The resulting SI is a complex number with the angle representing the phase shift between the two oscillatory events while the absolute value of SI represents the strength of the coupling between the two oscillatory events.

Phase-locking analyses
To find the SO phase at which the ripples/spindles preferentially occurred, the instantaneous phase of cortical excitatory neuron LFP filtered in the SO frequency range was estimated by Hilbert transform.

Statistical analysis
Statistical tests were performed with MATLAB and for circular analyses CircStat toolbox ( Berens, 2009 ) was used. Pearson correlation was used to compute the correlations. The comparison of the spindle and ripple properties between the isolated and coupled events were performed using Wilcoxon rank sum test. Rayleigh-test was used to test for non-uniformity of the phases. The codes will be available from the corresponding author on request.

Co-occurrence of SOs, spindles and ripples during SWS
We used two publically available datasets to investigate the coupling of SOs, spindles and ripples during SWS and validate the model predictions. The first dataset ( Fujisawa et al., 2008( Fujisawa et al., , 2015 which was used for all the analyses involving SOs/spindles consisted of 33 recording sessions (one rat), each with 64 recording sites in PFC (8 shanks) and 32 recording site (4 shanks) in CA1. The SOs/spindles and SWRs were detected from PFC and CA1 LFPs respectively. For all the analyses involving only SWRs, in addition to the first dataset, we used a second dataset ( Peyrache et al., 2015a( Peyrache et al., , 2015b consisting of 18 recording sessions (4 mice) each with 8 recording sites (1 shank) in CA1.
During SWS, in the absence of any stimuli, the cortical network generates synchronous slow oscillations between an active up and a quiescent down state ( Fig. 1 A). The duration of the up state was much larger than the duration of the down state which corresponded to large positive deflections of LFP with a width of few 100 ms when the MUA was very low. High frequency oscillations in beta and gamma frequency range were present during the up state of the SO indicating high synaptic activity ( Cox et al., 2014 ;Le Van Quyen et al., 2010Valderrama et al., 2012 ). Thalamocortical spindles which are waxing and waning oscillations at 11-15 Hz are another spontaneously generated oscillatory phenomenon that are present during SWS mostly at the active up state phase of SO. Spindles were detected in the cortical LFP as large amplitude oscillations in the spindle frequency range with a minimum duration of 0.4 s ( Fig. 1 A-C, mean ± SD: amplitude = 3.21 ± 0.94 SD, duration = 0.99 ± 0.51 s, n = 720). During SWS, the hippocampal LFP in CA1 showed large negative deflections with high frequency oscillations and large MUA activity corresponding to SWRs ( Fig. 1 A-C, mean ± SD: amplitude = 5.72 ± 1.9 SD, duration = 0.04 ± 0.02 s, n = 29,855). Similar to the previous findings ( Stark et al., 2014 ), ripple frequency was positively correlated with the sharp wave amplitude ( Fig. 1 D, r = 0.44, p = 3.23 ×1 0 −11 , n = 220 recording sites, Pearson correlation). This result indicates that the high frequency oscillations in CA1 are associated with the strong input from CA3.
Next we developed a hippocampo-cortico-thalamic model that can spontaneously generate SOs, spindles and SWRs and then quantified and compared the temporal relationship among these rhythms in the model and experimental data.

Generation of SWRs in a simple CA3-CA1 network
In our model, each of CA3 and CA1 networks was composed of one group of identical excitatory (E) neurons and one group of identical inhibitory (I) neurons Fig. 2 A). Each neuron was modeled by mean field dynamics to reduce the dimensionality and number of parameters. The two networks were connected by the long-range excitatory connections from the excitatory neurons of CA3 to both excitatory and inhibitory neurons of CA1. E-E, E-I, I-E recurrent connections with E-E connections subject to dendritic spike frequency adaptation (DSFA) were considered for the CA3 network. As we described previously ( Ghorbani et al., 2012 ), the recurrent connections in this simple E-I network with adaptation can spontaneously generate oscillations between active and quiescent states with a frequency controlled by the adaptation recovery time. The ratio of the duration of active to the duration of quiescent state was negatively correlated with the strength of the adaptation process, Δc CA3 , showing the level of increase of the adaptation variable of CA3 excitatory neuron when it receives input from other CA3 excitatory neuron ( Fig. 2 B, see Materials and methods for details). For large values of this parameter, the duration of the quiescent state was much larger than the duration of the active state and this active transient state with a duration of ∼∼ 100 ms represented sharp waves that were spontaneously generated in the CA3 network corresponding to large negative defections in the hippocampal LFP ( Fig. 2 C). The amplitude of the sharp wave was positively correlated with the magnitude of the synaptic connections from the excitatory to excitatory neurons in CA3 ( Fig. 2 D). CA1 network also was composed of identical groups of E and I neurons with only E-I, I-E and I-I connections ( Ramirez-Villegas et al., 2018 ;Stark et al., 2014 ). The large input from CA3 to CA1 neurons during the active state together with the recurrent E-I and I-E connections in the CA1 network, resulted in high frequency oscillations in the ripple frequency range in CA1 network ( Fig. 2 E). In this model, the dynamic of the synapses was considered to be very fast (synaptic time constant much lower than the membrane time constant) so that no rate equation was considered for the synapses (see Materials and methods for details). To investigate the underlying mechanism of the SWRs generation in the hippocampal network, we further solved analytically the rate equations for CA1 network subject to a depolarizing input current representing the excitatory input from CA3. For the range of hippocampal parameters for which the eigenvalues of the Jacobean matrix ( Eq. (9) -( (12) , Supplementary material) at the active fixed point were complex conjugates, the frequency of the oscillations around the active fixed point was given by Eq. (14) , Supplementary material. As was obtained both analytically ( Eq. (14) , Supplementary material) and with simulation ( Fig. 2 F), the high frequency oscillations in the frequency range of ripples can be obtained in this simple firing rate model by considering low values of the g CA1 e and g CA1 i corresponding to sharp dependence of the firing rate on the membrane potential of CA1 neurons ( Fig. 2 F). The frequency of the oscillations in CA1 also depended on the magnitude of the synaptic connections from excitatory to excitatory neurons in CA3 (Fig. S1). As a result, in agreement with the experimental observation ( Fig 1 D), the ripple frequency was positively correlated with the amplitude of the sharp waves corresponding to the magnitude of the excitatory input to CA1 network ( Fig. 2 G). Consistent with these results, the frequency of the oscillations around the active fixed point was found to be a monotonic function of the input current to CA1 network (Supplementary material, Eq. (14) ).
To investigate the effect of inhibition in CA1 on ripple generation we next computed the frequency of the oscillations in CA1 network versus the strength of the input from CA1 inhibitory neurons to CA1 excitatory neurons ( Fig. 2 H). Only for a limited range of this inhibitory input (around 0.8 to 1.8 times the default value of J CA1 ie reported in Table  S1), the CA1 network generated the oscillations in the ripple frequency range. For lower values, the oscillations in the CA1 network disappeared and for larger values, the CA1 network oscillated with frequencies lower than the ripple frequency range consistent with the analytical results ( Eq. (14) , Supplementary material). These results confirmed the importance of the inhibition in ripple generation. Moreover, the duration of the ripples was negatively correlated with the strength of the adaptation in CA1 network corresponding to Δc CA1 ( Fig. 2 I). Larger adaptation resulted in more decaying of the excitatory input received by CA1 excitatory neurons and hence lower duration of the ripples.
Taken together, we introduced a novel mechanism to generate SWRs in a simple mean field CA1-CA3 model subject to adaptation with instantaneous synaptic dynamics. The modeled SWRs, replicated the charac-   Table S1. E) Top, average TFR of CA1 excitatory neuron locked to ripple trough. Bottom, grand average ( ± SEM) ripple in the filtered signal from model data time-locked to ripple trough ( n = 305). F) The oscillation frequency in the CA1 network during the large CA3 inputs as a function of g CA1 e (Left) and g CA1 i (Right) representing the sharpness of the firing rate dependence on the membrane potential for CA1 excitatory and inhibitory neurons respectively. G) Relationship between sharp wave amplitude and ripple frequency. Each point corresponds to a different value of J CA3 ee . H) The ripple frequency versus the strength of the input from CA1 inhibitory neurons to CA1 excitatory neurons, α CA1 ie , which shows the factor multiplied to the default value of J CA1 ie reported in Table S1. The solid curve shows the quadratic fit to the simulation data points. I) Average duration of ripples decreases with increasing Δc CA1 . The solid curve shows the quadratic fit to the simulation data points.
teristics of SWRs during SWS especially the positive correlation between the sharp wave amplitude and ripple frequency. We further solved the rate equations analytically and found an estimate of the frequency of the oscillations in the CA1 network during the large excitatory input showing a monotonic increase of the frequency of oscillations as a function of the input current.

Spontaneously generated rhythms in the full model
To investigate the mechanism underlying the SO-spindle-ripple coupling, we developed a hippocampo-cortico-thalamic model. For the thalamocortical networks we used a simple mean field model consisting of two cortical and one thalamic networks that we previously introduced ( Hashemi et al., 2019 ). This model can spontaneously generate up and down states and spindles as a result of the recurrent connections between the neurons together with the adaptation and bursting properties in the cortical and thalamic networks respectively. As was discussed extensively in our previous paper ( Hashemi et al., 2019 ), the isolated cortical network in our model can produce broad frequency activity in high frequency (beta and gamma) range during the up state consistent with previous experimental reports in cortical slices ( Compte et al., 2008 ). In the thalamocortical network, fast spindles were generated in the thalamus due to bursting properties of the thalamic neurons so that the beginning and termination of them were controlled by the long-range cortico-thalamic input ( Hashemi et al., 2019 ).
In the full hippocampo-cortico-thalamic model, the CA3-CA1 and cortical neurons were connected by excitatory long-range recurrent connections that were much weaker than the short-range synaptic connections within each network resembling a small world network ( Fig. 3 A). In the absence of any input, the full model can generate SOs and spindles in the thalamocortical networks and SWRs in the CA1 network. Similar to the experimental results during SWS, the up state duration was much larger than the down state duration and high frequency oscillations were nested at the up state ( Fig. 3 B). Moreover, while the model was deterministic with no intrinsic or external noise, it generated irregular oscillations with variable duration and amplitude of spindles ( Fig. 3 C, mean ± SD: amplitude = 2.84 ± 0.54 SD and duration = 0.69 ± 0.17 s, n = 713) and SWRs ( Fig. 3 D, mean ± SD: amplitude = 7.30 ± 0.72 SD and duration = 0.0427 ± 0.004 s, n = 8693) as well as irregular up and down states characterized by deterministic noise ( Fig. 3 B) due to nonlinear interactions of dephased weakly coupled oscillators with different frequency of oscillations.
Overall, for the first time to our knowledge, we developed a minimal deterministic hippocampo-cortico-thalamic network that can spontaneously generate SOs, spindles and SWRs with irregular characteristics resembling the NREM sleep rhythms. We next used this simple model to compute and explain the precise temporal coordination of these rhythms which is expected due to long-range excitatory connections among different networks.

The coupling of thalamic spindles and hippocampal ripples
We next investigated the temporal relationship between thalamocortical spindles and hippocampal ripples in both experimental and model data. During SWS, the majority of the ripples co-occurred with the spindles ( Fig. 4 A, for 79.2% of ripples, spindles occurred within a ± 0.5 s interval around their maximum peak). Next we compared the properties of the isolated ripples with ripples co-occurring with spindles. The duration of ripples co-occurring with spindles was larger than the duration of the ripples occurring alone ( Fig. 4 B, p = 3.38 ×1 0 −7 , Wilcoxon rank sum, mean ± SD: duration = 0.04 ± 0.01 s, n = 957 for ripples coupled with spindles and duration = 0.03 ± 0.01 s, n = 250 for isolated ripples) while the amplitude of ripples coupled to the spindles was lower compared to the ripples occurring alone ( Fig. 4 B, p = 9.15 ×1 0 −6 , Wilcoxon rank sum, mean ± SD: amplitude = 6.9 ± 2.94 SD, n = 957 for ripples coupled with spindles and amplitude = 7.13 ± 2.87 SD, n = 250 for isolated ripples).
We next investigated the spindle-ripple coupling in the model data. In all the subsequent figures and discussion, the membrane potential of the excitatory neurons of the model data was multiplied by − 1 to enable easy comparison with the experimental LFP data. In a quantitative match to the experimental data, 72.5% of the ripples were co-occurred with spindles in the model data ( Fig. 4 C). Moreover, the duration of the ripples co-occurring with spindles was larger compared to the ripples occurring alone ( Fig. 4 D, p = 1.96 ×1 0 −7 , Wilcoxon rank sum, mean ± SD: duration = 0.04 ± 0.003 s, n = 6307 for ripples coupled with spindles and duration = 0.03 ± 0.002 s, n = 2386 for isolated ripples) while these isolated ripples had larger amplitude ( Fig. 4 D, p = 7.96 ×1 0 −3 , Wilcoxon rank sum, mean ± SD: amplitude = 6.67 ± 0.78 SD, n = 6307 for ripples coupled with spindles and amplitude = 7.41 ± 0.79 SD, n = 2386 for isolated ripples).
Further, during SWS, the majority of the spindles were coupled to the ripples (for 85.9% of spindles, ripples occurred within a ± 1 s interval around their maximum peak). We next compared the amplitude and duration of spindles that were coupled to ripples with the spindles uncoupled to ripples. Both duration and amplitudes of the spindles coupled to the ripples were larger compared to the spindles occurring alone. ( Fig. 4 E, p = 2.94 ×1 0 −6 and p = 4.06 ×1 0 −8 for duration and amplitude respectively, Wilcoxon rank sum, mean ± SD: duration = 1.03 ± 0.54 s and amplitude = 3.23 ± 0.9 SD, n = 619 for spindles coupled with ripples and duration = 0.81 ± 0.41 s and amplitude = 2.65 ± 0.75 SD, n = 101 for isolated spindles). Consistent with the experimental data, a large portion of the spindles in the model data were coupled to the ripples (98.6%). Comparison between spindles coupled and uncoupled with ripples revealed similar results to the experimental data so that the duration and amplitude of the spindles co-occurring with the ripples were larger than the isolated spindles ( Fig. 4 F, p = 3.83 ×1 0 −3 and p = 9.15 ×1 0 −3 for duration and amplitude respectively, Wilcoxon rank sum, mean ± SD: duration = 0.69 ± 0.15 s and amplitude = 8.95 ± 1.59 SD, n = 709 for spindles coupled with ripples and duration = 0.53 ± 0.09 s and amplitude = 7.39 ± 0.6 SD, n = 10 for isolated spindles).
To further explore the interaction between spindles and ripples, we employed phase amplitude coupling (PAC) analyses to quantify the modulation of hippocampal ripple amplitude with the phase of detected spindles in the cortical LFP. The PAC analysis during SWS revealed that ripples were coupled to the trough (active phase) of spindles with the mean SI angle close to 180° ( Fig. 4 G, SI phase = 174.40 ± 101.14°, SI strength = 0.52 ± 0.2, n = 1004). In a precise quantitative match with the experimental observation, in the model data also the ripples occurred mainly during the troughs of spindles corresponding to SI angle of 180°( Fig. 4 H, SI phase = 173.28 ± 53.74°, SI strength = 0.29 ± 0.16, n = 713).
Overall we found that in both experimental and model data, the ripples were coupled to the spindles and confirming the previous results ( Staresina et al., 2015 ;Varela and Wilson, 2020 ), the ripples mostly occurred at spindle troughs representing the active phase of the spindles. Previous findings also revealed that ripples with larger duration improved the memory while unlike the duration, the amplitude of the ripples were not associated with memory performance ( Fernández-Ruiz et al., 2019 ). We also found that the duration but not the amplitude of the ripples co-occurring with the spindles were larger compared to the isolated ripples. These results support the previous experimental observations indicating of the association between spindle-ripple coupling and memory performance ( Oyanedel et al., 2020 ;Siapas and Wilson, 1998 ).

The role of bidirectional connections in spindle-ripple coupling
We hypothesized that the spindle-ripple coupling is controlled by the long-range bidirectional connections between the hippocampal and cortical networks. Hence we next studied the spindle-ripple coupling in the model data as a function of the strength of the projections between the hippocampal and cortical networks. The percentages of the ripples cooccurring with the spindles did not depend substantially on the strength  LFPs during one spindle co-occurring with ripple events. The PFC and CA1 LFPs were filtered in the spindle (11-15 Hz) and ripple (100-250 Hz) frequency ranges respectively. B) Histogram of duration (Left) and amplitude (Right) of ripples co-occurring with spindles (red) and isolated ripples (blue) for the experimental data.
The mean values are shown by vertical dotted lines. C) Raw (Top) and filtered (Bottom) traces of the membrane potential of the cortical excitatory neuron (red) and CA1 excitatory neuron (black) during one spindle co-occurring with ripple events. The cortical and hippocampal traces were filtered in the spindle (12-16 Hz) and ripple (100-250 Hz) frequency ranges respectively. In this and all the subsequent figures and discussion, the membrane potential of the excitatory neurons of the model data was multiplied by − 1 to enable easy comparison with the experimental LFP data. D) The same as (B) for the model data. E) Histogram of duration (Left) and amplitude (Right) of ripples co-occurring with spindles (red) and isolated ripples (blue) for the experimental data F) Histogram of duration (Left) and amplitude (Right) of ripples co-occurring with spindles (red) and isolated ripples (blue) for the model data. G) Top, average of spindle peak-locked TFR of the hippocampal LFP (percentage change from pre-event baseline, n = 1004). Black curves show grand average filtered cortical LFP in the spindle frequency range (11-15 Hz) aligned to the spindle peak (time 0). Bottom, circular histogram of SI angles of the ripples relative to the spindle. SI angles were non-uniformly distributed ( p = 5.70 ×1 0 −5 , z = 9.75, Rayleigh test). H) Top, average of spindle peak-locked TFR of CA1 excitatory neuron ( n = 713) of the full model. Black curves show grand average filtered trace of cortical excitatory neuron in the spindle frequency range (12-16 Hz) aligned to the spindle peak (time 0). Bottom, circular histogram of SI angles of the ripples relative to the spindle in model data.
of long-range bidirectional connections between hippocampal and cortical networks. (Fig. S2).
The spindle-ripple coupling strength was an increasing function of the strength of the excitatory connection from the hippocampal to the cortical networks ( Fig. 5 A, r = 0.82, p = 0.002, Pearson correlation). For low values of the hippocampo-cortical excitatory connection, the coupling between the spindles and ripples was only due to the corticohippocampal connection, so that the depolarizing phase of the spindle (maximum at spindle trough, 180°) provided the required excitation for the generation of sharp wave ripples in the hippocampal network. As a result of these causal relationship, the phase of the coupling in this case was larger than 180°corresponding to the occurrence of the ripples after the spindle depolarizing phase ( Fig. 5 A Bottom Right). For larger values of hippocampo-cortical input, due to the influence of sharp wave ripples on the spindles, the phase of spindle-ripple coupling changed towards lower angles corresponding to the occurrence of ripples before the spindle trough. Ripples are accompanied by sharp waves which are large depolarizing deflections with a duration of the order of spindle cycles (one over the spindle frequency). Due to large hippocampo-cortical projection, cortex received a large depolarizing input at the time of the SWRs. When this depolarizing input to the cortex arrived at the time of the spindle occurrence, it can change the phase of the spindle towards more depolarizing phases which resulted in the increase of the spindleripple coupling strength and shift of the spindle-ripple coupling phase towards lower values ( Fig. 5 A Bottom Right). To validate this proposed mechanism, we omitted the hippocampal input to the thalamocortical network and applied a depolarizing input current (a 50 ms pulse with amplitude of 500 pA) resembling sharp waves to the cortical excitatory neurons during spindles. When the input current was applied at the hyperpolarizing phase of the spindle, it resulted in the phase shift to more depolarizing phases (Fig. S3). The applied input current at the depolarizing phase of the spindles did not change the phase of spindle oscillations considerably. These results explains how large enough hippocampo-cortical input at the time of SWR can change the phase of the spindle and increase the spindle-coupling strength with a mean spindle-ripple coupling phase before 180°.
The correlation between the strength of excitatory corticohippocampal connections and spindle-ripple coupling strength followed a non-monotonic inverted-U shape demonstrating that while for low and large values of cortico-hippocampal excitatory connection, the strength of the coupling between the thalamocortical spindles and hippocampal ripples was low, the maximum spindle-ripple coupling occurred at intermediate values of cortico-hippocampal excitatory connection strength ( Fig. 5 B, r = − 0.1, p = 0.75, Pearson correlation). At this intermediate value of the cortico-hippocampal strength, the ripples mostly occurred at 180°corresponding to the largest depolarizing time of the spindle ( Fig. 5 B Bottom Right) and the probability of the occurrence of the ripples at other phases (corresponding to less depolarization) was low since the cortico-hippocampal input was not large enough to induce the required depolarization essential for ripple generation. In this case since most of the ripples occurred at a consistent phase of the spindle (close to 180°), the spindle-ripple coupling was large. For larger values of corticohippocampal strength, the depolarizing input to the hippocampus was large enough even at less depolarized phases of the spindle. Hence in this case, the ripples can occur at a wider range of the spindle phase (more variable phases) resulting in the decreasing of the spindle-ripple coupling strength.
To further check the dependence of spindle-ripple coupling on the bidirectional connections between hippocampal and cortical networks, we next varied both of the connections simultaneously and computed the spindle-ripple coupling measures (Fig. S4). The coupling strength was mostly larger at a region corresponding to large values of hippocampo-cortical excitatory connection and intermediate values of cortico-hippocampal excitatory connections. Consistent with the proposed mechanism explained above, the mean angle at these regions with large coupling strength was around 180°. We next investigated the effect of long-range hippocampo-cortical and cortico-hippocampal excitatory to inhibitory inputs on the spindleripple coupling (Fig. S5). The spindle-ripple coupling strength increased by increasing the hippocampo-cortical excitatory to inhibitory input but it did not change considerably by increasing the cortico-hippocampal input. The spindle-ripple coupling phase also did not depend substantially on the cortico-hippocampal or hippocampo-cortical excitatory to inhibitory input.
These results shed light on our understanding of the mechanisms behind the role of bidirectional long-range connections between the cortical and hippocampal networks in the temporal coordination of spindle and ripples which was suggested to be important for successful memory consolidation ( Maingret et al., 2016 ;Siapas and Wilson, 1998 ). Our results indicating the importance of the hippocampal to cortical connections for spindle-ripple coupling is consistent with a recent finding that revealed the impairment of the modulation of ripples by spindles as a result of the silencing the projection from hippocampus to medial prefrontal cortex ( Binder et al., 2019 ). Our findings also indicate the reciprocal influence of thalamocortical spindles and hippocampal ripples.

Occurrence of the thalamic spindles at up state phase of the SO
We next investigated the temporal coordination between SOs and spindles in both experimental and model data. During SWS, the majority of the spindles co-occurred with SOs (for 92.9% of the spindles, SO peak occurred within a ± 1.5 s interval around their maximum peak). Both duration and amplitude of the spindles coupled to SO were larger compared to the spindles occurring alone. ( Fig. 6 A, p = 3.96 ×1 0 −4 and p = 1.22 ×1 0 −2 for duration and amplitude respectively, Wilcoxon rank sum, mean ± SD: duration = 1.04 ± 0.59 s and amplitude = 3.25 ± 0.94 SD, n = 669 for spindles coupled to SOs and duration = 0.78 ± 0.46 s and amplitude = 2.81 ± 0.76 SD, n = 51 for isolated spindles). In a quantitative match to the experimental data, the majority of the spindles in the model data were also coupled to the SOs (89%,). In agreement with the experiment, the duration and amplitude of the spindles co-occurred with the SO were also larger than the isolated spindles in the model data ( Fig. 6 B, p = 0.70 ×1 0 −6 and p = 1.13 ×1 0 −3 for duration and amplitude respectively, Wilcoxon rank sum, mean ± SD: duration = 0.69 ± 0.15 s and amplitude = 2.87 ± 0.55 SD, n = 635 for spindles coupled to SOs and duration = 0.56 ± 0.17 s and amplitude = 2.62 ± 0.79 SD, n = 78 for isolated spindles).
Moreover, PAC analyses revealed that during SWS, the power in the frequency range of thalamic spindles was dominant at the up state phase of the SO before the SO trough with the mean SI angle equal to 141.42 ± 50.5° ( Fig. 6 C, p = 2.49 ×1 0 −15 , z = 33.48, Rayleigh test. SI strength = 0.48 ± 0.22). Similarly, in the model data, the mean SI angle of SO_spindle coupling was 154.86 ± 11.72°so that spindles mainly occurred before the SO trough ( Fig. 6 D, SI strength = 0.76 ± 0.18).
We hypothesized that the increase in the duration and amplitude of the spindles coupled to the SO compared to the isolated spindles was due to the large cortical drive to the thalamus during the up state phase of the SO. To test this hypothesis, we next computed the amplitude and duration of the spindles as a function of the strength of long-range cortico-thalamic excitatory connections in the model data. Consistent with our hypothesis, both duration and amplitude of the spindles increased monotonically by increasing the strength of cortical excitatory projections to the thalamic network ( Fig. 6 E, r = 0.79, p = 2.53 ×1 0 −4 and r = 0.77, p = 5.15 ×1 0 −4 for duration and amplitude respectively, Pearson correlation). The percentages of the spindles co-occurring with SOs did not correlate significantly with the strength of the long-range corticothalamic connections (Fig. S2).
Overall, our model replicates the experimental observation showing large probability of SO and spindle co-occurrence and larger amplitude and duration of the spindles coupled to SOs compared to the isolated spindles that can be explained by long-range cortico-thalamic connections.   7. Temporal relationship between SOs and ripples. A) Raw (Top) and filtered (Bottom) traces of PFC (red) and CA1 (black) LFPs during one SO co-occurring with ripple events. The PFC and CA1 LFPs were filtered in the SO (0.5-2 Hz) and ripple (100-250 Hz) frequency ranges respectively. B) Histogram of the duration (Left) and amplitude (Right) of ripples co-occurring with SOs (red) and isolated ripples (blue) for experimental data. C) Raw (Top) and filtered (Bottom) traces of the membrane potential of the cortical excitatory neuron (red) and CA1 excitatory neuron (black) during one SO co-occurring with ripple events. The cortical and hippocampal traces were filtered in the SO (0.5-2 Hz) and ripple (100-250 Hz) frequency ranges respectively. D) Histogram of the duration (Left) and amplitude (Right) of ripples co-occurring with SOs (red) and isolated ripples (blue) for model data. E) The circular distribution of the SO phase at the time of the ripples co-occurring with SO for experimental data ( p = 6.47 ×1 0 −3 , z = 9.61, Rayleigh test, n = 657). F) The circular distribution of the SO phase at the time of the ripples co-occurring with SO for model data ( p < 1 0 −20 , z = 692.91, Rayleigh test, n = 5447).

Temporal coordination of SOs and ripples
We next investigated the temporal association between hippocampal ripples and cortical SOs. Unlike the strong spindle-ripple and SO-spindle coupling, only 54.4% of the ripples were coupled to SOs (SO peak occurred within a ± 0.5 s interval around ripple maximum peak, Fig. 7 A). The duration of ripples coupled to SOs was larger than the duration of the ripples occurring alone ( Fig. 7 B, p = 3.88 ×1 0 −10 , Wilcoxon rank sum, mean ± SD: duration = 0.04 ± 0.01 s, n = 657 for ripples coupled with SOs and duration = 0.03 ± 0.01 s, n = 550 for isolated ripples). However, the amplitude of the ripples co-occurring with SOs was lower than the isolated ripples ( Fig. 7 B, p = 1.22 ×1 0 −3 , Wilcoxon rank sum, mean ± SD: amplitude = 6.36 ± 2.52 SD, n = 657 for ripples coupled with SOs and amplitude = 7.18 ± 3.01 SD, n = 550 for isolated ripples). In a very good agreement with the experimental data, in the model data also a considerable percentage of the hippocampal ripples were not co-occurred with the cortical SOs ( Fig. 7 C, 37.3%). In addition, the duration of the ripples coupled to SOs was larger than the duration of the isolated ripples while the amplitude of the ripples co-occurring with SOs was lower than the amplitude of the ripples occurring alone ( Fig. 7 D, p = 3.88 ×1 0 −6 and p = 1.25 ×1 0 −5 for duration and amplitude respectively, Wilcoxon rank sum, mean ± SD: duration = 0.04 ± 0.002 s and amplitude = 6.69 ± 0.69 SD, n = 5447 for ripples coupled to SOs and duration = 0.03 ± 0.003 s and amplitude = 7.37 ± 0.89 SD, n = 3246 for isolated ripples).
To find the SO phase at which the ripples co-occurring with SO preferentially occurred, we next obtained the SO phase at the time of the ripples for both experimental and model data using Hilbert transform. The distribution of the SO phase at the time of the ripples was non-uniformly distributed with the mean phase of 206.53 ± 116.40°corresponding to the up state phase of the SO after the SO trough ( Fig. 7 E). However, the phases were less concentrated compared to the SO-spindle or spindleripple coupling phase confirming a weaker coupling between SOs and ripples (resultant vector length = 0.13, 0.22, 0.28 for SO-ripple, spindleripple and SO-spindle respectively). Consistently, in the model data, the distribution of the SO phases at the time of the ripples co-occurring with SO was also non-uniformly distributed revealing preferential occurrence of ripples at the up state phase of the SO slightly after the SO trough ( Fig. 7 F,mean angle = 196.44 ± 82.69°). However consistent with the results obtained during SWS, in the model data also the SO phases at the time of the ripples were less concentrated compared to SO-spindle and spindle-ripple coupling phase (resultant vector length = 0.35, 0.64, 0.98 for SO-ripple, spindle-ripple and SO-spindle respectively).
For both experimental and model data, the SO phases at the time of the spindle event (maximum peak) were dominantly in the upper left quadrant (90-180°), i.e., before the SO trough (Fig. S6, 108.97 ± 100.95°f or the experimental data and 137.53 ± 53.25°for the model data). Hence comparing with our previous results showing that ripples mostly occurred after the SO trough ( Fig. 7 E & F), we found that ripples mostly occurred at later phase of SO compared to spindles.
Next, we investigated the role of the strength of long-range excitatory connections between the hippocampal and cortical networks on the SO-ripple coupling in the model. The percentage of the ripples co-occurring with SOs was positively correlated with the strength of cortico-hippocampal connections but it was not correlated with the strength of hippocampo-cortical connections (Fig. S2).
The resultant vector length of the SO phases at the time of the ripples which measures the strength of the phase concentration, increased by increasing the strength of hippocampo-cortical excitatory connections ( Fig. 8 A). The increase of the hippocampo-cortical excitatory connection changed the mean angle from around 150 to 200°. The resultant vector length and mean phase also depended non-monotonically on the strength of the excitatory projections from the cortex to the hippocampus ( Fig. 8 B). The maximum resultant vector length occurred at an intermediate value of the cortico-hippocampal excitatory connections while for low and large values of the cortico-hippocampal excitatory connections the phases were less concentrated.
To further check the dependence of SO-ripple coupling on the bidirectional connections between hippocampal and cortical networks, we next varied both of the connections simultaneously and computed the SO-ripple coupling measures (Fig. S7). The resultant vector length was larger at large values of hippocampo-cortical excitatory connection and intermediate values of cortico-hippocampal excitatory connections.
The mean angle also depended more substantially on the corticohippocampal connection than hippocampo-cortical connection so that by increasing the cortico-hippocampal connection strength, independent of the hippocampo-cortical connection strength, it decreased from the phases corresponding to SO up-to-down transition (after SO trough) towards the phases corresponding to down-to-up transition (before SO trough).
Overall these results suggest that the long-range bidirectional excitatory connections between thalamic/hippocampal and cortical networks are the key ingredients responsible for the coupling among SOs, spindles and ripples. In addition, ripples are more strongly coupled to spindles than SOs.

Discussion
Using a minimal mean field model of hippocampo-cortico-thalamic network, we investigated the mechanism underlying the generation and coupling of the SWS rhythms. In our model, the isolated CA3-CA1 network can spontaneously generate ripples with characteristics in line with the experimental observations. Especially we showed analytically, by simulation and by CA1 LFP analyses during SWS that the ripple frequency was linked to the strength of the input from CA3 corresponding to the sharp wave amplitude. In the full model, the hippocampal ripples interacted with the SOs and spindles that were generated in the cortical and thalamic networks respectively. Consistent with the experiment, by employing phase amplitude coupling analyses we found that the spin-dles and ripples were strongly coupled so that ripples mostly occurred during troughs of the spindles. This large coupling in the model can be explained by the long-range excitatory connections especially the connection from the hippocampus to the cortex. Further, we showed that during SWS as well as in the model data, the duration of the ripples cooccurring with spindle was larger than the isolated ripples. Our model was also able to reproduce the experimentally observed SO-spindle as well as SO-ripple coupling and the differential properties of spindles and ripples co-occurring with SOs and the isolated ones.
We showed SWRs emerged in a simple mean field firing rate model with no synaptic rate equation. In line with previous studies ( Evangelista et al., 2020 ;Levenstein et al., 2019 ), in this model the sharp waves were generated in the CA3 network due to adaptation of CA3 excitatory neurons. Sharp dependence of the firing rate of CA1 neurons on the membrane potential was also the key ingredient to reproduce high frequency activity corresponding to ripples. As the CA1 model was simplified consisting of only 2 identical groups of neurons, we further solved analytically the corresponding dynamical system subject to excitatory inputs representing the CA3 input and showed that the frequency of the oscillations was a monotonically increasing function of the excitatory input. Our model further predicted that the changes in the amplitude of the SWR might be due to changes in the local excitatory to excitatory connection within CA3.

Temporal coordination of SOs, spindles and ripples
Our results that SWRs preferably occurred during up state phase of SO mostly after SO trough is consistent with the previous findings ( Helfrich et al., 2019 ;Maingret et al., 2016 ;Peyrache et al., 2011 ). However, studies also reported the main occurrence of SWRs before the trough of the SO corresponding to down-to-up transition phase of SO (Battaglia et al., 2004) or at SO peak corresponding to the down state ( Varela and Wilson, 2020 ).
Although previous studies have provided evidence for precise temporal coupling of ripples to spindles ( Clemens et al., 2011 ;Helfrich et al., 2019 ;Jiang et al., 2019 ;Mölle et al., 2006 ;Ngo et al., 2020 ;Siapas and Wilson, 1998 ;Sirota et al., 2003 ;Staresina et al., 2015 ;Varela and Wilson, 2020 ) the mechanism underlying this coupling remains to be elucidated. It has been shown that the inhibition of CA1 parvalbuminpositive cells disrupted the spindle-ripple coupling ( Xia et al., 2017 ). This observation cannot be tested in our mean field model with only one identical group of inhibitory CA1 neurons. A recent study reported that silencing the projections from hippocampus to medial prefrontal cortex  Fig. 7 ). p = 5.98 ×1 0 −28 , 1.47 ×1 0 −31 , 2.88 ×1 0 −19 , z = 62.15, 289.03, 42.34 for α CX−HP ee = 0.6, 0.9, 1.2 respectively, Rayleigh test. Bottom, the resultant vector length (Left) and mean angle (Right) of SO phases at the time of the ripples co-occurring with SOs as a function of α CX−HP ee . reduced the coupling of SWRs to spindle troughs ( Binder et al., 2019 ). Consistent with this experimental observation, we found that the projection from the hippocampus to the cortex contributes to the spindle-ripple coupling. Our model further predicts an inverted U-shape relationship between the spindle-ripple coupling and strength of the projection from cortex to the hippocampus.
While previous computational models did not account for the simultaneously generation and coupling of all the 3 main SWS rhythms, i.e. ripples, spindles and SOs, previous models of thalamocortical networks were developed to generate simultaneously both SOs and spindles ( Bazhenov et al., 2002 ;Bonjean et al., 2011 ;Cona et al., 2014 ;Destexhe and Sejnowski, 2003 ;Hashemi et al., 2019 ;Schellenberger Costa et al., 2016 ). In addition, hippocampo-cortical network was modeled to investigate the interaction of SOs and ripples ( Taxidis et al., 2013 ). However, a model to investigate the coupling of spindles and ripples is lacking. Recently a model of hippocampo-corticothalamic network was proposed to explain the interaction of SOs and SWRs ( Sanda et al., 2021 ). However, spindles were not generated and studied in this model. In our model, SOs, spindles and SWRs are generated by local interactions between excitatory and inhibitory neurons of cortical, thalamic and hippocampal networks respectively. In the full model, long-range multi-regional projections regulate the temporally precise interaction of the SOs, spindles and SWRs so that this simple model can simultaneously reproduce the experimentally observed coupling of ripples to spindle troughs as well as both ripples and spindles to the up state phase of SO.
Manipulating studies revealed the role of spindle-ripple coupling in successful memory consolidation ( Binder et al., 2019 ;Latchoumane et al., 2017 ;Maingret et al., 2016 ;Xia et al., 2017 ). Our results that the SO-spindle and spindle-ripple coupling were stronger than the SO-ripple coupling highlights the important mediating role of the spindles in the interaction between cortex and hippocampus required for memory consolidation.

Differential properties of coupled and isolated spindles and ripples
While our model is deterministic with no intrinsic or external noise, the variability in the duration and amplitude of the SOs, spindles and ripples emerges as a result of nonlinear interaction of thalamic, cortical and hippocampal oscillators indicating a chaotic dynamic with deterministic noise ( Cortes et al., 2013 ;Ghorbani et al., 2012 ;van Vreeswijk and Sompolinsky, 1996 ). We showed that in both experimental and model data the amplitude and duration of the spindles coupled to the SOs or ripples were larger compared to the isolated spindles. We further showed that in the model, the amplitude and duration of the spindles were controlled by the cortico-thalamic projections. The association between spindle amplitude and duration with the memory consolidation was reported in previous studies ( Fogel et al., 2007 ;Ngo et al., 2013 ). In line with these findings, we found that spindles with larger amplitude and duration were associated with stronger multi-regional interactions. We also found that in both experimental and model data, the majority of the spindles co-occurred with SOs and ripples.
Recent finding revealed lower amplitude of ripples coupled to spindles compared to the ripples uncoupled to the spindles ( Yang et al., 2019 ). We also found that the ripples coupled to the spindles or SOs had lower amplitude compared to the isolated ripples. In addition, our results revealed that in both experiment and model, the duration of the ripples co-occurred with spindles or SOs was larger than the isolated ripples. It has been recently reported that the prolongation of the ripples improved the memory performance ( Fernández-Ruiz et al., 2019 ). Our findings indicate that longer ripples are more strongly coupled to the spindles and SOs involving in the interaction between hippocampus and cortex which is critical for memory consolidation. These results confirm the results of a recent study that reported greater spindle power during ripples with longer duration in humans ( Ngo et al., 2020 ). Our model and experimental observation suggest that the large excitatory in-put from the cortex to the thalamus during SO up state results in longer spindles which in turn are linked to longer ripples. Longer duration of spindles and ripples provide longer windows of opportunity for crossregional synchronization and information transfer.

Conclusion
In conclusion, while previous work was mostly restricted to the modeling of hippocampus or cortex, we developed a full hippocampocortico-thalamic model to investigate the mechanism underlying the precise temporal coupling of SOs, spindles and ripples simultaneously. We found that large cortico-thalamic input during the active phase of SO induced spindles with larger duration and amplitude which in turn coordinated more strongly with hippocampal ripples due to bidirectional projections between cortex and hippocampus. Our model and experimental findings during SWS shed light on the mediating role of spindles in the communication between cortex and hippocampus by providing a key processing time window for synaptic plasticity required for memory consolidation.