Abstract
Sleep is critical for regulation of synaptic efficacy, memories, and learning. However, the underlying mechanisms of how sleep rhythms contribute to consolidating memories acquired during wakefulness remain unclear. Here we studied the role of slow oscillations, 0.2–1 Hz rhythmic transitions between Up and Down states during stage 3/4 sleep, on dynamics of synaptic connectivity in the thalamocortical network model implementing spike-timing-dependent synaptic plasticity. We found that the spatiotemporal pattern of Up-state propagation determines the changes of synaptic strengths between neurons. Furthermore, an external input, mimicking hippocampal ripples, delivered to the cortical network results in input-specific changes of synaptic weights, which persisted after stimulation was removed. These synaptic changes promoted replay of specific firing sequences of the cortical neurons. Our study proposes a neuronal mechanism on how an interaction between hippocampal input, such as mediated by sharp wave-ripple events, cortical slow oscillations, and synaptic plasticity, may lead to consolidation of memories through preferential replay of cortical cell spike sequences during slow-wave sleep.
SIGNIFICANCE STATEMENT Sleep is critical for memory and learning. Replay during sleep of temporally ordered spike sequences related to a recent experience was proposed to be a neuronal substrate of memory consolidation. However, specific mechanisms of replay or how spike sequence replay leads to synaptic changes that underlie memory consolidation are still poorly understood. Here we used a detailed computational model of the thalamocortical system to report that interaction between slow cortical oscillations and synaptic plasticity during deep sleep can underlie mapping hippocampal memory traces to persistent cortical representation. This study provided, for the first time, a mechanistic explanation of how slow-wave sleep may promote consolidation of recent memory events.
Introduction
Sleep constitutes approximately one-third of our lives. During sleep, human and animal brains are primarily decoupled from sensory input. Nevertheless, brain activity remains high and is manifested in various types of sleep-specific synchronized rhythms (Frank, 2012). It was proposed that sleep plays a critical role in memory consolidation, in which newly acquired information during wakefulness is transformed into long-term memory during sleep (Born, 2010). Specifically, slow-wave sleep (SWS) has been shown to facilitate consolidation of declarative, hippocampus-dependent memory (Gais et al., 2000; Marshall et al., 2006; Marshall and Born, 2007; Barnes and Wilson, 2014), whereas rapid eye movement (REM) sleep can benefit procedural, hippocampus-independent memory (Born, 2010).
SWS, often referred to as deep sleep, consists of stage 3/4 of non-REM (NREM) sleep. Slow-wave brain activity, observed during SWS in humans and animals (Steriade et al., 1993), is characterized by large-amplitude slow oscillations (0.2–1 Hz) in the electroencephalograph (Blake and Gerard, 1937), reflecting alternating periods of active (Up) and silent (Down) states in the intracellular recordings of thalamic and cortical neurons (Steriade et al., 2001). Previous studies have pointed to the importance of SWS in memory consolidation. In humans, it has been shown that enhancing SWS by transcranial electrical stimulation can enhance memory consolidation (Marshall et al., 2006). In rodents, during SWS, temporally ordered firing sequences related to a recent experience replayed in both hippocampus and neocortex synchronously (Ji and Wilson, 2007; Mehta, 2007). Such sequence replay has been proposed to be a neural substrate of memory consolidation (Barnes and Wilson, 2014) and is believed to result in synaptic changes in the neocortex responsible for integration of memory representations (Schwindel and McNaughton, 2011). Cognitive models based on experimental studies suggested a distributed role for memory representation, and the interaction between the hippocampus and cortex is proposed as a mechanism for consolidation (McClelland et al., 1995; Rogers and McClelland, 2006). Nevertheless, specific neuronal mechanisms of replay and how spike sequence replay leads to synaptic changes that underlie memory consolidation are hard to demonstrate in animals in vivo and are still poorly understood.
To address these questions, we examined how synaptic connectivity in a biologically realistic network model of the thalamocortical system were changed by and in return affected the spatiotemporal pattern of sleep slow oscillations (Bazhenov et al., 2002, 2011). Spike-timing-dependent plasticity (STDP) was implemented to the cortical network to regulate synaptic efficacy. Sparse external stimulation was applied on a few cortical pyramidal cells to mimic hippocampal input during SWS. We found that even weak and spatially localized input influenced the spatiotemporal pattern of slow oscillations and led to a persistent change of synaptic efficacy between neurons. These synaptic changes remained even after the input was removed and promoted the specific order of neuronal firing, similar to sequence replay. Our study proposes a neuronal mechanism of how the interaction between cortically generated slow waves and sparse external input, simulating input from hippocampal formation, may lead to reorganization of synaptic strength during stage 3/4 sleep, therefore consolidating recent memories.
Materials and Methods
Thalamocortical model
The network model (Fig. 1) incorporated a model of thalamus with thalamic relay (TC) and reticular (RE) neurons, as well as a model of cortical column with pyramidal neurons (PY) and inhibitory interneurons (IN; Bazhenov et al., 2002, 2011). All neurons were modeled based on the Hodgkin–Huxley kinetics. The units and description of parameters were summarized in Table 1.
The units and description of the parameters used in the model
Intrinsic currents: thalamus
TC and RE cells were modeled as a single compartment that included voltage- and calcium-dependent currents described by Hodgkin–Huxley kinetics (Bazhenov et al., 2002):
where Cm = 1u F/cm2 is the membrane capacitance, gL is the leakage conductance (0.01 mS/cm2 for TC cell and 0.05 mS/cm2 for RE cell), and EL is the reversal potential (−70 mV for TC and −77 mV for RE). Isyn is the sum of synaptic currents. Iint is the sum of active intrinsic currents that include a fast sodium current, INa, a fast potassium current IK (Traub and Miles, 1991), a low-threshold Ca2+ current, IT, for RE (Huguenard and Prince, 1992) and TC (Huguenard and McCormick, 1992) cells, and a potassium leak current, IKL, where IKL = gKL(V − EKL), gKL is the potassium leakage conductance, and EKL is the potassium reversal potential. A hyperpolarization-activated cation current, Ih (Huguenard and McCormick, 1992), was included in TC neurons. Therefore, the sum of active intrinsic currents for RE and TC are
All the voltage-dependent ionic currents Ij have the general form:
where gj is the maximal conductance, and Ej is the reversal potential. The dynamic of gating variables are described as
where x = m or h. QT is a temperature-related term. Usually, QT = Q((T−23)/10) = 2.9529, with Q = 2.3, and T = 36. The expressions of transition rates for all currents were given in a previous study (Bazhenov et al., 1998).
Fast sodium current INa.
The INa current is described by
where gNa is the maximal sodium conductance (90 mS/cm2 for TC, 100 mS/cm2 for RE). The expressions of transition rates for activation variable m and inactivation variable h are described by
Fast potassium current IK.
The IK current is described by
where gK is the maximal potassium conductance (90 mS/cm2 for TC, 100 mS/cm2 for RE). The expressions of transition rates are
Low-threshold calcium current IT.
The IT current is described by
where gT is the maximal low-threshold calcium conductance (2.3 mS/cm2 for TC and RE). The expressions of transition rates for RE cells (Huguenard and Prince, 1992; Destexhe et al., 1996b) are
where Qm = 5((T−24)/10) = 6.8986, and Qh = 3((T−24)/10) = 3.7372, with T = 36.
The expressions of transition rates for TC cells (Huguenard and McCormick, 1992; Destexhe et al., 1996a) are
where Qm = 3.55((T−24)/10) = 4.5738, and Qh = 3((T−24)/10) = 3.7372, with T = 36.
Hyperpolarization-activated cation current Ih.
The Ih current for TC cells is described by (Destexhe et al., 1996a)
where k = 2.2 and Eh = −40 mV. gh is the maximal hyperpolarization-activated cation conductance (0.017 mS/cm2 for TC). The fraction of the channels in the closed (C = 1-O-OL), opened (O), and locked (OL) forms were calculated by
where k1 = 7.9012 × 107 mm−4ms−1, k2 = 0.004 ms−1, k3 = 0.1 ms−1, and k4 = 0.001 ms−1 are constant rates. The calcium dependence was based on higher-order kinetics involving a regulation factor P. The binding of calcium molecules with an unbound form of the regulation factor P0 (=1 − P1) leads to the bound form P1. α and β are the voltage-dependent transition rates (Huguenard and McCormick, 1992) as follows:
with
Potassium leak current IKL.
The IKL current is described by
where gKL is the potassium leak conductance (0.005 mS/cm2 for RE and 0.03 mS/cm2 for TC), and EKL = −95 mV.
Calcium dynamics.
For both RE and TC cells, the calcium dynamics is described by
where [Ca]∞ is 2.4 × 10−4 mm, indicating equilibrium calcium concentration, A = 5.1819 × 10−5 mm cm2/(ms μA), and τ = 5 ms.
Intrinsic currents: cortex
The cortical PY and IN cells were two-compartment models described by Hodgkin–Huxley kinetics (Bazhenov et al., 2002):
where Cm is the membrane capacitance, gL is the leakage conductance of the dendrite compartment, and gSD = 1/(R*Sdend) and gDS = 1/(R*Ssoma) are the conductances between the axosomatic and dendritic compartment, where R = 10 MΩ is the resistance between compartments. VD and VS are the membrane potentials of dendritic and axosomatic compartments. In this model, we set the capacitance of axosomatic compartment to be 0 (Chen et al., 2012). This model was first proposed (Mainen and Sejnowski, 1996) as a reduction of a multicompartmental pyramidal cell model, based on the assumption that the current dynamics in the axosomatic compartment are fast enough to ensure that VS is always at equilibrium state, as defined by the equation above. IDint and ISint are the sums of active intrinsic currents in the dendritic and axosomatic compartments and can be described as follows:
In PY cells, IDint included a fast sodium current, INa, a persistent sodium current, INa(p), a slow voltage-dependent non-inactivating potassium current, IKm, a slow Ca2+-dependent K+ current, IKCa, a high-threshold Ca2+ current, IHVA, and a potassium leak current, IKL, in the dendritic compartment. ISint included a fast sodium current, INa, a persistent sodium current, INa(p), and a fast delayed rectifier potassium current, IK, in the axosomatic compartment. The persistent sodium current, INa(p), was included in the axosomatic and dendritic compartments of PY cells to increase bursting propensity. IN cells have the same intrinsic currents as those in PY except that INa(p) is not included. The ratio of dendritic area to axosomatic area was ρ = 165 for PY and ρ = 50 for IN.
Fast sodium current INa.
The INa current is described by
where gNa is the maximal sodium conductance (3000 mS/cm2 for soma, 0.8 mS/cm2 for dendrite). ENa = 50 mV is the reversal potential of sodium. The expressions of transition rates for fast sodium current are
and
Persistent sodium current INa(p).
The INa(p) current is described by
where gNa(p) is the maximal persistent sodium conductance (15 mS/cm2 for soma, 2.5 mS/cm2 for dendrite). The expressions of transition rates are
Slow voltage-dependent noninactivating potassium current IKm.
The IKm current is described by
where gKm is the maximal slow voltage-dependent noninactivating potassium conductance (0.02 mS/cm2). EK = −90 mV is the reversal potential of potassium. The expressions of transition rates are
Slow calcium-dependent potassium current IKCa.
The IKCa current is described by
where gKCa is the maximal slow Ca2+-dependent K+ conductance (0.3 mS/cm2). The expressions of transition rates are
where [Ca2+]i is the intracellular calcium concentration.
High-threshold calcium current IHVA.
The IHVA current is described by
where gHVA is the maximal high-threshold Ca2+ conductance (0.02 mS/cm2). ECa = 140 mV is the reversal potential of calcium. The expressions of transition rates are
Fast delayed rectifier potassium current IK.
The IK current is described by
where gK is the maximal potassium conductance (200 mS/cm2). The expressions of transition rates are
Calcium dynamics.
For each cell, the calcium dynamics is described by a simple first-order model:
where [Ca]∞ is 2.4 × 10−4 mm, indicating equilibrium calcium concentration, A = 5.1819 × 10−5 mm cm2/(ms μA) and τ = 5 ms.
Synaptic currents
The equations for all synaptic currents are given in previous studies (Timofeev et al., 2000; Bazhenov et al., 2002). GABAA, AMPA, and NMDA synaptic currents are given by
where gsyn is the maximal conductance, [O] is the fraction of open channels, and Esyn is the reversal potential (EAMPA = 0 mV, ENMDA = 0 mV, and EGABAA = −70 mV). Here f(V) = 1 for GABAA and AMPA receptors, but for NMDA receptor, it is a voltage-dependent sigmoid function f(V) = 1/(1 + exp(−(V − Vth)/σ)), where Vth = −25 mV, and σ = 12.5 mV (Traub and Miles, 1991; Chen et al., 2012). The maximal conductances for each synapse were gGABAA(RE−TC) = 0.15 μS, gGABAA(RE−RE) = 0.125 μS, gAMPA(TC−RE) = 0.35 μS, gAMPA(TC−PY) = 0.1 μS, gAMPA(TC−IN) = 0.05 μS, gAMPA(PY−PY) = 0.08 μS, gNMDA(PY−PY) = 0.006 μS, gAMPA (PY−IN) = 0.08 μS, gNMDA(PY−IN) = 0.005 μS, and gGABAA(IN−PY) = 0.25 μS. The maximal conductances between different cortical layers were gAMPA(PY−PY) = 0.05 μS and gNMDA(PY−PY) = 0.005 μS. The faction of open channel [O] is calculated according to the equation
where θ(x) is the Heaviside function and t0 is the time of receptor activation. The amplitude and duration of the neurotransmitter pulse were A = 0.5 and tmax = 0.3 ms. The rate constants were α = 1.1 ms and β = 0.19 ms for AMPA receptors, α = 1 ms and β = 0.0067 ms for NMDA receptors, and α = 10.5 ms and β = 0.166 ms for GABAA receptors.
In the intracortical connections, GABAA, and AMPA synaptic currents are modified by a short-term depression term D (Timofeev et al., 2000):
where D represents the amount of available “synaptic resources.” D was calculated according to the following interactive scheme (Timofeev et al., 2000):
where Δt is the time interval between nth and (n + 1)th spike, τ (700 ms for AMPA and GABAA) is the time constant of recovery of the synaptic resources, and U (0.073 for GABAA, 0.07 for AMPA) is the fraction of resources used per action potential.
GABAB synaptic current is given by the following (Destexhe et al., 1996a):
where [R] is the fraction of activated receptors, [G] is the concentration of G-proteins, gGABAB = 0.04 μS, and EK = −95 mV is the potassium reverse potential. The rate constants were K1 = 0.52 mm−1ms−1, K2 = 0.0013 ms−1, K3 = 0.098 ms−1, K4 = 0.033 ms−1, and K = 100 μm4.
Network geometry
The cortical column included three layers: layer II/III/IV, layer V, and layer VI. Each layer was organized in a one-dimensional chain of n PY neurons and n/4 INs, where n = 200 in this model. The radii of connections between cortical neurons were RAMPA(PY−PY) = 5, RAMPA(PY−IN) = 2, and RGABAA(IN−PY) = 5. We varied RAMPA(PY−PY) = 5–7 in layer V to confirm that layer V is the most likely site of Up-state initiation across layers for Figure 2. The thalamus model consisted of core and matrix systems (Jones, 2001). Each system consisted of n/4 TC and n/4 RE cells. The radii of connections between neurons within each system were RGABAA(RE−TC) = 6, RGABAB(RE−TC) = 6, RGABAA(RE−RE) = 4, and RAMPA(TC−RE) = 4. The core TC neurons projected to the middle layers of the cortex (layer II/III/IV) (RAMPA(TC−PY) = 5 and RAMPA(TC−IN) = 1), whereas the matrix TC neurons projected to the superficial layer (layer I) whose cell bodies were in the layer V (Jones, 2001) (RAMPA(TC−PY) = 20 and RAMPA(TC−IN) = 4). The matrix thalamocortical pathway had more diffuse connections than the core (Jones, 2001). The radii of connections from cortical layer VI to the core thalamus system were RAMPA(PY−TC) = 6 and RAMPA(PY−RE) = 4 and from cortical layer V to the matrix thalamus system were RAMPA(PY−TC) = 10 and RAMPA(PY−RE) = 5.
Mechanism of initiation and termination of Up states
Slow oscillation in the model was driven by intracortical dynamics. Spontaneous summation of miniature EPSPs (mEPSPs) through AMPA-mediated synapses among pyramidal neurons triggered Na+ spikes in some (randomly located) PY neurons, which then initiated spikes in their postsynaptic neighbors (Timofeev et al., 2000; Bazhenov et al., 2002; Chen et al., 2012). In the model, we set different state of the computer system to change the random seed for different initial conditions. The spread of activity within a cortical layer and between cortical layers was mediated by the lateral excitatory connections, which limited the degree of synchronization at active state onset.
The arrival times of spontaneous mEPSPs were modeled by Poisson processes (Stevens, 1993), with time-dependent mean rate μ = (2/(1 + exp(−(t − t0)/F))−1)/250 (Bazhenov et al., 2002), where t0 is a time instant of the last presynaptic spike (Timofeev et al., 2000). The mEPSP frequency (F) and amplitude (A) within cortical layers were FPY−PY = 50, FPY−IN = 80, APY−PY = 0.2 μS, and APY−IN = 0.05 μS. The mEPSP frequency and amplitude between different cortical layers were FPY−PY = 80 and APY−PY = 0.05 μS.
Different mechanisms have been proposed to be involved in termination of Up states. In our previous study of the slow oscillations based on recordings from isolated cortical slab (Timofeev et al., 2000; Bazhenov et al., 2002), we proposed that termination of Up states depends on (1) slow adaptation currents such as Ca2+-dependent K+ current (IKCa) and (2) synaptic depression. In another study, it has been proposed that termination may involve activation of the Na+-dependent K+ current (IKNa) (Compte et al., 2003). Finally, another model of slow oscillation proposed that “depolarization-activated potassium currents and synaptic depression terminate the Up state” without specifying the origin of the depolarization-activated potassium current (Hill and Tononi, 2005). More recently, we extended our model and proposed that activation of the inhibitory interneurons can be involved in Up-states termination (Chen et al., 2012); this model was able to explain higher synchrony of Up-state termination versus initiation as it was reported in vivo. Our current model incorporate all the properties mentioned above except for IKNa.
Spike-timing-dependent synaptic plasticity
Synaptic plasticity leading to enhancement or depression in synaptic strength is widely thought to underlie learning and memory storage in the brain. Among the various forms of synaptic plasticity, we used STDP, which adjusts the synaptic efficacy based on the relative timing of the presynaptic and postsynaptic spikes. The basic STDP is described by (Song et al., 2000)
where Δt, the spike time of postsynaptic neurons i minus the spike time of presynaptic neuron j, determines whether the synaptic weight is strengthened by long-term potentiation or weakened by long-term depression. A+ and A− determines the maximum amounts of synaptic modification. Here, we set A+ = A− = 0.0005, τ+ = τ− = 20 ms.
In this study, excitatory synaptic connections within cortex layer were implemented with STDP:
where gmax is the maximal synaptic conductance. If Δt is positive, the conductance of the excitatory-to-excitatory synapse is potentiated. If Δt is negative, the conductance is depressed. To limit the synaptic changes, we assumed that the synaptic efficacy stays within an interval [0, gmax]. In text, the synaptic weight was then normalized with gmax.
Spontaneously released glutamate vesicles can trigger mEPSCs. The amplitude of mEPSC is thought to reflect postsynaptic receptor density. Here, we assumed that synaptic plasticity affects the amplitude of mEPSC (AmEPSC) (Paré et al., 1997; Prange and Murphy, 1999):
where AmEPSC0 is an initial value of AmEPSC, and f = 0.01 is a factor representing the change of STDP on AmEPSC is slower than on gAMPA.
Data analysis
Analysis of synaptic weight.
The total input synaptic weight was measured as the sum of the strength of all incoming synaptic weights at each neuron, and the strength of recurrent connections was calculated by summing the product of the weights along each of the recurrent loops (incoming connection times outgoing connection). The spatial gradient of synaptic weight of each neuron was calculated by the difference of total input synaptic weight from the lower index neurons and the higher index neurons. The sign of spatial gradient of synaptic weight determines the direction of propagation. The positive value indicates that the wave prefers to propagate from the lower index neuron to the higher index neuron. The higher amplitude of spatial gradient of synaptic weight indicates higher preference to propagate in one direction.
Up-state initiation detection.
We first detected the timing of the onsets of Up state for each neuron, which is mainly based on bimodality of the membrane potential distributions (Volgushev et al., 2006) and the number of spikes during a short interval (e.g., 30 ms). We then selected Up states only if all the neurons are active simultaneously and their onsets should be separated by <300 ms, called clusters (Volgushev et al., 2006). We defined the neuron that actives first during an Up state as the initiation neuron.
Sequence detection and sequence replay measurement.
We used two methods to detect spike sequences. (1) We calculated the firing rate of each cell by the number of spikes during a moving window (with a window length 50 ms and window overlap 10 ms) divided by the duration of this interval and then smoothed it by a seven-point moving average. The firing sequences of neurons were determined by ordering the peaks of their smoothed firing rate during each Up state. (2) Because sequences are determined by the order of spike times, we computed a cross-correlation of Gaussian convoluted spike trains between the Up-state initiation neuron and all other neurons for each Up state. Based on the cross-correlation analysis, we identified time delays (lags of the cross-correlation peak) between neurons. This approach essentially gave us a “mean” or characteristic delay between spiking of any two neurons. Sorting these time delays provided the order of spiking of the neurons during one specific Up state (spike sequence).
We then used two methods to quantify sequence replay. (1) For a small subset neurons (e.g., four neurons), we calculated the probability of all possible sequences of these neurons before stimulation, during stimulation, and after stimulation. (2) For a large group of neurons, we applied a string distance (SD) method to quantify how similar two sequences are. The SD was calculated by SD (S1, S2) =
Statistical analysis.
The Kolmogorov–Smirnov test was used to determine whether the sample data have come from a normal distribution. The two-sample t test (for normal distribution) or the nonparametric Mann–Whitney U test (for non-normal distribution) was used for statistical analysis of differences between means from two samples when appropriate. One-way ANOVA was used for comparisons across the multiple groups, with Bonferroni's post hoc test applied to determine the significant group. Two-way ANOVA was used for comparisons across the different stimulation frequency groups with both periodic and Poisson stimulation (see Fig. 3C,D). For the histograms, a gray dotted line was plotted to represent the mean ± 2 STD (Standard Deviation). The significance of each bin was tested by z score, and an asterisk was put on the bin where its p value <0.01.
Computational methods
All model simulations were performed using a fourth-order Runge–Kutta integration method with a time step of 0.02 ms. Source C++ was compiled on a Linux server using the G++ compiler. All data processing was done with custom-written programs in MATLAB (MathWorks).
Results
Spontaneous slow oscillations in the thalamocortical network model
The network model included a three-layer cortical network (layers II/III/IV, V, and VI) with PY cells and INs, as well as a two-subsystem thalamic network (core and matrix; Jones, 2001; Bonjean et al., 2012) with TC and RE neurons (Fig. 1). The model spontaneously generated slow oscillations characterized by periodic transitions between active (Up) and silent (Down) states (Fig. 2A,B) with properties similar to recordings in vivo (Steriade et al., 2001). The Up states of the slow oscillation were initiated because of the summation of the spontaneous mEPSPs (minis) between cortical PY cells during Down states (Timofeev et al., 2000; Bazhenov et al., 2002). Although minis are rare events in each synapse (Stevens, 1993), occasional summation of minis across multiple synapses was sufficient to activate persistent Na+ currents and initiate spikes in one or few cortical cells, leading to Up-state initiation. Excitatory connections between cortical neurons (Bazhenov et al., 2002) and cortico-thalamo-cortical projections (Lemieux et al., 2014) mediated spread of activity over the entire cortical and corresponding thalamic networks. Active states were then terminated by the slow activation of the Ca2+-dependent K+ current, synaptic depression (Timofeev et al., 2000; Bazhenov et al., 2002), and the buildup of inhibition (see Materials and Methods; Chen et al., 2012). The network switched back to a silent state that was followed by initiation of the next Up state, leading to oscillations (Fig. 2A,B). We found that the global pattern of activity of PY cells was similar across cortical layers (Fig. 2A,B, top). The TC and RE cells also showed active and silent periods during corresponding Up and Down states in cortical neurons (Fig. 2A,B, middle and bottom). In TC cells, the active period usually consisted of a few spikes followed by subthreshold oscillations at 7–14 Hz and was commonly preceded by a brief hyperpolarization attributable to inhibition from RE cells in both the core and matrix subsystems. This spindle-like activity was similar to recordings in experiments with cats (Timofeev and Steriade, 1996).
Thalamocortical network model. The cortical column was modeled in three layers: layer II/III/IV, layer V, and layer VI. Each cortical layer was organized in a one-dimensional chain of PY cells and INs. The thalamus model consisted of core and matrix subsystems, including TC and RE neurons. The core TC neurons projected to the middle layers of the cortex (layer II/III/IV), whereas the matrix TC neurons projected to the superficial layer (layer I), making connections to the distal dendrites of layer V pyramidal neurons. The matrix thalamocortical pathway had diffuse connections to the cortex. Black filled circles and black bars represent excitatory and inhibitory connections between neurons, respectively.
Characteristic properties of slow oscillation in the model. A, Space–time plot representing characteristic activity of the network. The value of the membrane potential for each neuron is indicated with the color code. B, Voltage traces of individual neurons (PY in the layer II/III/IV and layer V, TC and RE in the core and matrix networks). C, Probability of Up-state initiation in three cortical layers. Layer V had higher probability to initiate an Up state. *p < 0.01, one-way ANOVA with Bonferroni's post hoc test. D, Effect of thalamocortical connectivity on slow oscillation. Down-state duration is shown in the model with different connection strength between thalamus and cortex. Increasing connection strength led to significant reduction of Down-state duration. *p < 0.01, one-way ANOVA with Bonferroni's post hoc test. Error bars in C and D indicate STD.
Cortical and thalamic networks contributed to specific features of the slow oscillations. Experimental studies (Sanchez-Vives and McCormick, 2000; Chauvette et al., 2010) revealed that Up states of the slow oscillations are more often initiated in layer V neurons. In our model, higher probability of Up-state initiation in layer V (Fig. 2C; p < 0.01, one-way ANOVA with Bonferroni's post hoc test) resulted from a higher number of synaptic inputs to the layer V pyramidal neurons, based on anatomical data of cortical connectivity (Thomson and Lamy, 2007). After initiation in layer V, slow oscillations propagated to other layers with a similar spatiotemporal pattern (Fig. 2A,B). Thalamic neurons are known to play a role in synchronization of active states of slow oscillation in vivo (Lemieux et al., 2014). In our model, the Down-state duration was significantly increased when the thalamus was removed (Fig. 2D, see “No conn.” bar, p < 0.01, one-way ANOVA with Bonferroni's post hoc test). We found that the cortical network alone generated slow oscillation activity that was similar to recordings in vitro in cortical slices and isolated cortical slabs (Sanchez-Vives and McCormick, 2000; Timofeev et al., 2000). Thus, thalamic network was necessary to obtain in vivo-like pattern of slow oscillation and to ensure that the synaptic reorganization described in this study occurs in the presence of the synchronized effect of matrix thalamic projections (Bonjean et al., 2012; Lemieux et al., 2014).
We next examined statistical features of the spatiotemporal pattern of the slow oscillations. In the model, spontaneous spike-independent minis mediated transitions from silence to active cortical states (Bazhenov et al., 2002; Chauvette et al., 2010), and thus any network site could potentially initiate an Up state (Fig. 3A, top and top middle). Initially, because of the same synaptic connectivity, the probability of Up-state initiation was nearly uniformly distributed across the network (Fig. 3A, bottom middle). Such distribution of Up-state initiation sites remained regardless of initial conditions (initial states of the random variables see Materials and Methods). Surprisingly, we also found that a network site that initiated an Up state at the previous cycle of oscillation was more likely to initiate it at the following cycle (Fig. 3A, bottom). Such history dependence resulted from the refractoriness of network. During an Up state, there was an accumulation of intracellular Ca2+, which led to a progressive increase in the hyperpolarizing Ca2+-dependent K+ current; IKCa then decreased during the Down state. This intrinsic current along with depression of synapses led to a refractory-like period affecting the generation of the next Up state. Because Up-state initiating neurons displayed slightly earlier termination compared with the rest of the network, these neurons were released from the refractory period earlier and thus had higher likelihood to generate the next Up state. Furthermore, the probability of synaptic minis accumulation increased with time after termination of the last Up state. This also facilitated initiation of the next Up state at the same location where Up state was initiated in the previous cycle. Our model suggests that intrinsic and synaptic properties determine minimal time of Down- to Up-state transition and ultimately minimal period of the slow oscillation. These results are consistent with in vivo in slab experiments, in which the period of slow oscillations depended on the cortical network size (Timofeev et al., 2000). In cortical slab with a reduced number of connections and limited size of the neuronal population, slow oscillation became irregular with an average period between Up states reaching up to 60 s (Timofeev et al., 2000). Interestingly, keeping isolated cortical slab in vivo for few days led (as we hypothesized) to homeostatic changes, and the period of slow oscillations reached again 1–3 s (Lemieux et al., 2014). Furthermore, it has been reported that Up states recurred with a periodicity of 2–5 s in vitro (Sanchez-Vives and McCormick, 2000).
External stimulation influences the spatiotemporal pattern of cortical slow oscillation. A, Representative examples of spontaneous network activity without external input. Top, Up-state initiation sites over the entire simulation time are indicated by black dots. The black triangular indicates the time interval expanded below. Top middle, Characteristic example of the network dynamics. Membrane voltage of pyramidal neurons in layer II/III/IV is indicated with a color code; white stars indicate site of Up-state initiation. Bottom middle, Probability of Up-state initiation across neurons. The dotted line represents mean ± 2 STDs. Bottom, Probability of Up-state initiation at the following cycle (i + 1) of oscillation as a function of the distance from the initiation site at the previous cycle (i). Because of the refractoriness properties of the network, probability of the next active state initiation was higher at the same network site that initiated activity at the previous cycle. B, Periodic (1 Hz) external stimulation (red arrow) was applied to pyramidal cells 101–103 in layer II/III/IV. Same panel description as in A. Note that Up-state initiation occurred mainly near the stimulation site (top). Spatiotemporal pattern of slow-wave propagation became organized (top middle). Probability of Up-state initiation around stimulation sites was increased significantly (bottom middle, *p < 0.01). Excess kurtosis of the Up-state initiation distribution (C) and probability of initiation at the simulation sites (D) averaged over five different random network initializations as a function of stimulation frequency for Poisson (blue) and periodic (black) stimulation. Two-way ANOVA test showed that there was a significant difference among stimulation frequencies (p=3.87 × 10−12 and p=1.09 × 10−24 for C and D, respectively) and no significant difference between Poisson and periodic distributed inputs (p = 0.8909 and 0.5733 for C and D, respectively). Error bars in C and D indicate STD.
The history effect (mediated by hysteresis) alone could only maintain repeatable pattern of Up-state initiation for up to two or three cycles, so in a long timescale, the Up-state initiations were uniformly distributed across the network without plasticity or external input (Fig. 3A, bottom middle). Nevertheless, this effect helped to maintain the stereotyped Up-state initiation pattern in the presence of even very sparse (in time) hippocampal input because it was necessary for reliable synaptic changes (see below).
External input influences the spatiotemporal pattern of slow oscillations
It has been proposed that declarative memory is consolidated by reactivation of neural representations of recent memories in the hippocampus and the neocortex during SWS, leading to formation of persistent cortical memory traces (Marshall and Born, 2007; Born, 2010). The key element of consolidation stage is cortical replay triggered by hippocampal sharp wave-ripple (SPW-R) events (Peyrache et al., 2009). Experimental studies have demonstrated that hippocampal neurons in CA1 project to multiple regions in the medial prefrontal cortex (Jay and Witter, 1991; Hoover and Vertes, 2007; Parent et al., 2010). In our study, we mimicked the hippocampal input to the neocortex with an external stimulation to a small group of pyramidal cells. Each input was a current step with 100 ms duration, which is similar to characteristic ripple duration in vivo (50–150 ms; Sullivan et al., 2011). We tested a broad range of inter-ripple intervals as found in experiments (Axmacher et al., 2008). Figure 3B shows the spatiotemporal pattern of cortical slow oscillation with a periodic (1 Hz) external input applied to a group of three pyramidal cells (cells 101–103) in layer II/III/IV of the cortical network. During stimulation, the pattern of Up-state propagation became organized (Fig. 3B, top middle), and the majority of active states were initiated around the stimulation site (Fig. 3B). This result is consistent with the pervious findings with a one-layer cortical model (Compte et al., 2003).
The probability of initiating an Up state at the simulation location was highest when the frequency of stimuli was ∼1 Hz, which was also the natural frequency of slow oscillation in our model (Fig. 3C,D). Lower stimulation frequencies led to less reliable initiation. However, the probability of Up-state initiation at the stimulation site remained above chance even for very low input frequencies (Fig. 3D). The excess kurtosis of the Up-state probability distribution peaked close to 0 for 1 Hz stimulation, indicating a close to normal distribution (Fig. 3C). A two-way ANOVA test showed that both periodic and Poisson distributed inputs produced similar results (p = 0.8909 and 0.5733 for Fig. 3C,D, respectively, compare blue and black lines), but there was a significant difference from the contribution of stimulation frequency (p < 0.01 for both Fig. 3C,D). These experiments suggest that even sparse hippocampal input could potentially influence the spatiotemporal pattern of slow oscillations in the neocortex.
STDP leads to the appearance of local “groups” of neurons initiating Up states
To explore effects of synaptic plasticity during slow oscillations, we first considered the condition when there was no external input to the network, and synaptic weights of the excitatory connections between all cortical neurons were allowed to change based on STDP rule (see Eq. 3). The synaptic changes were set to be very small: spiking activity during a single Up state affected synaptic weights by <1% of the initial value. In the presence of STDP, the main spatiotemporal properties of the slow oscillations, including duration and frequency of Up states, remained intact. However, because of the change of synaptic weights during spontaneous slow oscillations, the probability distribution of initiation sites deviated from the uniform distribution (Fig. 4, compared A,B). Some sites of the network (such as a group of neurons around sites 50 and 150 in the network shown in Fig. 4B) initiated Up states more often than other sites. In the model without external stimulation, such preferential initiation sites spontaneously emerged, and their location varied across simulation trials with different initial conditions.
Spontaneous slow oscillation in the presence of STDP. A, B, Representative examples of the network activity during early (A) and late (B) stages of simulation. Top, Up-state initiation sites are indicated by black dots. The black triangular indicates the time interval expanded below. Middle, Characteristic example of the network dynamics. Membrane voltage of pyramidal neurons in layer II/III/IV is indicated with a color code; white stars indicate the site of Up-state initiation. Bottom, Probability of Up-state initiation across neurons. The dotted line represents mean ± 2 STDs. The red curve represents the smoothed histogram by Gaussian filter. Synaptic reorganization leads to appearance of the local groups of neurons preferentially initiating Up states (near cells 50 and 150) during the late stage of simulation (*p < 0.01). C, The percentage change of synaptic weights at time t = 1000 s compared with the initial (t = 0) synaptic weights shows increase (red) or decrease (blue) of synaptic weights. The x-axis gives the index of the postsynaptic neuron; the y-axis represents the left (−) and right (+) 10 presynaptic neurons relative to the postsynaptic neurons at x-axis. D, At t = 5000 s, total input synaptic weight received by each neuron (top) and recurrent input strength and spatial gradient of synaptic weights (bottom). Note the group of neurons in the middle (neurons 80–120) with strongly asymmetric connectivity matrix (with connections from one direction of the network enhanced while the connections from the opposite direction reduced) that had high absolute amplitude of spatial gradient would have low probability to initiate the Up states. E, Time evolution of the total (net) input synaptic weight per neuron for 200 pyramidal neurons from layer II/III/IV. F, Individual synaptic weights of two neurons (indicated by stars in D) with low (left) and high (right) net synaptic weight. G, Time evolution of the size (neuron count) of the cell population that had strongly asymmetric connectivity matrix.
The emergence of the repetitive spatiotemporal pattern of slow oscillation was determined by synaptic changes. Indeed, in the model with STDP, the direction of Up-state propagation determined whether the synaptic connection increased or reduced after each wave (Up state) propagation. On average, connections increased from any neuron that had Up state earlier in time to any neuron that had it later in time, whereas the opposite connections decreased, thus leading to strongly asymmetric weight distribution (neuron 100 in Fig. 4C). Because of the refractoriness of the network, this effect could accumulate across cycles of slow oscillation mediating systematic changes of synaptic weights.
To reveal how distribution of synaptic weights influences initiation of Up states, we examined whether total synaptic weight per neuron, strength of recurrent excitatory connections, or spatial gradient of synaptic weights along the network (see Materials and Methods) could predict the site of Up-state initiation. Up states were initiated more often within groups of neurons with variable total synaptic strength (Fig. 4D, top). Among those neurons, the neurons with maximal total synaptic input had many of the incoming synapses reaching the upper limit (neuron 151 in Fig. 4F; brown star in Fig. 4D), and the neurons with minimal total synaptic input had many of the incoming synapses reaching the lower limit (neuron 154 in Fig. 4F; blue star in Fig. 4D). The neurons with maximal total synaptic input experienced high synaptic (minis) drive during Down states and therefore could more likely initiate a spike. Furthermore, because not all spikes could lead to an Up state, Up states were initiated mainly at the network sites in which recurrent excitatory connection strength between pyramidal cells was high (Fig. 4D, bottom, blue line); this increased the probability that an initial spike would trigger a wave of excitation. Finally, close to zero spatial gradient of the total synaptic weight per cell indicated locations of bidirectional wave propagations that was also a predictor of the initiation site (Fig. 4D, bottom, green line). [The sign of the spatial gradient at a specific location determined the preferential direction of wave propagation. For example, for neurons 80–120, the spatial gradient was positive (Fig. 4D, bottom, green), indicating that the wave propagation was primarily unidirectional and from the lower index neuron to the higher index neurons.] For the sample network shown in Figure 4D, we predicted that the group of neurons (neurons 50–70) with highly variable synaptic weights, high recurrent synaptic strength, and close to zero spatial gradient of synaptic weight was the main initiation site and another group (neurons 160–170) formed the secondary site. This was indeed consistent with the distribution of the probability of Up-state initiation shown in Figure 4B.
The spatiotemporal pattern of slow oscillation, once it was formed, remained throughout the simulation time (Fig. 4E). Why did this pattern preserve? We found that locations of Up-state initiation determined the dynamics of synaptic weights and vice versa. This mutual interaction helped preserve the pattern. Over time, the overall size of the regions of unidirectional wave propagation slowly increased (Fig. 4G), forming a spatiotemporal pattern with Up-state initiation occurring primarily in a few specific locations. This may suggest that this spatial pattern of the synaptic weights is an attractor state for the network.
External stimulation in the presence of STDP leads to input-specific synaptic weight reorganization
Next we applied a periodic (1 Hz) stimulation to few cortical pyramidal neurons (neurons 101–103) in layer II/III/IV. As before, the properties of the stimulation were designed to represent minimal characteristics of the hippocampal ripples. We observed that the probability of Up-state initiation at the stimulation sites increased (Fig. 5B). As we showed in the previous section, interaction of the slow waves and STDP led to characteristic synaptic changes; however, in presence of stimulation, these changes reflected input pattern to the network. Groups of cells with higher probability of Up-state initiation have been formed, and their location was not random but reflected location of the stimulation sites. This pattern remained even after stimulation was removed (Fig. 5, compare A,C). Importantly, this effect was independent on the initial state of the network (initial conditions of all the random variables; see Materials and Methods).
Pattern of slow oscillation in the presence of STDP and external stimulation. A–C, Representative examples of network activity before stimulation (A), during stimulation (B), and after stimulation (C). The site of stimulation is indicated with red arrow. Top, Up-state initiation sites are indicated by black dots. Middle, Characteristic example of the network dynamics. Membrane voltage of pyramidal neurons in layer II/III/IV is indicated with a color code; white stars indicate site of Up-state initiation. Bottom, Probability of Up-state initiation across neurons. The dotted line represents mean ± 2 STDs. The red curve represents the smoothed histogram by Gaussian filter. Note that the probability of Up-state initiation was very high around the stimulation site (cell 100) during stimulation (B, *p < 0.01), and it remained high after stimulation was terminated (C, *p < 0.01). D, The percentage change of synaptic weights at time t = 300 s (left), t = 1000 s (middle), and t = 3500 s (right) compared with the initial (t = 0) synaptic weights. The x-axis gives the index of the postsynaptic neuron, and the y-axis represents the left (−) and right (+) 10 presynaptic neurons relative to the postsynaptic neurons at x-axis. E, Total input synaptic weight received by each neuron at the end of stimulation (t = 3500 s). F, Recurrent input strength and spatial gradient of synaptic weight at t = 3500 s.
In these experiments, external stimulation was applied at the early stage of the network, so before the stimulation, synaptic weights were distributed almost uniformly across the network (Fig. 5D, left). During stimulation, Up states initiated mainly near the input sites (Fig. 5B). This led to stereotyped wave propagation and synaptic changes reflecting the propagation pattern as we described previously (Fig. 4). Synaptic connectivity matrix of many neurons became asymmetric, with connections from one side of the network enhanced whereas the connections from the other side reduced, except around stimulation sites. Such an asymmetric synaptic pattern already emerged during the early stage of stimulation (Fig. 5D, middle) and persisted after stimulation was removed (Fig. 5D, right).
Why was Up-state initiation more likely to occur around stimulation sites even after termination of the input? To answer this question, we analyzed distribution of the total synaptic weight per neuron, strength of recurrent excitatory connectivity, and spatial gradient of synaptic weights along the network. Although relatively large cell populations had highly variable synaptic weights (Fig. 5E), recurrent connectivity and spatial gradient of synaptic weights predicted initiation near site 100, where external input was applied during the stimulation phase (Fig. 5F).
We interpreted these results as formation of long-term memory representation in the cortical network after period of replay driven by hippocampal ripples. This effect can be generalized to more complex input patterns with multiple stimulation sites (Fig. 6A–C). We found that both stimulation sites displayed higher probability to initiate Up states during stimulation (Fig. 6B). After input termination, there was a relatively high probability of Up-state initiation near both stimulation sites, although one site (neurons 75–77) dominated Up-state initiation (Fig. 6C, bottom).
Model of the hippocampal input applied at two cortical network sites. A–C, Two inputs are applied simultaneously (at cells 75–77 and 125–127). Up-state initiation sites (top), representative examples of the network activity (middle), and distribution of the probability of initiation (bottom) before stimulation (A), during stimulation (B), and after stimulation (C). The dotted line represents mean ± 2 STDs, whereas *p < 0.01 from z score. Red arrows shows the sites of two simultaneous inputs. D, E, Two inputs are applied sequentially (first at cells 75–77 and then at cells 125–127). D, Up-state initiation sites (black dots). Green bar, Duration of the first stimulus; red bar, duration of second stimulus. E, Distribution of the probability of Up-state initiation after stimulation (*p < 0.01). The dotted line represents mean ± 2 STDs. The red curve represents the smoothed histogram by Gaussian filter.
We found that preexisting memory traces remained in the cortical network after new stimulation was applied (Fig. 6D,E). Once the first memory was imprinted into the network connectivity during the first (earlier) stimulation period, it remained strongly expressed after the second (later) stimulation (Fig. 6E). The second memory trace integrated with the first one, forming two network sites with a higher than average probability of Up-state initiation (Fig. 6E). Because of the small network size, only local intracortical synaptic connectivity and reduced dimensionality (one-dimensional network) of our model, these two memories shared some overlapping fragments of the spike sequences. In a large-scale multidimensional system, better approximating complexity of the animal and human brain, we expect that different and non-overlapping sequences would represent different memories. Our study suggests that the mechanisms of interaction between hippocampal inputs and cortical slow oscillation can maintain storage of multiple memories. Mapping of the new memories does not necessarily disrupt the previously stored memories. However, we also found that the exact outcome of this experiment depends on the strength of the earlier memory, which in our model depended on the duration of the first stimulation. We found that, if the first stimulation was long enough leading to sufficiently strong synaptic changes, then during stimulation at the second site, both network sites continued to initiate Up states simultaneously, resulting in a complex pattern with the interference between two stimulation sites (data not shown). Thus, our results suggest that the initial period of SWS is crucial for memory consolidation, consistent with experimental studies (Gais et al., 2000; Eschenko et al., 2008).
External stimulation in the presence of STDP leads to stimulation-specific cortical replay
In the previous section, we found that, after stimulation, cortical neurons at the stimulation site had higher probability to initiate an Up state (Fig. 5B,C). Here, we showed that such preferential Up-state initiation results in repetition of the firing sequences of neurons. First, we showed that, within an Up state, the relative timing of the neuronal spiking is preserved on average and similar to that at the Up-state onset. Using cross-correlation analysis, we found “mean” or characteristic delay between spiking of any two neurons over the entire length of any given Up state. We then compared these characteristic delays with those at the Up-state onset and found that the order of cell spiking is primarily preserved from an Up-state initiation to its entire duration. To illustrate this result, we showed a typical Up-state spike pattern in Figure 7A. Spikes in different neurons formed “stripes” during an Up state following the same pattern as that at the onset of the Up state. Characteristic delay between an Up-state initiating neuron and any other neuron increased with distance from the Up-state initiating neuron, reflecting wave propagation pattern (data not shown).
Spike sequence replay. A, Spike raster plot of a typical Up state. B, Distribution of the SDs between every sequence and all the other sequences from the subset of Up states initiated at the same network location, repeated independently for each network location (blue), and distribution of SDs between every sequence and all the other sequences with randomly shuffled spike times from the same subset of Up states, repeated independently for each network location (red). There was a significant difference between the means of the two distributions (*p < 0.01, Mann–Whitney U test). C, SD between template sequence (sequence 101, 102, 100, 103, 99, 104, 98 … for initiation site 101) and spike sequences from all the Up states as a function of time. Stimulation was applied at sites 101–103; red bar indicates stimulation time. Blue dotted line represents mean SD between the template sequence and sequences extracted from Up states before stimulation. D, Mean SD between a template spike sequence and all the spike sequences before, during (from 300 to 3500 s, red bar), and after stimulation, as shown in C. Mean SDs during and after stimulation were significantly different from that before stimulation (*p < 0.01, one-way ANOVA with Bonferroni's post hoc test). Error bars indicate standard error. E, Representative examples of four neurons (a = 20, b = 35, c = 50, d = 65) spike trains (top), smoothed firing rate (middle), and probability of each sequence (bottom) before, during, and after stimulation. The dotted line represents mean ± 2 STDs. Note the high probability of the “dcba” sequence during (*p < 0.01) and after (*p < 0.01) stimulation.
Next, we demonstrated that all the Up states initiated from the same neuron result in similar spike sequences. We used string distance (SD) to measure similarity between sequences (for details, see Materials and Methods). For each network site (neuron), we detected spike sequences from the Up states initiated at that particular network location. Within these sequences, we then calculated SD between every sequence and all the other sequences (so, essentially, between all the pairs of sequences). This was repeated independently for each Up-state initiation site. Distribution of all the SDs calculated this way (Fig. 7B, blue) was significantly lower than distribution of the SDs between every sequence and all the other sequences of the subset with randomly shuffled spike times, also repeated for each initiation site (Fig. 7B, red). This suggests that all Up states initiated by the same neuron result in similar spike sequences.
We then compared SDs before, during, and after the stimulation (for details, see Materials and Methods). We first artificially generated a template sequence representing an “ideal” sequence initiated at a specific network site (for example, sequences 101, 102, 100, 103, 99, 104, 98 … for initiation site at 101). We then computed the SD between this template sequence and every other sequence generated by the network as a function of time (Fig. 7C). The mean value of SD before stimulation represented baseline (Fig. 7C, blue dotted line). Because of the random location of sites of Up-state initiation before stimulation, the SD was significantly higher (Fig. 7D, p < 0.01) and showed large fluctuations around the mean (Fig. 7C). The SD was significantly reduced during stimulation and displayed only minimal fluctuations (Fig. 7C,D). Importantly, we found that the SD remained relatively lower after stimulation was removed (Fig. 7D), which corroborates our prediction that the change of synaptic weights during stimulation promotes specific sequence replay after stimulation.
Finally, we examined sequence replay by a subset of neurons in which sequences were detected based on the firing rate. We selected four cortical neurons and assigned them letters (a = 20, b = 35, c = 50, d = 65) to test repetition of specific letter combinations, similar to the previous studies (Ji and Wilson, 2007). Figure 7E (top) illustrates the characteristic spiking activity of four selected neurons during slow oscillations. The firing sequences of these neurons were determined by ordering the peaks of their smoothed firing rate (Fig. 7E, middle; Lee and Wilson, 2002). We then computed the probability of each sequence that could be formed by four selected neurons (Fig. 7E, bottom) before, during, and after the stimulation. Before stimulation, initiation of Up states was random, which was reflected in the low probability of any specific sequence, e.g., “dcba” sequence (Fig. 7E, left column). During stimulation, the pattern of slow oscillations reflected Up-state propagation from the stimulation site. The neurons fired sequentially in the direction of Up-state propagation; therefore, the probability of the sequence “dcba” became high (Fig. 7E, middle column). When stimulation persisted long enough, we observed that the probability of “dcba” sequence remained high even after stimulation was removed (Fig. 7E, right column). This suggests that external input reorganized synaptic weights to increase the probability of specific firing sequences. When stimulation was applied at the very different network sites (e.g., around site 150), the network formed different dominant sequences (“abcd” and “abdc”) that also remained after the stimulus was removed.
In all of the above simulations, we assumed that memory traces are initially stored in the hippocampus, and we presented evidence that the input from the hippocampus to the cortical network during slow oscillations may govern formation of corresponding memory traces in the cortical network. We also tested the possibility that memory traces are embedded to the cortical network during the state preceding sleep (e.g., wakefulness), and we asked the question how these traces would be modified by slow oscillations. We found that spike sequences of cortical neurons reflecting those embedded synaptic patterns replayed during ongoing slow oscillations and that led to further enhancement of characteristic synaptic patterns (data not shown).
Global (such as transcranial) stimulation augments the effect of synaptic changes associated with local (such as hippocampal) input
In our previous simulations, the local input targeted only a few cortical neurons. This type of input was designed to model the interaction between the hippocampus and neocortex during sleep. In the study by Marshall et al. (2006), transcranial electrical stimulation led to global modulation of the neural membrane potentials to boost slow oscillations at the early stage of NREM sleep. To test the effect of such stimulation in our model, we applied global input (with a frequency of 1 Hz) to the entire network to modulate slow oscillation. At the same time, we applied a local stimulation, as described in the previous sections, to the neurons 101, 102, and 103 to mimic input from the hippocampus. We found that initiation of the Up states was more likely to occur near the stimulation site when both local and global stimuli were applied simultaneously (Fig. 8B) compared with the local stimulation alone (Fig. 8A). Initiation of the Up states became more synchronized because of the global input (Fig. 8B, middle top). As a result, synaptic weights increased more rapidly (Fig. 8, bottom), and after the stimulation, the probability of Up-state initiation was higher around the stimulation site (Fig. 8, middle bottom) when the global stimulation was coupled with the local input. These findings suggest that global external stimulation, when co-occurring with the local input (such as normal projections from hippocampus to the cortex during NREM sleep), can accelerate synaptic changes associated with sleep-related memory consolidation, consistent with experimental results obtained previously (Marshall et al., 2006). In the “closed-loop” hippocampo-cortical network in vivo, the timing of the hippocampal SPW-Rs is influenced by the cortical activity, and those hippocampal events are more likely to occur during transition from a Down to Up state of slow oscillations (Isomura et al., 2006; Mölle et al., 2006). Our study suggests that transcranial electrical stimulation can enhance synchrony of Up states and promote sequence activation in the hippocampo-cortex loop.
Effect of global network stimulation. Local input (at neurons 101, 102, 103) was applied without (A) and simultaneously with (B) global modulation of the network activity. Top, Up-state initiation sites (black dots). Red bar represents the time of stimulation. Top middle, Representative example of spatiotemporal patterns of Up and Down states. Note the higher synchrony of slow oscillation in the presence of global stimulation (right). Middle, Distribution of Up-state initiation sites during the entire 1000 s simulation time. Bottom middle, Distribution of Up-state initiation sites after stimulation. Bottom, The percentage change of synaptic weights at the end of stimulation (t = 500 s) compared with initial (t = 0) synaptic weights. The x-axis gives the index of the postsynaptic neuron, and the y-axis represents the left (−) and right (+) presynaptic neurons relative to the postsynaptic neurons at the x-axis. The dotted line represents mean ± 2 STDs, whereas *p < 0.01 from z score. The red curve represents the smoothed histogram by Gaussian filter.
Discussion
In this study, we showed that the spatiotemporal pattern of sleep slow oscillations in the neocortex could be influenced by the sparse synaptic input from the hippocampus. In the presence of synaptic plasticity, this leads to the input-specific changes in synaptic connections between cortical neurons and formation of a persistent locus of preferred Up-state initiation at specific cortical sites. This pattern remains even without hippocampal input, which manifests consolidation of the memory traces in the neocortex by hippocampal input. Our study thus suggests a putative neuronal mechanism for memory consolidation during SWS.
Interaction between slow oscillations and synaptic plasticity
Synaptic plasticity is believed to be a cellular mechanism of learning and information storage in the brain. STDP is an experimentally well-characterized form of plasticity that adjusts synaptic weights based on the relative timing of the presynaptic and postsynaptic spikes (Song et al., 2000). Our modeling study supports the idea that thalamocortical activity during SWS can be well suited to promote specific and precise synaptic changes by STDP.
We observed mutual interaction between synaptic changes and spatiotemporal pattern of slow oscillations. First, the pattern of slow oscillations determined the change of synaptic connections between cortical neurons. During slow oscillations, the transition from a Down to Up state of oscillations occurs in a form of traveling waves (Massimini et al., 2004; Volgushev et al., 2006), possibly initiated at many independent network locations. In the presence of STDP, synaptic connections between neurons changed according to the relative order of spiking: synaptic connections in the direction of wave propagation potentiated, and those in the opposite direction depressed. Thus, the repetitive pattern of slow oscillation led to characteristic changes of the network connectivity.
Second, the changes of synaptic connections can in return affect pattern of slow oscillations. We found that, after a long enough period of slow oscillations, the synaptic connectivity matrix of many neurons became strongly asymmetric (inputs from one direction were higher than those from the opposite direction). These neurons also had a low probability of excitatory recurrent connections and were less likely to initiate an Up state. In contrast, synaptic weights of a few neurons located near stimulation sites exhibited a complex pattern with either high or low total (net) synaptic weight. The neurons receiving high total synaptic input also had high recurrent connectivity and close to zero spatial gradient of synaptic weights. These neurons had the highest probability for initiating Up states. Our results are consistent with in vivo finding in cats indicating that the Up states preferentially spread from one or few local foci during slow oscillations (Volgushev et al., 2006). Previous computational studies also reported that, when STDP was included to the cortical model of Up and Down states, the network self-organized; strengthening of connections led to reactivation of precise temporal sequences (Kang et al., 2008). This study predicts that, within a cortical architecture, the interaction between slow oscillations and synaptic plasticity may form local clusters of neurons involved to Up-state initiation, and location of those clusters is altered by input from other brain areas. These complex synaptic associations between neurons could potentially contribute to the formation of assembly sequences for reactivation (O'Neill et al., 2010).
In the model, the global synaptic dynamics leads to bimodal distribution in the long term as reported in many studies of STDP. We implemented a hard boundary (Gerstner et al., 1996; Kempter et al., 1999; Song et al., 2000; Drew and Abbott, 2006) in this study to prevent STDP from runaway synaptic dynamics. Such synaptic dynamics can be potentially prevented by including heterosynaptic plasticity (Chen et al., 2013). However, the characteristic changes in synaptic weights and spatiotemporal dynamics of slow oscillations emerged early in our simulations and were qualitatively similar for the rest of the simulation. This suggests that consolidation of memory and replay of learned sequences does not require extreme changes in synaptic weights.
Although previous studies have explored the role of STDP during sleep (Esser et al., 2007; Olcese et al., 2010), these works mainly focused on the global synaptic weight dynamics (synaptic homeostasis during sleep). In contrast, our study is focused on the local synaptic reorganization mediated by the spatiotemporal pattern of sleep slow oscillations. We found that sleep slow-wave activity leads to the change of the connectivity between cortical neurons, and, in return, changes of the synaptic weights affect the pattern of slow oscillations. Ultimately, our study aims to explain cortical sequence generation that requires a fine-tuning of specific connections rather than global scaling of synaptic weights. The novelty of our research is in linking local synaptic dynamics to specific patterns of sleep activity and ultimately to sequence replay.
Interaction between the neocortex and hippocampus to consolidate memory
It has been proposed that declarative memory is consolidated by interaction between the hippocampus and neocortex during sleep, leading to long-term memory storage in the neocortex (Marshall and Born, 2007; Born, 2010; Schwindel and McNaughton, 2011). Hippocampal SPW-R events have been hypothesized to provide active communication between the hippocampus and neocortex during SWS (Buzsáki, 1996). Indeed, the low cholinergic activation, such as during SWS, permits information flow from the hippocampus to the neocortex (Hasselmo, 1999). Several studies supported the hypothesis that SPW-R events occur mainly in the neocortical transition from Down to Up states (Battaglia et al., 2004), and suppression of hippocampal ripple results in poor performance in memory tasks (Ego-Stengel and Wilson, 2010). In our study, we have mimicked hippocampal SPW-R input to the neocortex by sparse external cortical stimulation. Duration of the stimuli events matched the typical length of SPW-R events (Sullivan et al., 2011). We have found that such input promotes Up-state initiation at specific network locations and induces long-term changes in synaptic conductances, leading to the sequence replay. These results implicated the importance of the hippocampal SPW-R events in transferring memory information to the neocortex.
There are evidences suggesting that spindles occur during the Up states of slow oscillations (Mölle et al., 2011), and spindles are temporally aligned with hippocampal ripples (Siapas and Wilson, 1998; Sirota et al., 2003; Staresina et al., 2015). Long-term synaptic changes induced by spindles were studied in rat somatosensory layer V pyramidal cells, suggesting that spindles can support long-term storage of memory-related information in the neocortical networks (Rosanova and Ulrich, 2005). These experimental studies imply that precise coordination of the ripples, spindles, and slow oscillations plays an important role in memory consolidation.
Replay of sequences arises from reactivation of spatiotemporal pattern of slow oscillations
During SWS, the brain reorganizes its synaptic connectivity to store new information acquired during awake. Some studies have suggested that the primary effect of sleep is homeostatic downregulation of synaptic efficacy (Olcese et al., 2010). Discovery of replay, a phenomenon when temporally ordered firing sequences occurred during recent awake are replayed during sleep-promoting memory consolidation (Ji and Wilson, 2007; Mehta, 2007), suggested a more complex picture of synaptic changes during sleep. Although some synaptic weights may indeed decay, others, corresponding to recent awake experience, would enhance to form long-term memories. In our model, we observed sequential firing of neurons during slow oscillations, triggered by interaction between the input and the intrinsic dynamics of the thalamocortical network. Although input to the model was limited by short transient stimulation at one or few cortical sites, our study predicts that more complex stimulation patterns, such as representing sequences of firing of hippocampal neurons, would trigger characteristic responses (firing sequences) in the neocortical network. We found that reorganization of cortical connectivity attributable to STDP increased the likelihood of specific sequence replay even after hippocampal input was removed. This suggests a mechanistic explanation for the consolidation of specific memories during SWS, whereby the memory traces become independent of the hippocampus and are stored in the neocortex.
Model limitations and future directions
Our simplified one-dimensional thalamocortical model revealed the basic principles of how interaction between hippocampal SPW-Rs, cortical slow oscillations, and synaptic plasticity during sleep may lead to memory-specific synaptic changes, which in return promote reliable spike sequences. In our model, (1) hippocampal input was limited by one or two cortical sites and uniquely determined pattern of Up-state propagation, (2) the spatiotemporal pattern of the Up state then determined the relative spike timing of neuronal firing within an Up state–cortical spike sequence, and (3) two nearby cortical initiation sites triggered overlapping but non-identical spike sequences. In a biological system, many additional factors may affect spike sequences during memory replay, including the following: (1) memory-specific changes of the intracortical connectivity during training phase in awake; (2) the presence of the plastic long-range intracortical connectivity; and (3) existence of the multiple cortical input sites for a SPWR event, leading to a complex spatiotemporal pattern of slow oscillation. In such complex settings, the relationship between hippocampal memory events and cortical spike sequences can be nonlinear with abrupt change of spike sequences as a function of changing hippocampo-cortical input.
Summary
The results of our study present, for the first time, a mechanistic explanation of how characteristic brain activity during stage 3/4 sleep may promote consolidation of recent declarative memory events. It makes predictions that can be tested experimentally, including specific interventions to suppress or augment memory consolidation processes.
Footnotes
This work was supported by Office of Naval Research/Multidisciplinary University Research Initiative Grant N000141310672 and National Institutes of Health Grants R01 EB009282 and R01 MH099645.
The authors declare no competing financial interests.
- Correspondence should be addressed to Maxim Bazhenov, Department of Cell Biology and Neuroscience, University of California at Riverside, Riverside, CA 92521. bazhenov{at}salk.edu
