Selection of stimulus parameters for enhancing slow wave sleep events with a Neural-field theory thalamocortical computational model

Slow-wave sleep cortical brain activity, conformed by slow-oscillations and sleep spindles, plays a key role in memory consolidation. The increase of the power of the slow-wave events, obtained by auditory sensory stimulation, positively correlates to memory consolidation performance. However, little is known about the experimental protocol maximizing this effect, which could be induced by the power of slow-oscillation, the number of sleep spindles, or the timing of both events’ co-occurrence. Using a mean-field model of thalamocortical activity, we studied the effect of several stimulation protocols, varying the pulse shape, duration, amplitude, and frequency, as well as a target-phase using a closed-loop approach. We evaluated the effect of these parameters on slow-oscillations (SO) and sleep-spindles (SP), considering: (i) the power at the frequency bands of interest, (ii) the number of SO and SP, (iii) co-occurrences between SO and SP, and (iv) synchronization of SP with the up-peak of the SO. The first three targets are maximized using a decreasing ramp pulse with a pulse duration of 50 ms. Also, we observed a reduction in the number of SO when increasing the stimulus energy by rising its amplitude. To assess the target-phase parameter, we applied closed-loop stimulation at 0º, 45º, and 90º of the phase of the narrow-band filtered ongoing activity, at 0.85 Hz as central frequency. The 0º stimulation produces better results in the power and number of SO and SP than the rhythmic or aleatory stimulation. On the other hand, stimulating at 45º or 90º change the timing distribution of spindles centers but with fewer co-occurrences than rhythmic and 0º phase. Finally, we propose the application of closed-loop stimulation at the rising zero-cross point using pulses with a decreasing ramp shape and 50 ms of duration for future experimental work. Author summary During the non-REM (NREM) phase of sleep, events that are known as slow oscillations (SO) and spindles (SP) can be detected by EEG. These events have been associated with the consolidation of declarative memories and learning. Thus, there is an ongoing interest in promoting them during sleep by non-invasive manipulations such as sensory stimulation. In this paper, we used a computational model of brain activity that generates SO and SP, to investigate which type of sensory stimulus –shape, amplitude, duration, periodicity– would be optimal for increasing the events’ frequency and their co-occurrence. We found that a decreasing ramp of 50 ms duration is the most effective. The effectiveness increases when the stimulus pulse is delivered in a closed-loop configuration triggering the pulse at a target phase of the ongoing SO activity. A desirable secondary effect is to promote SPs at the rising phase of the SO oscillation.

Humans spend about one-third time of their life sleeping. This behavior has paramount 2 importance for the process of learning, as it contributes to the consolidation of 3 memories [1,2]. During the non-rapid eye movement (NREM) phase of sleep, the brain's 4 electrical activity is characterized by the occurrence of events recognizable on the 5 cortical activity measured by electroencephalogram (EEG) registers. These events are 6 the slow oscillations (SO), single high-amplitude cortical oscillations lasting 0.8 to 2 activity [20,21]. In this way, mechanistic insights about the generation of SP and SO 48 activity can guide the search for a better stimulation that maximizes its efficacy. In this 49 work, we present a systematic study of possible stimulation protocols and the more 50 efficient features of the stimuli in promoting these SWS events' co-occurrence. 51 Specifically, in this work we use a model based on the Neural Field Theory 52 (NFT) [22,23] with corticothalamic loops between two cortical neural populations: 53 excitatory or pyramidal(e) and inhibitory or interneurons(i ), and two thalamic 54 populations: reticular nucleus(r ) and relay nuclei(s). This model is appropriate for the 55 proposed exploration as it simulates the cerebral structures that originate SO and SP, 56 namely the cortico-thalamic loops. Moreover, this model has been previously tuned to 57 reproduce the electrical brain activity in different states: from arousal to slow-wave 58 sleep and including the different stages of NREM sleep [24,25]. In this NFT model, 59 sensory stimulation can be introduced as spatial-localized perturbations on the noisy 60 input neural population, obtaining an EEG-like signal as the model's output. From this 61 signal, we can assess the efficacy of the stimulation paradigms and implement a 62 closed-loop feedback mechanism that can fine-tune the stimuli delivery timing. 63 With the purpose of simulating a single-modality sensory stimulation, we searched 64 Furthermore, Fig 1A shows excerpts of the cortical excitatory population activity for 97 N2 and N3 stages, and we traced a trajectory between these two stages using linear 98 interpolation on the connection strengths. Other examples of activity trace are shown 99 within the linear interpolation and next to variations of the steady-state input value. 100 Dominant spindle activity is appreciable using the N2-stage parameters indicated 101 in [24]. The spindle activity decreases as the simulation point slides along the 102 interpolated trajectory towards the N3 stage, which shows slow-wave sleep dynamics 103 but not noticeable spindle activity. 104 We are interested in an intermediate state between N2 and N3. For this, we 105 concentrated on having energy both in N2 and N3 frequency bands. We obtained this 106 condition in the second half of the linear interpolation between N2 and N3. Within this 107 segment, we arbitrarily selected the point at 2/3 on the N2 to N3 trajectory. Fig 1B   108 shows the linearized model spectrums of N2 and N3 stages and three extra points in the 109 interpolation trajectory, including the 2/3 position in the blue-line. 110 The increase of the steady-state input value (φ (0) n ) raises the model response in the 111 Z-axis (Fig 1A), indicating lightness of NREM sleep [25,27]. Besides, it also increases 112 the power and frequency of the SP-band peak (Fig 1C), then a higher steady-state input 113 value facilitates SP occurrence. 114 We selected φ Parameters of the model used in all the simulations.
The variation of the input, expected with the addition of stimuli signal, has a high 121 effect in the relay nuclei activity φ   Using the time-averaged scalograms, we calculated the spectrums shown in Fig 2D. 145 We define the power difference index, I (SO,SP ) , as the area between the corresponding 146 STIM condition and SHAM spectrums inside the frequency band of interest, normalized 147 by the sum of the power of both cases. Then, I (SO,SP ) quantifies the differences with 148 respect to the SHAM case for the frequency bands of interest (see Methods). We characterized the efficiency of different stimulation signals to enhance the 152 occurrence of slow-wave events. In a first step, we tested stimulus pulse features with 153 the open-loop stimulation (STIM-R and STIM-P). We looked at power differences at 154 the frequency bands of interest and the number of occurrences of SWS events while 155 changing the shape, the energy, and the duration of the stimulus pulse. In the following, 156 we will evaluate each of these parameters' impact on the co-occurrence of SP and SOs. 157

158
With STIM-P and STIM-R, the shape, duration, and pulse energy change the temporal 159 and frequency domain characteristics of u(t). The stimulation frequency also modifies 160 them in the case of STIM-P. We tested for the six shapes in Table 2, and shown in  The co-occurrence of SP and SO is evaluated using: (i) the conditional probability 169 of counted co-occurrences given the number of detected spindles, P (C|SP ); and (ii) the 170 probability of slow oscillations, P (SO). In Fig 3E, a positive result means that 171 P (C|SP ) is higher than P (SO); in other words, that the probability of spindles 172 co-occurring with SOs is higher than chance. Given this, both ramp shapes (decreasing 173 and increasing) have positive results for the co-occurrence of events.    always considering that lower duration translates to higher amplitude in order to keep 196 the same pulse energy, which in turn must be enough to overpass the intrinsic activity. 197 Nevertheless, optimal results are found without need to rise over a six-fold value of the 198 intrinsic power per second.

199
When exploring the stimulation frequency using a decreasing ramp of 0.1 seconds, 200 we found no statistical differences between frequencies in the 0.5-1.25 Hz range (see S1 201 SWS events activity variation with closed-loop stimulation 208 Next, we asked whether the phase of the ongoing activity could influence the power and 209 the number of slow-wave sleep events at the moment of the stimulus onset. To answer 210 this, we implemented a closed-loop stimulation protocol. Here, the stimulus pulse is 211 applied when the inverse-notch filtered signal arrives at a particular target phase (see 212 Fig 5A and Methods). Our phase detector used a fixed frequency of the ongoing 213 activity, and 0.85 Hz has been used in other closed-loop implementation as the central 214 frequency of a phase-locked loop (PLL) stage [28,29]. 215 We applied the decreasing ramp stimulation pulses of 0.1 s duration and E p = 40 (a.     Timing of SWS events occurrence 236 A higher occurrence, or higher power of slow-wave events, may not be enough to 237 enhance memory consolidation [30]. The relative timing of occurrence between SPs and 238 SOs also play a role, with SPs being expected to be nested in the up-phase of a slow   Table 3. Kolmogorov-Smirnov test p-values from searching for dissimilarity on the distribution of time-delay from the down-peak of the coincident slow oscillations and spindles. The mark '*' indicates p < 0.05, and '**' indicates p < 0.01. This research asks which stimulation pattern is effective in raising the power of 255 slow-wave sleep events, consequently increasing the SP and SO occurrences. Moreover, 256 we were also interested in maximizing the co-occurrence of SP with the up-peak of SO. 257 Our most significant result is that a pulse of 0.05 seconds with decreasing ramp  The pulse energy must be higher than the intrinsic activity. Following our results, 264 the preferred energy value is 60 (a. u.). However, the translation of this energy value to 265 a physical measurement of sensory stimulation is difficult. Thus, the useful 266 recommendation is to use an energy value around six times the intrinsic EEG activity 267 power per second. The pulse shape and its onset amplitude reshape the SWS events response 274 The shape of the stimulation pulses has small relevance in the previous 275 neuro-stimulation works [32]. One probable cause is the technical difficulty of 276 performing the stimulation following a particular shape-wave with precision in short 277 duration times. For example, the triangular shape in Fig 3A and the Gaussian shape in 278 Fig 3B may be indistinguishable from each other if the stimulation device has a small 279 sampling rate. At the cellular level, the work of Pyragas et al. [33] demonstrates that 280 using Pontryagin's maximum principle (a control theory formulation), a bipolar 281 bang-off-bang waveform is optimal for the entrainment of a neuron to a specific 282 frequency. This waveform type is not suitable for sensory stimulation, because a sensory 283 stimulus affects excitatory and inhibitory neuronal populations at the same time [34]. shapes produce more beneficial and notorious changes for both SWS events. In the 294 same direction, small duration pulses with higher initial amplitude for keeping the same 295 pulse energy showed a higher increase in the power of both SWS events bands and a 296 higher N SP /min. in Fig 4B. Furthermore, the rising ramp (time-mirrored version of the 297 decreasing ramp) results in Fig 3D have the lowest performance of the tested shapes.

298
These results suggest high relevance of the onset amplitude and the intrinsic phase of 299 the stimulus pulse. Secondary effects of these stimuli, such as non-desired arousal, will 300 have to be experimentally determined. Mean-field models of the brain activity relay in two assumptions: The average 304 population dynamics are similar to one neuron dynamics, and the average is enough 305 statistic to characterize the population activity [22,36]. Further, the output activity of a 306 mean-field neural population depends on the dynamics (see Methods Eq (1)) and 307 statistics of the incoming signals (Eq 4). In addition to the u(t) dynamics, the input 308 activity's statistics could play a role on the output dynamics of the model.

309
As in previous works of the Neural Field Theory [23][24][25]27], we used Gaussian white 310 noise as input, which provides the linear impulse response of the model. The noisy 311 input activity represents the net contribution from other sources as the basal ganglia, 312 brainstem, and peripheral nervous system [25]. At the minute scale, variations from the 313 average input value come from these region's activity rather than modifications to the 314 connection's strengths.In our simulations, we used a permanent average value for the 315 noisy input. 316 We superimpose a pulse train in this approach, modifying the input mean value, and 317 its variance. Still, the steady-state results in Fig 1C do not predict all the dynamic 318 effects produced by adding a pulse train at the model's input.

319
The most notable difference from the expected results with the steady-state physiologically-related topography in future works. Consequently, the results reported 335 here are not related to a particular sensory type; they are general in how the increase of 336 activity in a delimited region could modify the whole system's output.

337
More work is needed to consider non-linear effects that surpass the expected 338 response using the linear approximation to explain the obtained results with stimulation. 339 The use of modern control techniques [39], and the integration of information theory Non-invasive methods of neural system stimulation operate at large spatial scales, and 363 they can be by direct application of electrical or magnetic fields or by the use of sensory 364 stimuli to modify neural activity [32,44]. The non-invasive techniques are under 365 active-research because of the potential clinical applications in humans. Several 366 works [12,13] tried to enhance the SWS activity to achieve memory improvements.

367
There were applications of olfactory [45,46], auditory [3,28], and tactile stimuli [47] with 368 similar or even better results than direct stimulation [48][49][50][51][52][53][54]. stimulus at a precise phase of the ongoing cortical activity [3] obtains better behavioral 371 results than rhythmic stimulation [31,55]. Another result of the same authors show that 372 February 1, 2021 12/23 the application of more than two consecutive pulses does not improve physiological 373 outcomes [56]. These results suggest the presence of mechanisms that prevent an 374 over-driving of SO activity, which could be present in our case given a similar increase 375 of I (SO) but lesser N SO /min. with STIM-R (see S1 Fig). A recent work [30] shows that 376 this closed-loop stimulation pattern does not always enhance memory consolidation 377 even with the increase of the power of the SWS events.

378
Our results suggest that the shape, energy, and application phase of the stimulation 379 pulses modifies the power and occurrence of SWS events. On the framework of EEG, 380 these stimuli parameters are often neglected in sensory event-related potential (ERP) 381 analysis, given the confounds with attention in awake experimental studies [57] and the 382 incomplete knowledge about the codification pathways of sensory information. However, 383 these stimulation characteristics are under research for neural implants, and direct 384 stimulation [58]. Finally, the selected stimulation parameters need experimental 385 verification to correlate with memory task performance.

401
Large-scale brain model: Neural Field Theory 402 Neural field theory represents the the activity of a whole population of neurons as its 403 mean firing rate. Spatial propagation can then be calculated to obtain the 404 spatio-temporal evolution of neuronal activity [22].
where α is the inverse of the decaying time of the impulse response, and β is the inverse 416 of the rising time of the impulse response. Both parameters shape the 417 post-synapto-dendritic response. ν ab is the strength of the connection from population b 418 to a, with a, b ∈ b = {e, i, r, s}, and its sign determines whether the connection is 419 excitatory (positive connection strength) or inhibitory (negative connection strength). 420 τ ab is the propagation delay between the populations a and b. Then, the model uses a 421 sigmoid function to convert the soma voltage into the firing rate population response, 422 Q a , as follows where Q max is the maximum firing response, and the sigmoid function approximates the 424 normal cumulative probability distribution of the firing threshold voltage with mean θ 425 and standard deviation σ where σ ρ = σ √ 3/π.

426
The spatial propagation of the firing response only applies to the excitatory 427 population e, the one with axons long enough for making it relevant [23]. It is expressed 428 as where γ e = v e /r e is the temporal damping coefficient for pulses propagating at a 430 velocity v e through axons with a characteristic range r e . For the other populations, 431 φ a (r, t) = Q a (r, t).

432
The expressions Eq (1)-Eq (3) govern the dynamics of each neural population. Their 433 parameters are physiologically related values [23] and are listed in Table 1. Note that 434 the sub-indexes of the parameters follow the order: 'input to'-'output from'.  There are two intra-cortical loops with gains G ee and G ei related to the X-dimension, one intra-thalamic loop with gain G sr G rs related to the Y-dimension, and the gains G es G se and G es G re G sr are cortico-thalamic gain-loops related to the Z-dimension. (B) Diagram of the dynamics of a population of the Neural Field Theory. Particularly, the diagram shows as example the excitatory population e of the eirs-NFT model. (C) Complete simulation diagram including the modules of spatial coupling of the sensory stimulation input and the spatial averaging at the output. The grid in the cortico-thalamic activity block shows the nodes distribution, and the application location of stimuli.

437
To calculate the spectrum of the system's linear approximation with diverse sets of 438 connections strengths, we used the steady-state solutions of the Eq (1) and Eq (3). The 439 steady-state solution of the corticothalamic model comes from setting all temporal and 440 spatial derivatives to zero. Defining φ ν es S ν se φ (0) e + ν sr S ν re φ (0) e + (ν rs /ν es )

444
The linear approximation of the sigmoid function in Eq (2) uses its Taylor's serie where ρ a is the slope of the sigmoid function evaluated at V Taking the Fourier transform of Eq (1), Eq (3), and Eq (6), and rearranging and 451 eliminating φ r and φ s , leads to the system's spectrum with output φ e , and input 452 φ n [24,25].
where G ab = ρ a ν ab is the gain of the connection ab, and L is an abbreviation for L(ω), 454 the reciprocal of the Fourier transform of the operators applied to V a (r, t) in Eq (1) and 455 expressed as The linear approximation in Eq (7)  The eirs-NFT model uses spatio-temporal variables, but our stimuli input is a 465 single-channel time series. To give the stimuli a spatial representation, the signal goes 466 through the sensory input module at the left of Fig 7C. We formed a receptive field 467 using a Difference of Gaussians (DoG) spatial kernel: The spatial filter is centered at an interior node µ r (x = 7, y = 7) of the square sheet of 469 256 nodes with positive standard deviation σ E = 1 and negative standard deviation 470 σ I = 2. The chosen center and standard deviations values avoid the grid boundaries.

471
As shown in Fig 7C, before entering the thalamo-cortical model, we added Gaussian 472 white noise, v(r, t), forming the potential field of the nth population, φ n (r, t).

473
Simulation procedure 474 We solved the model dynamics in Eq (1)-Eq (3) using the Euler's integration method 475 with a time step h = 1e − 4 seconds, and a uniform sheet of N = 256 equidistant nodes 476 in a square cortex of 0.5 meters per side.

477
Each simulation lasts 910 seconds. We performed five simulations for every tested 478 stimulation pattern u(t), each with a different random number generator seed, giving a 479 total of 4550 seconds per stimulation case. The stimulation started after the first 5 480 seconds and stopped before the last 5 seconds. The initial conditions of neural 481 populations variables come from a previous 6-second simulation with initial conditions 482 equal to zero, and without stimulation. 483 We stored the output cortical activity φ e (r, t)), and the input signal φ n (r, t) with a 484 sampling frequency of 100 samples per second. 485 We obtained an EEG-like signal, x(t), spatially averaging the output activity for all 486 time points, and then we subtracted its temporal mean value 487 where r is the grid position, N the number of nodes, andφ e (t, t) the average of the time the signal has a negative peak below -40 µV, and a peak-to-peak amplitude higher 498 to 75 µV. As an alternative, we searched for events with time-lengths between 0.8 and 2 499 seconds (frequency: 0.5 -1.25 Hz) [3], and we kept the higher amplitude oscillations .

500
The detection of SOs is based on the search for high amplitude oscillations in a 501 delta-band sub-region. First, we filtered the signal in the SO frequency band (0.5-1.25 502 Hz, Chebyshev type I filter, fourth-order, 1e-6dB allowed ripple, zero-phase lag).

503
Second, we used the z-score of the filtered signal for searching the zero-crossing points. 504 Later, we applied a threshold to the filtered signal, and we searched for negative peaks 505 below the arbitrary low value of -1e-06 s -1 . Using the zero-crossing points, we calculated 506 the peak-to-peak amplitude for the single oscillations that underpass the negative peak 507 threshold. Finally, slow oscillations are confirmed if the peak-to-peak amplitude is 508 higher than 1.25 times the baseline peak-to-peak amplitude average.

509
Detection of sleep spindles 510 Spindles are bursts of oscillations in the 9-16 Hz band. To detect them, we normalized 511 each simulated output x(t) by the mean and variance of their corresponding b(t) signal. 512 The results are single-channel time series expressed as z-scored EEG registers for each 513 simulation. Then, we computed the root mean square (RMS) value of the z-score x(t) 514 filtered in the SPs frequency band (9-16 Hz, Chebyshev type I filter, fourth-order, 515 1e-6dB allowed ripple, zero-phase lag) and applied an amplitude threshold for the 516 automatic detection of SPs [18,59]. The RMS value calculated at each sample uses a 0.2 517 s window length, and we also smoothed it with a Hamming window of the same time 518 length [60]. SPs are then detected every time the signal overpasses a detection threshold 519 of 1.25 standard deviations of the RMS value of the filtered baseline signal.

520
Counting of co-occurrent events 521 We defined co-occurrence of events every time SOs and SPs overlap at least 250ms. The 522 counting process of co-occurrences uses tagged events to avoid duplication. We also 523 extracted the time percentage of occurrence of SOs, SPs, and co-occurrences. From here, 524 we use the identifier N SO for the number of slow oscillations, N SP for spindles, and N C 525 for co-occurrence of events.

526
Power increases in the event's frequency bands 527 We define the scalogram power difference index (I (j) ), which relates the power of the 528 oscillations in the jth frequency band as Eq (8) 529 where S x (s, t) is the STIM register scalogram, S b (s, t) is the baseline or SHAM register 530 scalogram. s SO comprises the scales inside the 0.5-1.25 Hz spectrum related to slow 531 oscillations and the s SP the range of scales correspondent to the frequency range of the 532 sleep spindles of 9-16 Hz. T is the total simulation time. 533 We used the Morlet wavelet in the scalogram computation. The mother wavelet has 534 central frequency ω 0 = 15, using scaling factors, that result in a frequency resolution of 535 300 scales related to a starting frequency at 0.  The closed-loop stimulation requires the on-line detection of instantaneous phase of the 539 ongoing activity, to apply stimuli at the target phase. For this, we filtered x(t) with an 540 IIR inverse-notch filter for ω 0 = 2πf 0 h with poles z = {0.9999 exp(±jω 0 ), 0}, and zeros 541 z = {−1, −1, 1}, considering a scaling factor of 1e − 4. In addition, we extracted the 542 envelope of the filtered signal at each sample to normalize itself. The envelope detector 543 (see Fig 5A in Results) has three stages: (i) We detected the peaks of the rectified signal 544 using its absolute value. (ii) We exponentially decreased the peak value, with a decay 545 rate of 12.5 s -1 , until the detection of a new peak or the reaching of the limit of a 546 sample counter. (iii) We normalized the filtered signal with the current envelope value 547 at each sample.

555
The stimulus pulse generator is triggered when the phase arrives inside a hysteresis 556 window of ε=±9 degrees (±5%) around the desired target-phase. The trigger waits for 557 the overpass of the phase value of 1.9 or 342 degrees before resetting for the next pulse 558 applying. 559 We applied the stimuli at a specific target-phase of a particular frequency in the 560 closed-loop condition. The signal goes through a narrow band-pass filter and we 561 detected the instantaneous phase with the arc-tangent with an online π/2-shifted phase 562 filtered version of the narrow band signal. We normalized the amplitude of the signal 563 before the arc-tangent by the envelope amplitude. The envelope value came from a peak 564 detector, where the peak value decreases exponentially with a decay rate of 12.5 s -1 565 until the detection of a new peak.

566
Statistical tests 567 We tested for normality of each set of results for each stimulation pattern with the 568 Shapiro-Wilk test.

569
To compare between pairs of stimulation patterns following a normal distribution, 570 we used the Welch's t-test to is the Welch's t-test. On the contrary, when at least one 571 pattern does not follow a normal distribution, we used the Wilcoxon signed-rank test. 572 We used the Kolmogorov-Smirnoff test for the similarity between the distributions of 573 time delays between the down-peak of the slow oscillation and the center of the 574 coincident spindle.

575
All the statistical analysis were implemented using the library scipy.stats 1.5.2 in 576 Python 3.7.